text:
# Fluctuation Cumulant Behavior for the Field-Pulse Induced Magnetisation-Reversal Transition in Ising Models

## I Introduction

The response of a pure magnetic system to time-dependent external magnetic fields has been of current interest in statistical physics mabkc1 ; korn1 ; mish1 . These studies, having close applications in recording and switching industry, have also got considerable practical importance. These spin systems, driven by time-dependent external magnetic fields, have basically got a competition between two time scales: the time-period of the driving field and the relaxation time of the driven system. This gives rise to interesting non-equilibrium phenomena. Tóme and Oliveira first made a mean-field study tome  of kinetic Ising systems under oscillating field. The existence of the dynamic phase transition for such a system and its nature have been thoroughly studied using extensive Monte-Carlo simulations. Later, investigations were extended to the dynamic response of (ferromagnetic) pure Ising systems under magnetic fields of finite-time duration pul1 . All the studies with pulsed field were made below $`T_c^0`$, the static critical temperature (without any field), where the system gets ordered. A ‘positive’ pulse is one which is applied along the direction of prevalent order, while the ‘negative’ one is applied opposite to that. The results for the positive pulse case did not involve any new thermodynamic scale pul1 . In the negative pulse case, however, interesting features were observed pul1 : the negative field pulse competes with the existing order, and the system makes a transition from one ordered state characterised by an equilibrium magnetisation $`+m_0`$ (say) to the other equivalent ordered state with equilibrium magnetisation $`m_0,`$ depending on the temperature $`T`$, field strength $`h_p`$ and its duration $`\mathrm{\Delta }t`$. This transion is well studied in the limit $`\mathrm{\Delta }t\mathrm{}`$ for any non-zero value of $`h_p`$ at any $`T<T_c^0`$. This transition, for the general cases of finite $`\mathrm{\Delta }t`$, is called here ‘magnetisation-reversal’ transition. Some aspects of this transition has been recently studied extensively mish1 ; pul2 .

## II Model and the Transition

The model studied here is the Ising model with nearest-neighbour interaction under a time-dependent external magnetic field. This is described by the Hamiltonian:

$$H=\frac{1}{2}\underset{\{ij\}}{}J_{ij}S_iS_jh(t)\underset{i}{}S_i,$$
(1)
where $`J_{ij}`$ is the cooperative interaction between the spins at site $`i`$ and $`j`$ respectively, and each nearest-neighbour pair is denoted by {…}. We consider a square lattice. The static critical temperature is $`T_c^0=2/\mathrm{ln}(1+\sqrt{2})2.269\mathrm{}`$ (in units of $`J/K_B`$). At $`T<T_c^0`$, an external field pulse is applied, after the system is brought to equilibrium characterised by an equilibrium magnetisation $`m_0(T)`$. The spatially uniform field has a time-dependence as follows:

$$h(t)=\{\begin{array}{cc}h_p& t_0tt_0+\mathrm{\Delta }t\\ 0& \mathrm{otherwise}.\end{array}$$
(2)
Typical time-dependent (response) magnetisation $`m(t)`$ (= $`<S_i>,`$ where $`<\mathrm{}>`$ denotes the thermodynamic ‘ensemble’ average) of the system under different magnetic field $`h(t)`$ are indicated in the Fig. 1. The time $`t_0`$ at which the pulse is applied is chosen such that the system reaches its equilibrium at $`T`$ ($`<`$ $`T_c^0`$ ). As soon as the field is applied, the magnetisation $`m(t)`$ starts decreasing, continues until time $`t+\mathrm{\Delta }t`$ when the field is withdrawn. The system relaxes eventually to one of the two equlibrium states (with magnetisation $`m_0`$ or $`+m_0`$). At a particular temperature $`T`$, for appropriate combinations of $`h_p`$ and $`\mathrm{\Delta }t,`$ a magnetisation-reversal transition occurs, when the magnetisation of the system switches from one state of equilibrium magnetisation $`m_0`$ to the other with magnetisation $`m_0`$. This reversal phenomena at $`T<T_c^0`$ is simple and well studied for $`\mathrm{\Delta }t\mathrm{}`$ for any non-zero $`h_p`$. We study here the dynamics for finite $`\mathrm{\Delta }t`$ values. It appears that generally $`h_p\mathrm{}`$ as $`\mathrm{\Delta }t0`$ and $`h_p0`$ as $`\mathrm{\Delta }t\mathrm{}`$ for any such dynamic magnetisation-reversal transition phase boundary at any temperature $`T`$ ($`<`$ $`T_c^0`$). In fact, a simple application of the domain nucleation theory gives $`h_p\mathrm{ln}\mathrm{\Delta }\mathrm{t}`$ = constant along the phase boundary, where the constant changes by a factor $`1/(d+1)`$, where $`d`$ denotes the lattice dimension, as the boundary changes from single to multi-domain region pul1 .

FIG. 1. Typical time variation of the response magnetisation $`m(t)`$ for two different field pulses $`h(t)`$ with same $`\mathrm{\Delta }t`$ for an Ising system at a fixed temperature $`T`$. The magnetisation-reversal here occurs due to increased pulse strength, keeping their width $`\mathrm{\Delta }t`$ same. The transition can also be brought about by increasing $`\mathrm{\Delta }t`$, keeping $`h_p`$ fixed. The inset indicates the typical phase boundaries (where the field withdrawal-time magnetisation $`m_w=0`$) for two different temperatures (sequential updating; note that for random updating the phase boundaries shift upwards).

A mean field study of the problem gives a qualitative understanding of the diverging time and length scales developed near the transition boundary (in the $`h_p`$ \- $`\mathrm{\Delta }t`$ plane at a fixed $`T<T_c^0`$). Mean-field approximation for the dynamics gives the equation of motion for the average magnetisation $`m_i`$ as

$$\frac{dm_i}{dt}=m_i+\mathrm{tanh}\left(\frac{_jJ_{ij}m_j+h(t)}{T}\right).$$
(3)
This equation, linearised near the magnetisation-reversal transition point, gives, for uniform magnetisation,

$$m(t)=\frac{h_p}{\mathrm{\Delta }T}\left(\frac{h_p}{\mathrm{\Delta }T}m_0\right)\left[\mathrm{exp}\left\{\frac{\mathrm{\Delta }T}{T}(tt_0)\right\}\right]$$
(4)
as a solution of eqn. (3), for $`t_0tt_0+\mathrm{\Delta }t`$. Here $`\mathrm{\Delta }T=T_c^{mf}T`$, where $`T_c^{mf}J(q=0)`$ is the static critical temperature in the mean-field approximation and $`J(q)`$ is the Fourier transform of the interaction $`J_{ij}`$. Due to application of the field $`h_p,`$ $`m(t)`$ decreases in magnitude from $`m(t_0)m_0`$ to $`m(t_0+\mathrm{\Delta }t)m_w`$ at the time of withdrawal of the pulse. Due to absence of fluctuation here, magnetisation relaxes back to its original value $`m_0`$ if $`m_w`$ is positive, or to a value $`m_0`$ if $`m_w`$ is negative. In the $`t>t_0+\mathrm{\Delta }t`$ regime, where $`h(t)=0,`$ the magnetisation (starting from $`m_w`$ at $`t=t+\mathrm{\Delta }t`$) relaxes back to its final equilibrium value $`\pm m_0`$, with a relaxation time mish1 ; pul1

$$\tau \frac{1}{(T_c^{mf}T)}\mathrm{ln}\left(\left|\frac{m_0}{m_w}\right|\right).$$
(5)
It diverges at the magnetisation-reversal transition boundary, where $`m_w`$ vanishes. The prefactor gives the divergence of $`\tau `$ at the static mean field transition temperature, and is responsible for critical slowing down phenomena at the static transition point ($`h=0`$). The other factor gives the diverging time scale, at any temperature below the static transition temperature, where magnetisation reversal occurs or $`m_w`$ vanishes due to appropriate combination of $`h_p`$ and $`\mathrm{\Delta }t.`$ The solution of the susceptibility $`\chi (q)`$ gives mish1

$$\chi (q)\mathrm{exp}\left(q^2\xi ^2\right),$$
(6)
where the correlation length is given by

$$\xi \left[\frac{1}{(T_c^{mf}T)}\mathrm{ln}\left(\left|\frac{m_0}{m_w}\right|\right)\right]^{\frac{1}{2}}.$$
(7)
Here, too, the prefactor in $`\chi `$ gives the usual divergence at $`T_c^{mf},`$ while the other factor gives the divergence at the magnetisation-reversal transition point. Incorporating fluctuations, extensive Monte-Carlo simulation studies have also convincingly demonstrated pul2  that the fluctuation in the order parameter $`|m_w|`$ and in the internal energy of the system grows with the system size and diverges at the magnetisation-reversal transition boundary, where $`m_w`$ vanishes.

## III Monte-Carlo study and the results

Here the Monte-Carlo study has been carried out in two-dimensions (square lattice) with periodic boundary conditions. Spins are updated following Glauber dynamics. The updating rule employed here are both sequential as well as random. In sequential updating rule one Monte-Carlo step consist of a complete scan of the lattice in a sequential manner; while in random updating a Monte-Carlo step is defined by $`N`$ (= $`L^2`$ ) random updates on the lattice, where $`N`$ is the total number of spins in a lattice of linear size $`L`$. Studies have been carried out at temperatures below the static critical temperature ($`T_c^02.27`$). The system is allowed to evolve from an initial state of perfect order to its equilibrium state at temperature $`T`$. The time $`t_0`$ is chosen to be much larger than the static relaxation time at that $`T`$, so that the system reaches an equilibrium state with magnetization $`+m_0(T)`$ before the external magnetic field is applied at time $`t=t_0`$. The field pulse of strength $`h_p`$ is applied for duration $`\mathrm{\Delta }t`$ (measured in Monte Carlo steps or MCS). The magnetisation starts decreasing from its equilibrium value $`m_0`$. The average value of the magnetisation $`m_w`$ at the time of withdrawal of pulse is noted. The phase boundary of this dynamic transition is defined by appropriate combination of $`h_p`$ and $`\mathrm{\Delta }t`$ that produces the magnetisation reversal by making $`m(t_0+\mathrm{\Delta }t)m_w=0`$ from a value $`m(t_0)=m_0,`$ i.e, $`m_w`$ changes sign across the phase boundary. The phase boundary changes with $`T.`$ The bahavior of different thermodynamic quantities are studied across the phase boundary. These quantities are averaged over $`100020000`$ different initial configurations of the system. The fluctuations over the average value are also noted.

Here we study the behavior of the reduced fourth order cumulant $`U`$ binheer  near the magnetisation reversal transition. This is defined as

$$U=1\frac{m_w^4}{3m_w^2^2},$$
(8)
where $`m_w^4`$ is the ensemble average of $`m_w^4`$. $`m_w^2`$ is similarly defined. The cumulant $`U`$ here behaves somewhat differently, compared to that in static and other transitions: Deep inside the ordered phase $`m_w1`$ and $`U2/3.`$ For other (say, static) transitions the order parameter ($`m_w`$) goes to zero with a Gaussian fluctuation above the transition point, giving $`U0`$ there. Here, however, due to the presence of the pulsed field, $`\left|m_w\right|`$ is non-zero on both sides of the magnetisation-reversal transition. Hence $`U`$ drops to zero at a point near the transition and grows again after it.

The universality class of the dynamic transition in Ising model under oscillating field has been studied extensively by investigating korn1  the critical point and the cumulant value $`U^{}`$ at the critical point, where the cumulant curves cross for different system sizes ($`L`$). In that case, of course, the variation of $`U`$ (at any fixed $`L`$) is similar to that in the static Ising transitions ($`U=2/3`$ well inside the ordered phase and $`U0`$ well within the disordered phase). In fact, $`U^{}`$ value in this oscillatory field case was found to be the same as that in the static case, indicating the same universality class korn1 . We observe different bahaviour in the field pulse induced magnetisation-reversal transition case.

We observe two kinds of distinct behavior of the cumulant $`U.`$ Typically, for low temperature and low pulse-duration region (see the inset in Fig. 1) of the magnetisation-reversal phase boundary, the cumulant crossing for different system sizes ($`L`$) occur at $`U^{}0.42`$ to $`0.46`$ (see Fig. 2). As mentioned already, we have checked these results for both sequential and random updating. Specifically, for $`T=0.5`$ and $`\mathrm{\Delta }t=5,`$ (see Fig. 2c) we find the transition point value of $`h_p2.6,`$ to be smaller than the value ($`1.9`$) for sequential updating. However, the value of $`U^{}`$ at this transition point is again very close to about $`0.44.`$ This indicates that updating rule does not affect the universality class ($`U^{}`$ value), as long as the proper region of the phase boundary is considered. For relatively higher temperature and pulse-duration region of the phase boundary, the crossing of $`U`$ for different $`L`$ values occur for $`U^{}0_+`$. This is true for both sequential (Fig. 3a, b) and random (Fig. 3c) updating. It may be noted that the phase boundary changes with the updating rule, as the system relaxation time (which matches with the pulse width at the phase boundary) is different for sequential and random updating binheer .

FIG. 2. Behavior of $`U`$ near the transition, driven by (a) $`T`$ at a fixed value of $`h_p`$ (=1.9) and $`\mathrm{\Delta }t`$ (=$`5`$) with sequential updating, (b) $`h_p`$ at a fixed value of $`T`$ (=0.5) and $`\mathrm{\Delta }t`$ (=$`5`$) with sequential updating, and (c) $`h_p`$ at a fixed value of $`T`$ (=0.5) and $`\mathrm{\Delta }t`$ (=$`5`$) with random updating, for different $`L`$, averaged over 1000 to 20000 initial configurations. The fluctuations are smaller than the symbol size. The insets show the typical behavior of the magnetisation $`m_w`$ at the time of withdrawal of the field pulse by varying (a) $`T`$ at a fixed $`h_p`$ and $`\mathrm{\Delta }t`$, for $`L=100`$ and $`800`$, (b) $`h_p`$ at a fixed $`T`$ and $`\mathrm{\Delta }t`$, for $`L=100`$ and $`400`$, (c) $`h_p`$ at a fixed $`T`$ and $`\mathrm{\Delta }t`$, for $`L=50`$ and $`200`$; $`m_w=0`$ at the effective transition point. (d) Finite size scaling study in this parameter range: the effective $`T_c`$ or $`h_p^c`$ values (see the insets), where $`m_w=0`$, are plotted against $`L^{1/\nu }`$ with $`\nu ^1=0.7`$. The values of the cumulant crossing points in (a), (b), (c) are taken to correspond the respective transition points for $`L\mathrm{}`$.

It might be noted that in the low temperature and $`\mathrm{\Delta }t`$ regions, there seems to be significant finite size scaling of the transition ($`m_w=0)`$ point (see the insets of Fig. 2a, b, c). In fact, in Fig. 2d, the finite-size scaling analysis of those data is presented. For the other cases, there seems to be no significant finite size effect on the transition point (cf. insets of Fig. 3a, b, c), indicative of a mean-field nature of the transition in this range. It may be noted that to compare the finite size effects, we normalise the parameters $`T`$ or $`h_p`$ by their ranges required for full magnetisation reversal. In fact, this weak finite size effect for high $`T`$ and $`\mathrm{\Delta }t`$ regions did not lead to any reasonable value for the fitting exponent in the scaling analysis.

FIG. 3. Behavior of $`U`$ near the transition, driven by (a) $`h_p`$ at a fixed value of $`T`$ (=2.0) and $`\mathrm{\Delta }t`$ (=$`5`$) with sequential updating, (b) $`T`$ at a fixed value of $`h_p`$ (=0.5) and $`\mathrm{\Delta }t`$ (=$`10`$) with sequential updating, and (c) $`h_p`$ at a fixed value of $`T`$ (=1.5) and $`\mathrm{\Delta }t`$ (=$`5`$) with random updating, for different $`L`$, averaged over 1000 to 6000 initial configurations. The fluctuations are smaller than the symbol size. The insets show the typical behavior of the magnetisation $`m_w`$ at the time of withdrawal of the field pulse by varying (a) $`h_p`$ at a fixed $`T`$ and $`\mathrm{\Delta }t`$, for $`L=50`$ and $`400`$, (b) $`T`$ at a fixed $`h_p`$ and $`\mathrm{\Delta }t`$, for $`L=50`$ and $`200`$, (c) $`h_p`$ at a fixed $`T`$ and $`\mathrm{\Delta }t`$, for $`L=50`$ and $`200`$; $`m_w=0`$ at the transition point.

For the static transition of the pure two-dimensional Ising system, $`U^{}0.6107`$ binheer ; blote ; bincu . For low temperature (and low $`\mathrm{\Delta }t`$) regions of the magnetisation-reversal phase boundary, the observed values of $`U^{}`$ (in the range 0.42 - 0.46) are considerably lower than the above mentioned value for the static transition. There is not enough indication of finite-size effect in the $`U^{}`$ value either (cf. korn1 ). This suggests a new universality class in this range. Also, the finite-size scaling study for the effective transition points here (see Fig. 2d) gives a correlation length exponent value ($`\nu 1.4`$) larger than that of the static transition. For comparatively higher temperatures (and high $`\mathrm{\Delta }t`$), the $`U^{}0_+`$ at the crossing point. Such small value of the cumulant at the crossing point can hardly be imagined to be a finite-size effect; it seems unlikely that one would get here also the same universality class and $`U^{}`$ value will eventually shoot up to $`U^{}0.44`$ (for larger system sizes), as for the other range of the transition. On the other hand, such low value of $`U^{}`$ might indicate a very weak singularity, as indicated by the mean field calculations mish1  mentioned in the introduction. In fact, even for the static transition, as the dimensionality increases, and the singularity becomes weaker (converging to mean field exponents) with increasing lattice dimension, the cumulant crossing point $`U^{}`$ decreases ($`U^{}0.61`$ in $`d=2`$ to $`U^{}0.44`$ in $`d=4`$) bincu . We believe the mean field transition behavior here, as mentioned earlier, is even weaker in this dynamic case as reflected by the value $`U^{}0_+`$, corresponding to a logarithmic singularity (as in eqns. (5) and (7)).

## IV Summary and conclusions

The universality class of the dynamic magnetisation-reversal transition, induced by a competing pulse, in an Ising model on a square lattice, below its static ordering temperature, is studied here using Monte Carlo simulations. Both sequential and random updating have been used. The phase boundary at any $`T`$ ($`<T_c^0`$) is obtained first in the $`h_p\mathrm{\Delta }t`$ plane. They of course depend on the updating rule. The phase boundaries obtained compare well with the nucleation theory estimate $`h_p\mathrm{ln}\mathrm{\Delta }\mathrm{t}`$ = constant along the boundary pul1 . The mean-field theory applications mish1 ; pul1  indicated time and length (eqns. 5 and 7 respectively) scale divergences at these phase boundaries. Extensive Monte-Carlo studies for the fluctuations in the order parameter $`|m_w|`$ and internal energies etc. showed prominent divergences along the phase boundaries pul2 . Fourth order cumulant ($`U`$) of the order parameter distribution is studied here for different system sizes (upto $`L=800`$) around the phase boundary region. The crossing point of the cumulant (for different system sizes) gives the transition point and the value $`U^{}`$ of the cumulant at the transition point indicates the universality class of the transition. In the low temperature and low pulse width range, the $`U^{}`$ value is found to be around 0.44 (see Figs. 2a, b, c). The prominent discripancy with the $`U^{}`$ value ($`0.61`$) for the static transition in the same model in two dimensions indicates a new universality class for this dynamic transition. Indeed, the finite-size scaling analysis (Fig. 2d) suggests a different (larger) value of the correlation length exponent also. For comparatively higher temperatures and higher pulse widths, the $`U^{}`$ values are very close to zero (see Fig. 3a, b, c), and the transitions here seem to fall in a mean-field-like weak-singularity universality class similar to that obtained earlier mish1 , and indicated by eqns. (5) and (7). Here, the finite size effects in the order parameter and the transition point are also observed to be comparatively weaker (see insets of Fig. 3).

## V Acknowledgements

The authors are grateful to Bernard Nienhuis for some of his insightful questions regarding the transition, and to Arkajyoti Misra for useful comments and discussions.


format: text

text:
# A Novel Path Sampling Method for the Calculation of Rate Constants

## I INTRODUCTION

The calculation of rate constants of activated processes dominated by rare events, chemical reactions being a prime example, is still one of the major computational challenges. As transition rates depend exponentially on the activation barrier height, the expectation time for an event can exceed current computer capabilities by many orders of magnitude. As a result most chemical reactions can not be simulated by direct molecular dynamics (MD) methods, except those with very low activation energies. The conventional way to tackle this time scale problem is based on transition state theory (TST) and separates the problem in two steps Keck62 ; Anderson73 ; Bennet77 ; DC78 . The first step is the calculation of the free energy barrier as function of a reaction coordinate, the second stage is the calculation of the transmission coefficient by sampling fleeting trajectories departing from the top of the barrier. If the reaction coordinate is well chosen, the top of the free energy barrier corresponds to points in phase space close to the true transition state, and the transmission coefficient will have a reasonable value. However, in high dimensional complex systems the choice of reaction coordinate can be extremely difficult and usually requires detailed a priori knowledge of the transition mechanism. Consequently, an intuitively chosen but wrong reaction coordinate can result in a very low transmission coefficient, and hence a statistically inaccurate or immeasurable rate constant.

Chandler and collaborators TPS98 ; TPS98a ; TPS99 ; Bolhuis02 ; Dellago02  devised a method for which no prior knowledge of the system is needed. This method, called transition path sampling (TPS), gathers a collection of trajectories connecting the reactant to the product region by employing a Monte Carlo (MC) algorithm. The resulting path ensemble can be used to elucidate reaction mechanisms, transition states and reaction coordinates. The TPS method has been successfully used on such diverse systems as cluster isomerization, auto-dissociation of water, ion pair dissociation and on isomerization of a dipeptide, as well a reactions in aqueous solution (see Ref. Bolhuis02  for an overview). Just as in the conventional case mentioned above, an additional second simulation is needed to determine the rate constant within TPS. This simulation combines the path sampling method with the umbrella sampling technique to estimate the probability to reach the product state from the initial reactant state. The final macroscopic rate constant is given by a plateau in the time derivative of a correlation function TPS99 . In case of two distinct stable states this plateau region should always exist at times longer than the typical molecular relaxation time. However, when reaction pathways are complex and exhibit multiple recrossings, these typical molecular relaxation timescales can be relatively long. In that case the TPS rate constant calculation is computationally expensive, as the path length must exceed these timescales.

In this paper we improve the efficiency of the TPS rate constant calculation on several points by introducing an alternative scheme for calculating reaction rates, named transition interface sampling (TIS). The first of these improvements is allowing the path length to vary, so that by a well chosen definition of the stable states we can limit the length of each path to the strict minimum. Secondly, the new method is based on the effective positive flux through dividing surfaces or interfaces and is consequently much less sensitive to multiple recrossings or diffusive barrier crossings. Thirdly, the number of different types of Monte Carlo moves is reduced, making the implementation of the algorithm conceptually simpler.

This paper is organized as follows: In Sec. II we briefly describe the existing algorithms and present the theoretical derivation for the TIS rate constant expression. The implementation of the algorithm is discussed in Sec. III. We illustrate the algorithm on a diatomic molecule in a fluid of repulsive particles and make a quantitative comparison to the original TPS calculation in Sec. IV. We end with concluding remarks in Sec. V.

## II THEORY

### II.1 Transition state theory and the calculation of rate constants

Consider a dynamical system in which transitions can take place between two stable states $`A`$ and $`B`$. If the barrier between $`A`$ and $`B`$ is sufficiently high, the system will show exponential relaxation for which the forward and backward rate constants $`k_{AB}`$ and $`k_{BA}`$ are well defined and can be expressed in terms of microscopic properties. Measuring these rate constants by computer simulation is traditionally done by the two stage Bennett-Chandler (BC) procedure based on the principles of TST Bennet77 ; DC78 . The first step is the calculation of the reversible work or free energy to bring the system from stable state $`A`$ to the transition state. This free energy $`F(\lambda )`$ has to be calculated as a function of a suitably chosen reaction coordinate $`\lambda `$. This $`\lambda `$ can be a complex function of all particle coordinates $`r`$ and momenta $`p`$: $`\lambda =\lambda (x)`$, with $`x\{r,p\}`$. The maximum in $`F(\lambda )`$ defines the transition state dividing surface $`\lambda ^{}`$ FrenkelSmit ; Chandlerbook . By convention, the system is in $`A`$ if $`\lambda (x)<\lambda ^{}`$ and in $`B`$ if $`\lambda (x)>\lambda ^{}`$.

The main assumption in TST is that any trajectory coming from $`A`$ and crossing the transition state dividing surface $`\lambda (x)=\lambda ^{}`$ will remain at the $`B`$ side of the dividing surface for a long time. The reaction rate can therefore be expressed as the positive flux through the multidimensional dividing surface $`\lambda ^{}`$

$`k_{AB}^{TST}`$ $`=`$ $`\underset{\mathrm{\Delta }t0}{lim}{\displaystyle \frac{1}{\mathrm{\Delta }t}}{\displaystyle \frac{\theta \left(\lambda ^{}\lambda (x_0)\right)\theta \left(\lambda (x_{\mathrm{\Delta }t})\lambda ^{}\right)}{\theta \left(\lambda ^{}\lambda \right)}}`$ (1)
$`=`$ $`{\displaystyle \frac{\dot{\lambda }(x_0)\delta \left(\lambda (x_0)\lambda ^{}\right)\theta \left(\dot{\lambda }(x_0)\right)}{\theta \left(\lambda ^{}\lambda \right)}}`$
$`=`$ $`\dot{\lambda }(x_0)\theta \left(\dot{\lambda }(x_0)\right)_\lambda ^{}{\displaystyle \frac{e^{\beta F(\lambda ^{})}}{_{\mathrm{}}^\lambda ^{}e^{\beta F(\lambda )}𝑑\lambda }},`$
where $`x_t`$ specifies the set of coordinates and momenta of the system at time $`t`$, the dots denote derivatives with respect to time $`t`$, the brackets $`\mathrm{}`$ denote equilibrium ensemble averages and $`\theta (x)`$ and $`\delta (x)`$ are the Heaviside step-function and the Dirac delta function respectively. In the last equality of Eq. (1) the connection to the reversible work $`F(\lambda )`$ is made, and $`\beta =1/k_\mathrm{B}T`$, where $`k_\mathrm{B}`$ is Boltzmann’s constant and $`T`$ is the temperature. The subscript $`\lambda ^{}`$ to the ensemble brackets, indicates that the ensemble is constrained to the top of the barrier $`\lambda ^{}`$.

We consider the system to be completely deterministic and thus we can write $`x_t=f(x_t^{},tt^{})=f(x_0,t)`$, in which $`f`$ is the time-propagator function. Evaluation of the function $`f(x,t)`$ requires integrating the equations of motion over the time interval $`t`$ starting with configuration $`x`$. Nevertheless, the equations derived in this paper are still valid when applied to stochastic dynamics.

Even when the TST assumption is accurate, it can be extremely difficult to find a proper reaction coordinate for which recrossings do not occur. As a result a wrong choice for the reaction coordinate will give a much lower free energy barrier than the real activation free energy and will correspondingly overestimate the rate constant. Figure 1 illustrates that Eq. (1) overcounts trajectories. One can correct for this overcounting by multiplying the TST rate constant with the transmission coefficient $`\kappa (t)`$ to obtain the true rate constant

$$k_{AB}(t)=k_{AB}^{TST}\kappa (t).$$
(2)
The calculation of the time dependent transmission coefficient $`\kappa (t)`$ constitutes the second part of the two stage BC procedure Bennet77 ; DC78 . $`\kappa (t)`$ belongs to the approximate dividing surface $`\lambda ^{}`$ Keck62 ; Anderson73 ; Bennet77 ; DC78  and can be determined by taking an ensemble average of many short trajectories starting from the dividing surface:

$`\kappa (t)`$ $`=`$ $`{\displaystyle \frac{1}{\dot{\lambda }(x_0)\theta \left(\dot{\lambda }(x_0)\right)_\lambda ^{}}}\dot{\lambda }(x_0)\theta \left(\lambda (x_t)\lambda ^{}\right)_\lambda ^{}.`$ (3)
After a short molecular time $`t_{\mathrm{mol}}`$ the trajectories are committed to a stable state and $`\kappa (t)`$, and hence $`k_{AB}(t)`$, become constant: the transmission coefficient $`\kappa `$, and the rate constant $`k_{AB}`$, respectively. It is however important to start sufficiently close to the true transition state dividing surface. Otherwise the transmission coefficient will be extremely low, making an accurate estimate of the rate constant problematic or even impossible. In many cases, in particular for complex condensed matter systems, a sufficiently close reaction coordinate is difficult to find and requires considerable a-priori knowledge about the system.

### II.2 Transition path sampling

Transition path sampling (TPS) is developed to overcome the difficulties mentioned above TPS98 ; TPS98a ; TPS99 ; Bolhuis02 ; Dellago02 . Its main advantage is that no prior knowledge of the transition state is needed. The rate constant in TPS is expressed as the time derivative of a general time correlation function.

$$k_{AB}^{TPS}(t)=\frac{\mathrm{d}}{\mathrm{d}t}C(t),C(t)=\frac{h_A(x_0)h_B(x_t)}{h_A(x_0)},$$
(4)
in which $`h_A(x)`$ and $`h_B(x)`$ are the characteristic functions defined by:

$`h_A(x)`$ $`=`$ $`1,\text{ if }xA,\text{ else }h_A(x)=0`$
$`h_B(x)`$ $`=`$ $`1,\text{ if }xB,\text{ else }h_B(x)=0.`$ (5)
In case of a single order parameter $`\lambda (x_t)`$ describing the transition, the phase space regions $`A`$ and $`B`$ are defined by $`\lambda _A`$ and $`\lambda _B`$: $`x_tA\text{ if }\lambda (x_t)<\lambda _A`$ and $`x_tB\text{ if }\lambda (x_t)>\lambda _B`$. Knowledge of the precise location of the dividing surface $`\lambda ^{}`$, $`\lambda _A<\lambda ^{}<\lambda _B`$, is not required in TPS. Therefore, the order parameter $`\lambda `$ does usually not correspond to the reaction coordinate.

The microscopic expression for the rate constant in Eq. (4) is time dependent, while the phenomenological rate constant is not. However, just as the transmission coefficient $`\kappa (t)`$ becomes a constant, the time dependent function $`k_{AB}^{TPS}(t)`$ reaches a plateau after a molecular timescale $`t_{\mathrm{mol}}`$. The phenomenological rate constant is equal to the plateau value: $`k_{AB}=k_{AB}^{TPS}(T)`$. This plateau region should always exist for times $`T`$ between the molecular timescale and the characteristic reaction time: $`t_{\mathrm{mol}}<Tt_{\mathrm{rxn}}`$. In other words, $`T`$ is larger than the timescale to commit to one of the stable states, but much shorter than the expectation time $`t_{\mathrm{rxn}}`$ of a completely new reactive event. If we take $`\lambda _A=\lambda _B=\lambda ^{}`$ and the limit $`t0+`$, Eq. (4) transforms into the expression for the positive reactive flux or, equivalently, the TST equation (1). For $`t>0`$, however, the reactive flux measured by Eq. (4) no longer consists of purely positive contributions. The final rate constant is a sum of positive and negative fluxes, and thus the overcounting of trajectories in Eq. (1) is circumvented. (See fig.(1)).

We can rewrite the time dependent rate constant of Eq. (4) into TPS99 :

$$k_{AB}^{TPS}(t)=\frac{\dot{h}_B(t)_{A,H_B(T)}}{h_B(t^{})_{A,H_B(T)}}C(t^{}),$$
(6)
where $`H_B(T)=\mathrm{max}_{0<t<T}h_B(x_t)`$ and $`\mathrm{}_{A,H_B(T)}`$ denotes an average on the ensemble of paths of fixed length $`T`$ starting in $`A`$ and entering $`B`$ at least once TPS99 . These ensemble averages are evaluated using a Monte Carlo procedure employing the *shooting* and *shifting* moves TPS98a . The two factors in Eq. (6) have to be evaluated separately. First, a path sampling simulation is performed to compute $`h_B(t)_{A,H_B(T)}`$ in the interval $`[0,T]`$. The path length $`T`$ must be long enough for the time derivative to display a plateau. Subsequently, one chooses a $`t^{}`$ in interval $`[0,T]`$ and computes $`C(t^{})`$ using the path sampling in combination with an umbrella sampling technique TPS99 . A drawback of the TPS rate constant calculation is that the function $`k_{AB}^{TPS}(t)`$ can be strongly oscillatory because of recrossings and will reach a plateau only after a relatively long time. The path length in TPS must exceed the typical timescale of these oscillations, and consequently, in that case TPS is computationally costly.

### II.3 Transition Interface Sampling

Just as the BC and the TPS rate constant algorithms, the TIS method is based on a flux calculation. In contrast to these schemes, however, TIS measures the effective positive flux 111Here, effective means that the recrossings through the interfaces are not being counted, instead of a conditional general flux as in Eq. (4) or Eq. (1). This implies that only positive terms contribute to the rate, allowing for faster numerical convergence. A flux is normally defined through a hypersurface in phase space defined by an order parameter, the reaction coordinate. But, similar to the TPS case, we do not want to suffer from a bad choice of reaction coordinate. Therefore, instead of using a single dividing surface, we introduce a series of interfaces through which we measure this flux. We then derive an expression that relates the flux through a certain interface to the flux through an interface which is closer to $`A`$ to replace the expensive TPS umbrella sampling procedure.

In order to formulate a proper flux, we have to divide the entire phase space into two complementary regions called overall states $`𝒜`$ and $``$. These states do not only depend on the position at the time of consideration but also on its past behavior. Overall state $`𝒜`$ covers all phase space points lying inside stable region $`A`$, which constitutes the largest part, but also all phase space points that visit $`A`$, before reaching $`B`$ when the equations of motion are integrated backward in time. Similarly, state $``$ comprises stable state $`B`$ and all phase points, coming directly from this state in the past, i.e. without having been in $`A`$. It is useful to generalize the characteristic functions in Eq.(5) for an arbitrary phase space region $`\mathrm{\Omega }`$

$`h_\mathrm{\Omega }(x)`$ $`=`$ $`1,\text{ if }x\mathrm{\Omega },\text{ else }h_\mathrm{\Omega }(x)=0.`$ (7)
For each phase point $`x`$ and each phase space region $`\mathrm{\Omega }`$ we can determine the minimum (first entrance) times $`t_\mathrm{\Omega }^b(x)`$ and $`t_\mathrm{\Omega }^f(x)`$ needed to reach $`\mathrm{\Omega }`$ starting from configuration $`x`$ by integrating the equations of motion backward and forward in time, respectively:

$`t_\mathrm{\Omega }^b(x)`$ $``$ $`\mathrm{max}\left[\{t|h_\mathrm{\Omega }(f(x,t))=1t0\}\right]`$
$`t_\mathrm{\Omega }^f(x)`$ $``$ $`+\mathrm{min}\left[\{t|h_\mathrm{\Omega }(f(x,t))=1t0\}\right],`$ (8)
where the $`\mathrm{min}`$ and $`\mathrm{max}`$ function return respectively the lowest and highest value of their arguments. In addition, it is useful to define for each phase point $`x`$ and each set of two non-overlapping phase space regions $`\{\mathrm{\Omega }_1,\mathrm{\Omega }_2\}`$ the following characteristic functions:

$`\overline{h}_{\mathrm{\Omega }_1,\mathrm{\Omega }_2}^b(x)`$ $`=`$ $`\{\begin{array}{cc}1\text{ if }h_{\mathrm{\Omega }_1}\left(f(x,t_{\mathrm{\Omega }_1\mathrm{\Omega }_2}^b(x))\right)=1,\hfill & \\ 0\text{ otherwise}\hfill & \end{array}`$
$`\overline{h}_{\mathrm{\Omega }_1,\mathrm{\Omega }_2}^f(x)`$ $`=`$ $`\{\begin{array}{cc}1\text{ if }h_{\mathrm{\Omega }_1}\left(f(x,+t_{\mathrm{\Omega }_1\mathrm{\Omega }_2}^f(x))\right)=1,\hfill & \\ 0\text{ otherwise}\hfill & \end{array}`$ (9)
In words, these functions measure whether a trajectory reaches $`\mathrm{\Omega }_1`$ before $`\mathrm{\Omega }_2`$ or not. As the system is ergodic, each phase space region will be visited in finite time and thus $`\overline{h}_{\mathrm{\Omega }_1,\mathrm{\Omega }_2}^b(x)+\overline{h}_{\mathrm{\Omega }_2,\mathrm{\Omega }_1}^b(x)=\overline{h}_{\mathrm{\Omega }_1,\mathrm{\Omega }_2}^f(x)+\overline{h}_{\mathrm{\Omega }_2,\mathrm{\Omega }_1}^f(x)=1`$ for any $`x`$. Using these definitions the characteristic functions for the overall states $`𝒜`$ and $``$ are given by

$`h_𝒜(x)=\overline{h}_{A,B}^b(x),h_{}(x)=\overline{h}_{B,A}^b(x).`$ (10)
These states together span the complete phase space, as the system can never stay in the intermediate region between $`A`$ and $`B`$ forever. The overall states $`𝒜`$ and $``$ do not sensitively depend on the definition of stable state $`A`$ and $`B`$ as long as it is reasonably. Of course, the stable regions should not overlap, each trajectory between the stable states must be a true rare event for the reaction we are interested in. In addition, the probability that after this event the reverse reaction occurs shortly thereafter must be as unlikely as an entirely new event. In other words, the system must be committed to the stable states. Therefore, a reasonable definition of $`A`$ and $`B`$ requires that they should lie completely inside the basin of attraction of the respective two states 222In general, $`A`$ and $`B`$ are defined in phase space, but for most practical cases configuration space might be enough. (see also Ref. Dellago02 ). Special care has to be taken with this condition for processes which show many recrossings between state $`A`$ and $`B`$ before settling down. Such processes can occur in solution or in dilute gasses. For instance, for organic reactions in aqueous solution, a rare specific hydrogen bonded network can lower the bond-breaking barrier and initiate the reaction. If the lifetime of those rare solvation structures is high, a sudden reverse reaction can occur as the barrier for the backward reaction is also lowered by the same amount Strnad ; Okuno ; Erp . A similar phenomenon can happen in dilute gasses for which rare spontaneous fluctuations in the kinetic energy are the main driving force. A particle moving from one state to another due to a very high kinetic energy as result of sequence of collisions can cross the potential energy barrier several times before it will dissipate its energy by a new collision and relax into one of the stable states (see e.g. Refs. CHO ; Sun ). These problems can in principle be solved by an adequate choice of the stable state definitions. For instance, the definition can depend explicitly on the presence of certain hydrogen bonds or on kinetic energy terms.

With our definition of overall states $`𝒜`$ and $``$ we can write down our rate equation in the spirit of Eq. (4):

$$k_{AB}=\frac{h_𝒜(x_0)\dot{h}_{}(x_0)}{h_𝒜(x_0)},$$
(11)
where the dot denotes the time derivative taken at $`t=0`$. This rate expression does not depend on time although the evaluation of the characteristic functions still requires integration of the equations of motion. The transition from $`𝒜`$ into $``$ takes place when the system coming from $`A`$ will cross the interface $`\lambda _B`$ for the first time (see Fig. 2). After this event the system will stay in $``$. Eq. (11) counts therefore only the first crossing through interface $`\lambda _B`$ and is hence equivalent to the effective positive flux expression

$`k_{AB}`$ $`=`$ $`{\displaystyle \frac{h_𝒜(x_0)\dot{\lambda }(x_0)\theta (\dot{\lambda }(x_0))\delta (\lambda (x_0)\lambda _B)}{h_𝒜(x_0)}}`$
$`=`$ $`\underset{\mathrm{\Delta }t0}{lim}{\displaystyle \frac{1}{\mathrm{\Delta }t}}{\displaystyle \frac{h_𝒜(x_0)\theta (\lambda _B\lambda (x_0))\theta (\lambda (x_{\mathrm{\Delta }t})\lambda _B)}{h_𝒜(x_0)}}.`$
Note the similarity with Eq. (1). Strictly speaking $`\theta (\lambda _B\lambda (x_0))`$ is redundant in Eq. (II.3) as $`h_𝒜(x_0)=0`$ if $`\theta (\lambda _B\lambda (x_0))1`$. The last expression in Eq. (II.3) is most suitable for a numerical approach with $`\mathrm{\Delta }t`$ as the time step in a molecular dynamics simulation. Evaluation of Eq. (II.3) requires counting all phase space points which at $`t=0`$ are just about to cross interface $`\lambda _B`$ in one time step and will enter region $`A`$ before $`B`$ when integrating backward in time starting from $`x_0`$. Unfortunately, Eq. (II.3) is not very efficient from a computational point of view because only a very small fraction of phase points close to interface $`\lambda _B`$ actually belong to $`𝒜`$, leading to poor statistics. We can enhance the statistical accuracy by relating the flux through $`\lambda _B`$ to the flux through an interface closer to $`A`$. We therefore introduce a set of $`n`$ non-intersecting interfaces $`\lambda _1,\lambda _2,\lambda _3,\mathrm{}\lambda _n`$, each interface $`\lambda _i`$ closer to $`A`$ than the next interface $`\lambda _{i+1}`$ (see Fig. 3). We define the corresponding phase space regions $`\mathrm{\Omega }_{\lambda _i}\{x|\lambda (x)>\lambda _i\}`$. In this way $`\mathrm{\Omega }_{\lambda _B}`$ is equivalent to our stable state $`B`$, while $`\mathrm{\Omega }_{\lambda _A}`$ is the phase space outside stable state $`A`$. By introducing the following definition

$`\mathrm{\Phi }_{A,\lambda _i}(x_0)`$
$`\underset{\mathrm{\Delta }t0}{lim}{\displaystyle \frac{1}{\mathrm{\Delta }t}}\overline{h}_{A,\mathrm{\Omega }_{\lambda _i}}^b(x_0)\theta (\lambda _i\lambda (x_0))\theta (\lambda (x_{\mathrm{\Delta }t})\lambda _i)),`$ (13)
Eq. (II.3) reduces to

$$k_{AB}=\mathrm{\Phi }_{A,\lambda _B}/h_𝒜.$$
(14)
where $`\mathrm{\Phi }_{A,\lambda _i}`$ denotes the effective positive flux through interface $`\lambda _i`$. The rate constant is thus equal to the effective positive flux through interface $`\lambda _B`$ with the condition the trajectories came directly from $`A`$. Note again that $`\mathrm{\Phi }_{A,\lambda _A}/h_𝒜`$ is equal to the TST rate expression in Eq. (1) in case $`\lambda ^{}=\lambda _A=\lambda _B`$. The effective flux $`\mathrm{\Phi }_{A,\lambda _i}`$ can now be related to the effective flux $`\mathrm{\Phi }_{A,\lambda _{i1}}`$ through an interface $`\lambda _{i1}`$ closer to $`A`$ by (see Appendix A)

$$\mathrm{\Phi }_{A,\lambda _i}(x_0)=\overline{h}_{\mathrm{\Omega }_{\lambda _i},A}^f(x_0)_{\mathrm{\Phi }_{A,\lambda _{i1}}}\times \mathrm{\Phi }_{A,\lambda _{i1}}(x_0),$$
(15)
where $`\mathrm{}_{\mathrm{\Phi }_{A,\lambda _{i1}}}`$ denotes the ensemble average over all phase space points $`x_0`$ for which $`\mathrm{\Phi }_{A,\lambda _{i1}}(x_0)0`$. The factor $`\overline{h}_{\mathrm{\Omega }_{\lambda _i},A(x_0)}^f_{\mathrm{\Phi }_{A,\lambda _{i1}}}𝒫(\lambda _i|\lambda _{i1})`$ is the conditional probability that a trajectory, coming from $`A`$, passes $`\lambda _i`$, given the fact that it has passed the interface $`\lambda _{i1}`$ at an earlier time. By recursively substituting Eq. (15) into Eq. (14) the rate constant can be expressed as a product of conditional probabilities:

$`k_{AB}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Phi }_{A,\lambda _1}}{h_𝒜}}{\displaystyle \underset{i=1}{\overset{n1}{}}}\overline{h}_{\mathrm{\Omega }_{\lambda _{i+1}},A}^f_{\mathrm{\Phi }_{A,\lambda _i}}\overline{h}_{B,A}^f_{\mathrm{\Phi }_{A,\lambda _n}}`$
$``$ $`{\displaystyle \frac{\mathrm{\Phi }_{A,\lambda _1}}{h_𝒜}}{\displaystyle \underset{i=1}{\overset{n1}{}}}𝒫(\lambda _{i+1}|\lambda _i)𝒫(\lambda _B|\lambda _n)`$
$`=`$ $`{\displaystyle \frac{\mathrm{\Phi }_{A,\lambda _1}}{h_𝒜}}\overline{h}_{B,A}^f_{\mathrm{\Phi }_{A,\lambda _1}}{\displaystyle \frac{\mathrm{\Phi }_{A,\lambda _1}}{h_𝒜}}𝒫(\lambda _B|\lambda _1).`$
This expression is the central equation for TIS. Instead of just calculating the individual terms in the product of Eq. (II.3) we can equivalently determine a continuous crossing probability function $`𝒫(\lambda |\lambda _1)`$ for $`\lambda `$ between $`\lambda _1`$ and $`\lambda _B`$. This is reminiscent of umbrella sampling where a free energy difference is usually estimated as a function of a continuous parameter $`\lambda `$ FrenkelSmit . When calculating the ensemble average for $`𝒫(\lambda _i|\lambda _{i1})`$ we can also evaluate $`𝒫(\lambda |\lambda _{i1})`$ for interfaces $`\lambda `$ between $`\lambda _{i1}`$ and $`\lambda _i`$ by dividing the phase space into a finer grid of sub interfaces (see fig. (3)). In this way we acquire useful information without significant extra cost, and, in addition, a measure for the convergence of the ensemble averages. The final monotonically decreasing crossing probability function $`𝒫(\lambda |\lambda _1)`$ can be obtained by matching the histograms from the different ensemble simulations. Techniques commonly applied in umbrella sampling such as overlapping windows between two successive ensemble averages and the use of biasing functions can also be employed here.

## III THE TRANSITION INTERFACE SAMPLING ALGORITHM

Inspection of Eq. (II.3) clearly shows that the TIS rate constant calculation is also a two step procedure. The first step, the effective flux $`\mathrm{\Phi }_{A,\lambda _1}/h_𝒜`$ can be computed by simply running a MD simulation starting with a configuration in $`A`$ and counting the number of effective crossings. For an interface $`\lambda _1`$ close enough to stable state $`A`$ one can obtain a statistically accurate value.

The second part of the calculation consists of evaluating the product of the $`𝒫(\lambda _{i+1}|\lambda _i)`$ ensemble averages for the different interfaces $`\lambda _i`$ in Eq. (II.3). Here we need to sample all paths from region $`A`$ to either $`A`$ or $`\mathrm{\Omega }_{\lambda _{i+1}}`$ that exhibit at least one crossing with interface $`\lambda _i`$. The Monte Carlo moves in TIS are very similar to the shooting move used in the TPS algorithm. The main difference is that the backward and forward integration is abandoned as soon as the edge of either $`A`$ or $`\mathrm{\Omega }_{\lambda _{i+1}}`$ is reached. If the new path is accepted there is only one phase point $`x`$ along this path for which $`\mathrm{\Phi }_{A,\lambda _i}(x)0`$, defining phase space point $`x_0`$. The shifting moves that were required in the original TPS implementation to enable proper sampling and improve statistical accuracy are here unnecessary.

To bootstrap the sampling procedure we first generate an initial path that starts in $`A`$, then crosses the interface $`\lambda _i`$ and finally ends in either $`A`$ or $`\mathrm{\Omega }_{\lambda _{i+1}}`$ (see for more details on initial path generation Ref. Dellago02 ). The phase space point $`x_0`$ is then defined as the first crossing point of this path with interface $`\lambda _i`$. Further, let $`\tau =\text{int}(t/\mathrm{\Delta }t)`$ be the discrete time slice index, and $`\tau ^b\text{int}(t_A^b(x_0)/\mathrm{\Delta }t)`$ and $`\tau ^f\text{int}(t_{A\mathrm{\Omega }_{\lambda _{i+1}}}^f(x_0)/\mathrm{\Delta }t)`$ the forward and backward terminal time slice indices, respectively. Including $`x_0`$, the initial path then consists of $`N^{(o)}=\tau ^b+\tau ^f+1`$ time slices. With these definitions in mind is the TIS algorithm as follows:

1. From the current path with length $`N^{(o)}`$ choose a random time slice $`\tau `$, with $`\tau ^b\tau \tau ^f`$.
2. Change all momenta of $`x_{\tau \mathrm{\Delta }t}`$ by adding small random displacements $`\delta p`$ from a Gaussian distribution. Make sure the total momentum is conserved Dellago02
3. In case of a constant energy (NVE) simulation, rescale the new momenta to the old energy value and continue with step 4. In case of constant temperature (NVT) accept the new momenta (else reject the whole TIS move) with a probability FrenkelSmit :

$`\text{min}[1,\mathrm{exp}\left(\beta \left(E(x_{\tau \mathrm{\Delta }t}^{(o)})E(x_{\tau \mathrm{\Delta }t}^{(n)})\right)\right)].`$
Here, $`E(x)`$ is the total energy of the system at phase space point $`x`$.
4. Integrate equations of motion backward in time by reversing the momenta at time slice $`\tau `$, until reaching either $`A`$ or $`\mathrm{\Omega }_{\lambda _{i+1}}`$. Reject in case of $`\mathrm{\Omega }_{\lambda _{i+1}}`$ else continue with the next step.
5. Integrate from time slice $`\tau `$ forward until reaching either $`A`$ or $`\mathrm{\Omega }_{\lambda _{i+1}}`$. Reject if the entire trial path does not cross interface $`\lambda _i`$, else continue with the next step.
6. Accept the trial path with a probability

$`\text{min}[1,{\displaystyle \frac{N^{(o)}}{N^{(n)}}}],`$
where $`N^{(n)}`$ is the length of the new path. If accepted, replace the old path with the new one.
7. Reassign $`x_0`$ to be the first crossing point with $`\lambda _i`$ and sample the value of $`\overline{h}_{\mathrm{\Omega }_{\lambda _{i+1}},A}^f(x_0)`$ to measure $`𝒫(\lambda _{i+1}|\lambda _i)`$.
8. Repeat from step 1.

As usual in Monte Carlo schemes, any rejection along this route implies counting the old path again in the ensemble average. The acceptance probabilities at step 3 and step 6 are required to satisfy the detailed balance condition (see e.g. Ref FrenkelSmit ).

Instead of generating a complete path and then accepting or rejecting accordingly to the probability at step 6, it is more efficient to determine a maximum path length in advance. Before embarking on the time consuming fourth and fifth step, we first take a uniform random number $`\alpha `$ between $`0`$ and $`1`$ and determine the maximum allowed path length by:

$$N_{\text{max}}^{(n)}=\text{int}(N^{(o)}/\alpha ).$$
(17)
In this way we can directly stop the integration and reject the TIS move as soon the path length $`N^{(n)}`$ exceeds the maximum $`N_{\text{max}}^{(n)}`$. In the course of the TIS simulation the path-length fluctuates. This also means that the average path length becomes automatically shorter when changing from ensemble average $`𝒫(\lambda _{i+1}|\lambda _i)`$ to ensemble average $`𝒫(\lambda _i|\lambda _{i1})`$ closer to $`A`$.

The algorithm presented here does not require shifting moves because there is only one unique $`x_0`$ phase point along each pathway. However, one could consider the use of path-reversal moves as they have negligible computational cost and can sometimes facilitate ergodic sampling Dellago02 .

## IV NUMERICAL RESULTS

### IV.1 The model

We tested the TIS algorithm on a simple diatomic bistable molecule immersed in a fluid of purely repulsive particles. Such a system has previously been used in illustrating TPS rate constant calculations TPS99  and is therefore a good starting point for a comparison between the two methods. The system consists of $`N`$ particles in 2 dimensions with interactions given by a pairwise Lennard-Jones (LJ) potential truncated and shifted at the minimum, often referred to as the Weeks-Chandler-Andersen (WCA) potential WCA

$$V_{WCA}(r)=\{\begin{array}{cc}4ϵ[(r/\sigma )^{12}(r/\sigma )^6]+ϵ\hfill & \text{if }rr_0\hfill \\ 0\hfill & \text{if }r>r_0,\hfill \end{array}$$
(18)
where $`r`$ is the interatomic distance, and $`r_02^{1/6}\sigma `$. Throughout this section reduced units are used so that $`ϵ`$ and $`\sigma `$, respectively the LJ energy and length parameters, as well as the mass of the particles are equal to unity. The LJ unit of time $`(m\sigma ^2/ϵ)^{1/2}`$ is therefore also unity. In addition, two of the $`N`$ particles are interacting through a double well potential

$$V_{dw}(r)=h\left[1\frac{(rr_0w)^2}{w^2}\right]^2.$$
(19)
This function has two minima separated by a barrier of height $`h`$ corresponding to the two stable states of the molecule: a compact state for $`r=r_0`$ and extended state for $`r=r_0+2w`$. For a high enough barrier, transitions between the states become rare and the rate constant is well defined. Hence, this system provides a useful test case for the TPS and TIS algorithms.

The system is simulated at a constant energy $`E`$ in a simulation square box with periodic boundary conditions. The total linear momentum is conserved and is set zero for all trajectories. The equations of motion are integrated using the velocity Verlet algorithm with a time step $`\mathrm{\Delta }t=0.002`$. As in Ref. TPS99  we focus here on the computation of the rate constant for the isomerization reaction of the dimer from the compact state to the extended state. In the following section we describe general simulation details. In Section IV.3 we discuss the results for a system with a high enough barrier to avoid recrossings. Subsequently, we reproduce the simulations from Ref. TPS99  in Section IV.4. These results do show recrossings, and we discuss the consequences for TPS and TIS.

### IV.2 Methodology

The TPS rate constant calculation evaluates the two factors in Eq. (6) separately as explained in Sec. II.2. The first term in Eq. (6) is the ratio between the plateau value of the reactive flux correlation function $`\dot{h}_B(T)_{A,H_B(T)}`$ and the correction $`h_B(t^{})_{A,H_B(T)}`$. The second term $`C(t^{})`$ requires an umbrella sampling simulation in the form of a series of window calculations. An order parameter is chosen to define the characteristic functions of the stable states and to partition phase space in windows for the umbrella sampling. Besides shooting and shifting Monte Carlo moves to generate new paths in the transition path sampling we also employ a *diffusion* move that shifts the path by one time slice in arbitrary direction. This move is computationally very cheap but increases the statistics of the correlation functions. In all our simulations we therefore set the percentages for shooting, shifting and diffusion to 5%,10% and 85%, respectively. The parameters involved are always gaged such that the acceptance ratio is around 40% for shooting and shifting moves, ensuring an optimum efficiency of the sampling TPS99 .

The TIS method involves a direct determination of the flux and the calculation of the crossing probability functions $`𝒫(\lambda _i|\lambda _{i1})`$ between a series of successive interfaces as given by Eq. (II.3). The flux term in Eq. (II.3) is computed by means of a straightforward MD simulation starting in state $`A`$ and counting the number of effective positive crossings through interface $`\lambda _1`$, i.e. when the trajectory is directly coming from $`A`$. The second term in Eq. (II.3) is computed using the TIS algorithm of Sec. III. The basic requirement is a definition of a set of interfaces partitioning the phase space. Between these interfaces we defined a finer grid of sub-interfaces to construct the crossing probability function $`𝒫(\lambda |\lambda _1)`$. As in the TPS calculation we adjusted the momentum displacement for the shooting move to give an acceptance of about 40%.

Many parameters are involved in the two methods and to compare the relative efficiency we measured the CPU-time required for an arbitrary fixed error of $`2.5\%`$ for each step in both the TPS and TIS calculations under the same computational conditions (1Ghz AMD Athlon). In both methods the final rate constant consists of a product of factors which have to be calculated independently. For each factor we performed $`M`$ simulation blocks of $`N`$ Monte Carlo cycles and adjusted $`N`$ such that after $`M`$ block averages the relative standard deviation of each term in Eq. (4) and (11) was 2.5%. The total CPU-time is given by summing the individual 2.5% error CPU-times for each factor. The final error in the rate constants is obtained by the standard propagation rules using all simulation results (i.e. not only the ones for the 2.5% error CPU time calculation).

### IV.3 System with High Energy Barrier

This system had a total number of particles $`N=25`$, and a total energy $`E=25`$. The square simulation box was adjusted to give a number density of 0.7. The barrier height was $`h=15`$ and the width-parameter $`w=0.5`$, so that the minima of $`V_{dw}(r)`$ were located at $`r1.12`$ and $`r2.12`$ while the top of the barrier was at $`r1.62`$. In the TPS rate calculation we defined stable states $`A`$ and $`B`$ as $`r<r_A=1.5`$ and $`r>r_B=1.74`$, respectively. We computed the correlation function $`h_B(t)_{A,H_B(T)}`$ using TPS with a fixed path length $`T=2.0`$. The correlation function is shown in Fig. 4 together with its time derivative, the reactive flux. The latter function clearly displays a plateau. Next, we chose four different $`t^{}=0.1,0.3,1.0,2.0`$ and performed umbrella sampling simulations using 8 windows to calculate $`C(t^{})`$. In each window we measured the probability to find the path’s end point $`r(t^{})`$ at a certain value of $`r`$. These probability histograms were rematched and normalized. The final probability functions are shown in Fig. 5. Integration of the area under the histogram belonging to region B leads to $`C(t^{})`$ and finally to the rate constant TPS99 . In Table 1 we give the values of the different contributions to the rate constant given by Eq. 6, together with the rate constant. We report the average relative computation time needed to reach the $`2.5\%`$ error (see Sec. IV.2) in Table 2.

For the TIS calculations we use the same order parameter $`r`$ and the same definition for region $`B`$, i.e. interface $`\lambda _B`$ is set at $`r=1.74`$. Stable state $`A`$ was defined by setting $`\lambda _A=\lambda _1`$ at $`r=1.24`$. This interface is closer to the basin of attraction than the TPS stable state definition but yields a higher flux term $`\mathrm{\Phi }_{A,\lambda _1}/h_𝒜`$ and gives better statistics. Note that the different definition of stable state $`A`$ does not change the final rate constant, as the overall state $`𝒜`$ does not sensitively depends on this definition. The flux term is calculated by straightforward NVE MD. As $`\lambda _A`$ is equal to $`\lambda _1`$ every positive crossing of this interface is counted in the flux because all trajectories must by default come directly from A. The conditional crossing probabilities $`𝒫(\lambda _{i+1}|\lambda _i)`$ in Eq. (II.3) are calculated for $`n=5`$ interfaces between the stable states (see fig. 6). Between these interfaces we impose a finer grid to obtain the entire crossing probability function. The results for each stage and the final rate constant are shown in Table 1. The rate constants of both methods agree within the statistical accuracy, showing that the TIS method is correct. In Table 2 we give the relative computation time to reach the $`2.5\%`$ error for each term.

In comparing both methods we have to recall that the efficiency of TPS depends strongly on the choice of $`t^{}`$. On the one hand the umbrella calculation of $`C(t^{})`$ is much faster for low values of $`t^{}`$. But on the other hand the error in the correction term $`h_B(t^{})_{A,H_B(T)}`$ increases for lower $`t^{}`$. As a result there is an optimum $`t^{}`$ for the error/CPU-time ratio, in this case approximately at $`t^{}=0.3`$. Even for this optimized situation the TIS calculation is about two times faster. One could object that the correlation function in Fig. 4 has reached a plateau for $`t=1.5`$ already, reducing the TPS computation time by a factor $`3/4`$. But the choice for a path length $`T=1.5`$ can not be taken without a-priori knowledge. The first term in Eq. (6) implicitly depends on the path length $`T`$. Changing $`T`$ would alter the ensemble and might result in a different shape of the flux correlation function. We did not check this in detail, but we believe that $`T`$ cannot be chosen much smaller without introducing systematic errors. Furthermore, we emphasize here that we put much more effort in optimizing the TPS algorithm by tuning $`t^{}`$, the windows, the ratio between shooting, shifting and diffusion moves than we did for TIS.

Figure 7 shows the histograms of path lengths for each TIS ensemble calculation and shows why TIS is faster. Sampling paths of fixed length with TPS results in spending unnecessary computation time inside the initial and final stable regions A and B. In the TIS algorithm instead every path is adapted to its minimum length. Bringing the interface in closer to $`A`$ reduces these transition times. TIS optimizes itself during the simulation.

### IV.4 System with Low Energy Barrier

In order to compare with previous results, we adopted the parameters from Ref. TPS99 . The total number of particles was $`N=9`$, the total energy was $`E=9`$ and the square simulation box was adjusted for a number density of 0.6. The barrier height is $`h=6`$ and the width-parameter is $`w=0.25`$. Minima are at $`r1.12`$ and $`r1.62`$, while the top of the barrier is at $`r1.37`$. This barrier is much lower than in the previous section resulting in more frequent transitions. An approximate rate constant could even be achieved by straightforward MD simulations.

For the TPS calculations we defined the stable states $`A`$ and $`B`$ by $`r<r_A=1.30`$ and $`r>r_B=1.45`$, respectively TPS99 . Using standard TPS simulation we computed the correlation function $`h_B(t)_{A,H_B(T)}`$ with a total path length $`T=2`$ (shown in Fig. 8). Next, we measured the probability histograms to find the paths end point at a certain order parameter value $`r`$ for four different times $`t^{}=0.1,0.4,0.8,2.0`$, using five windows TPS99  (see Fig. 9). As described in the previous section, matching the probability histograms and subsequent integration leads to $`C(t^{})`$. The resulting final rate constants, shown in Table 3, are comparable with the results of Ref. TPS99 , but more accurate. We will discuss these values after giving the results of TIS.

Figure 10 shows that fast recrossings can occur for a low barrier, implying that $`r`$ alone is not sufficient as an order parameter to define the stable states in the simulations. Apparently, this does not effect the TPS results much, but it is very important for TIS because of the assumption that stable region B is really stable and recrossings do not take place. To ensure the stability of the TIS stable states we chose a new order parameter that not only depends on the inter-atomic distance $`r`$ in the dimer but also on a kinetic term, given by $`\dot{r}`$. The stable states can then be defined by

$`E_d(r,\dot{r}){\displaystyle \frac{\dot{r}}{4}}+V_{dw}(r)`$
$`xA\text{ if }r<1.37\text{ and }E_d(r,\dot{r})1.5`$
$`xB\text{ if }r>1.37\text{ and }E_d(r,\dot{r})1.5,`$ (20)
where $`E_d`$ is the sum of the kinetic and potential energy of the dimer that has a reduced mass of $`1/2`$. In the $`\{r,\dot{r}\}`$ plane these stable states form a D-shape and an inverse D-shape regions for $`A`$ and $`B`$ respectively (see Fig. 11). Crossing the interface $`\lambda _A`$ or $`\lambda _B`$ implies that the vibrational energy is decreased below the threshold, $`E_d=1.5`$. This threshold is made low enough to make fast recrossings to the other state unlikely. However, if we would have chosen it too low the paths would have become very long. We evaluated the crossing probability function in Eq. (II.3) for $`n=3`$ interfaces. The entire crossing probability function was obtained by partitioning the phase space in sub-interfaces of the form $`r=\lambda `$ and $`E_d(r,\dot{r})=\lambda `$ as shown in Fig. 11. Note that in TIS multidimensional or multiple order parameters can be used in one simulation without a problem. This is more difficult in TPS, where a proper mapping of the complete phase space is required. Figure 12 shows the final rematched crossing probability. The monotonically decreasing function tends to reach a plateau on approaching the last interface. The last two values are not exactly equal but differ by 0.03%, indicating that a small fraction of the paths crossing the one but last sub-interface still succeed to return to $`A`$ without crossing $`\lambda _B`$. This difference is comparable with the chance of a new independent transition (given by the rate constant). Note that without the kinetic energy definition for the stable states Eq. (20), the final crossing probability and thus the rate constant would have been overestimated by a factor $`5/4`$.

For the effective flux $`\mathrm{\Phi }_{A,\lambda _1}/h_𝒜`$ calculation we performed MD simulations as described in Sec IV.2. In contrast to the high barrier case, $`\lambda _1`$ is not equal $`\lambda _A`$, and not all positive crossings with $`\lambda _1`$ are effective crossings. We counted only the first crossing when the system left region $`A`$ and waited until the system fell back to region $`A`$ before counting a new crossing. As the MD trajectory sometimes displayed a spontaneous transition to region $`B`$, we stopped the simulation and started again by replacing the system in a randomized configuration of $`A`$. Table 3 shows the final values and the corresponding errors of these calculations. The relative computation time for each term is detailed in table 4.

If we compare the final results of table 3 we see that the efficiency of TIS is more than nine times better than the TPS efficiency for $`t^{}=2`$, and more than two times better than TPS value for $`t^{}=0.8`$. But the TPS $`t^{}=0.1`$ and $`t^{}=0.4`$ efficiencies are about 20 % better than TIS. When we compare the rate constants, however, we notice that the TPS results for different $`t^{}`$ do not agree. Among the TPS rate constants only the $`t^{}=2`$ case is consistent with the TIS result. We believe that the $`t^{}=0.1`$ and $`t^{}=0.4`$ results suffer from systematic errors. For instance, for the shorter paths the TPS simulations might not be completely ergodic. Another explanation might be that a path length of $`T=2`$ is too short to allow convergence of the reactive flux. In the TIS calculation the average path length in the three interface simulations, from the closest to $`B`$ to the closest to $`A`$, is, respectively, 7.4, 4.3, and 0.63; much longer than the TPS path length (see Fig.13). It is therefore surprising that the TPS approach with the simple stable state definition and very short paths still gives approximately the right rate constant. And indeed, when we computed the TPS correlation function with the TIS state definitions Eq. (20), we found that the path length had to be at least $`T=20`$ to see a plateau. We think that TPS works even with the simple state definitions and the short paths because both positive flux and negative flux terms contribute to Eq. (6). The TPS algorithm collects many paths of which some are not real transitions, but fast recrossings. The cancellation of positive and negative terms of these fast recrossing paths ensure the (almost) correct final outcome. In TIS each path must be true transition event and contributes as a positive term in rate equation (II.3), enhancing the convergence. This explains that the CPU time for the TIS calculation despite the much longer paths is still comparable with TPS one for low $`t^{}`$. We note that the path ensemble using the more strict stable state definition is of course more useful in the analysis of the reaction mechanism.

For a more accurate comparison of the computation time we must keep the systematic errors lower than the statistical errors. In other words, we have to make sure that the results are converged. To test the convergence of the flux correlation function in TPS we can derive the following equality from Eq. (6) :

$$\frac{h_B(t^{})_{A,H_B(T)}}{h_B(t^{\prime \prime })_{A,H_B(T)}}=\frac{C(t^{})}{C(t^{\prime \prime })}.$$
(21)
This equation is valid for any $`t^{},t^{\prime \prime }<T`$ if $`T`$ is large enough. We found that the equality does not hold for the system with the low barrier, indicating that $`T`$ is too low in the TPS calculation. Further examination of the flux correlation function $`h_B(t^{})_{A,H_B(T)}`$ reveals that the apparent plateau has in fact a small positive slope. Calculations for higher values of $`T`$ suggest that one has to increase the path length at least to $`T=8`$ to convergence to a plateau. With this in mind we think that the TIS computation is about a factor five more efficient than the TPS algorithm for the model system with the low barrier.

## V CONCLUSION

We developed a novel method, named transition interface sampling, for the calculation of rate constants based on transition path sampling concepts. Just as the original transition path sampling, the new method enables the calculation of rate constants of transitions between stable states separated by high free energy barriers without prior knowledge of the reaction coordinate. The new algorithm is different in spirit from the rate constant calculation that was introduced in RefsTPS98a ; TPS99 . In TPS the time correlation function $`h_A(x_0)h_B(x_t)/h_A`$ is determined for a single time using an umbrella sampling scheme followed by a calculation of the reactive flux prefactor in a separate path sampling simulation. The path length used in this simulation has to be long enough for the plateau to be reached. The TIS method advocated here calculates the flux correlation directly by measuring the fluxes through a number of different interfaces and relating the flux through one interface to the next one. The big advantage of a flux instead of a correlation function is that trajectories going through the interfaces all contribute to the rate whereas in TPS there are recrossings to be counted. In addition, the new method improves the original TPS method on a several other points. Once the interface is reached, the integration of motion can be stopped instead of going all the way to region $`B`$. In this way the TIS algorithm adapts itself to the optimal path length. One does not have to optimize the new method as much as TPS, where one has to find the optimal $`t^{}`$ value and a proper balance between shooting and shifting. Besides being faster, the concept of calculating a flux comes natural with the rate constant definition, and implementation of the algorithm is hence simpler. Also, multidimensional or even discrete order parameters can easily be implemented in TIS. In the illustrative example we showed that we can obtain an increase in efficiency of at least a factor of two to five with respect to the TPS method used in Ref. TPS99 .

However, one has to be more careful in the definitions of the stable states, meaning that stable states have to be really stable. In TPS the choice of stable states is a bit more flexible as the final rate constant consists of cancellation of positive and negative terms. In sec IV.4 we showed how this problem for TIS can be solved by defining stable regions that explicitly depend on kinetic energy terms.

In summary, we believe that the TIS algorithm can make the rate constant calculation of many processes feasible that were hitherto difficult to obtain. For instance, chemical reactions in solution, isomerization of clusters and conformational transitions in biomolecules. In a future publication we will report on these, more complex, applications.

###### Acknowledgements.

T.S.v.E acknowledges NWO-CW (Nederlandse Organisatie voor Wetenschappelijk Onderzoek, Chemische Wetenschappen). P.G.B acknowledges support from the FOM (Stichting Fundamenteel Onderzoek der Materie). We thank Evert Jan Meijer for carefully reading the manuscript.

## Appendix A FLUX RELATION

In this appendix we show how the effective flux $`\mathrm{\Phi }_{A,\lambda _i}`$ can be related to the effective flux $`\mathrm{\Phi }_{A,\lambda _{i1}}`$ through an interface $`\lambda _{i1}`$ closer to $`A`$. If at time $`t=0`$ a trajectory passes interface $`\lambda _i`$ while having started in $`A`$ some time earlier, there must always be an unique time when it passed interface $`\lambda _{i1}`$ for the first time. Therefore we can write:

$`\mathrm{\Phi }_{A,\lambda _i}(x_0)=\mathrm{\Phi }_{A,\lambda _i}(x_0){\displaystyle _0^{t_{A\mathrm{\Omega }_{\lambda _i}}^b(x_0)}}dt\mathrm{\Phi }_{A,\lambda _{i1}}(x_t)`$ (22)
and hence,

$`\mathrm{\Phi }_{A,\lambda _i}(x_0)=`$
$`{\displaystyle _0^{\mathrm{}}}dt\mathrm{\Phi }_{A,\lambda _{i1}}(x_t)\mathrm{\Phi }_{A,\lambda _i}(x_0)\theta (t_{A\mathrm{\Omega }_{\lambda _i}}^b(x_0)t)`$
$`={\displaystyle _0^{\mathrm{}}}dt\mathrm{\Phi }_{A,\lambda _{i1}}(x_0)\mathrm{\Phi }_{A,\lambda _i}(x_t)\theta (t_{A\mathrm{\Omega }_{\lambda _i}}^b(x_t)t)`$
$`=\mathrm{\Phi }_{A,\lambda _{i1}}(x_0){\displaystyle _0^{\mathrm{}}}dt\mathrm{\Phi }_{A,\lambda _i}(x_t)\theta (t_{A\mathrm{\Omega }_{\lambda _i}}^b(x_t)t)`$
$`=\mathrm{\Phi }_{A,\lambda _{i1}}(x_0){\displaystyle _0^{t_{A\mathrm{\Omega }_{\lambda _i}}^f(x_0)}}dt\mathrm{\Phi }_{A,\lambda _i}(x_t)`$
$`=\mathrm{\Phi }_{A,\lambda _{i1}}(x_0)\overline{h}_{\mathrm{\Omega }_{\lambda _i},A}^f(x_0).`$ (23)
The one but last equation follows because for each phase point $`x`$ and phase space region $`\mathrm{\Omega }`$ it can be shown that $`t>t_\mathrm{\Omega }^f(x)t_\mathrm{\Omega }^b(f(x,t))t\theta \left(t_\mathrm{\Omega }^b(f(x,t))t\right)=0`$. We rewrite the last expression of Eq. (23) as a different ensemble average:

$`\mathrm{\Phi }_{A,\lambda _{i1}}(x_0)\overline{h}_{\mathrm{\Omega }_{\lambda _i},A}^f(x_0)=`$
$`={\displaystyle \frac{\mathrm{\Phi }_{A,\lambda _{i1}}(x_0)\overline{h}_{\mathrm{\Omega }_{\lambda _i},A}^f(x_0)}{\mathrm{\Phi }_{A,\lambda _{i1}}(x_0)}}\times \mathrm{\Phi }_{A,\lambda _{i1}}(x_0)`$
$`\overline{h}_{\mathrm{\Omega }_{\lambda _i},A}^f(x_0)_{\mathrm{\Phi }_{A,\lambda _{i1}}}\times \mathrm{\Phi }_{A,\lambda _{i1}}(x_0),`$ (24)
where $`\mathrm{}_{\mathrm{\Phi }_{A,\lambda _{i1}}}`$ denotes the ensemble average over all phase points $`x_0`$ for which $`\mathrm{\Phi }_{A,\lambda _{i1}}(x_0)0`$. The last equality gives rise to Eq. (15).


text:
# A new infrared band in Interstellar and Circumstellar Clouds: C4 or C4H?

## 1 Introduction

Long carbon chain radicals have been detected in interstellar and circumstellar clouds through their pure rotational spectrum at radio wavelengths (C5H, C6H, C7H, C8H, H2C3, H2C4: Cernicharo et al. 1986, 1987, Guélin et al. 1987, 1997, Cernicharo and Guélin 1996). Polyynes such as C2H2, C4H2, C6H2 and single aromatic species such as benzene, have been detected in the protoplanetary nebula CRL618 (Cernicharo et al. 2001a&b) supporting the idea that much complex carbon-rich molecules could be formed in space. However, the mechanisms allowing the growth of carbon-rich molecules are still poorly known and the full set of molecules that could participate in the chemical reactions leading to the formation of large complex carbon-rich species has yet to be identified. Among these ”building blocks” the Cn chains have attracted the interest of laboratory spectroscopists and astronomers. C3 has been observed at optical wavelengths in the atmospheres of cool stars (see, e.g., Zuckerman et al. 1976), and it has been identified in the envelope of IRC+10216 through its $`\nu _3`$ antisymmetric stretching mode in the mid-infrared by Hinkle et al. (1988) who derived $`\chi `$(C3)$``$10-6. This key species has recently been found in the ISM (Cernicharo et al. 1996, 1997a, 2000, Giesen et al. 2001) through its $`\nu _2`$ low frequency bending mode (Schmuttenamer et al. 1990) and in the diffuse interstellar clouds through its electronic transition near 4052 Ȧ (Maier et al. 2001; Roueff et al. 2002). It has been suggested that C3 could be involved in the formation of the diffuse interstellar bands (Douglas, 1977; Clegg and Lambert 1982). ISM observations in Cernicharo et al. (2000) put a lower limit of 10-8 to the C3 abundance. In addition, C5 has been detected in the C-rich star IRC+10216 with an abundance ratio C3/C$`{}_{5}{}^{}`$10 (Bernath et al. 1989). Less is known about numbered carbon chains, but C4 could well be a very abundant molecule both in C-rich evolved stars and in the ISM.

In this Letter we report the detection of a blend of lines at 57.5 $`\mu `$m that we tentatively assign to the $`\nu _5`$ bending mode of C4. The feature has been observed in Sgr B2, IRC+10216, CRL618, CRL2688, and NGC7027. We discuss the possibility that the carrier could be another different carbon chain : C4H, although photodetachement experiments suggest that its lowest energy bending mode could be at higher frequencies (Taylor et al. 1998).

## 2 Observations and Results

We have searched for the emission/absorption of the low energy bending modes of polyatomic molecules toward interstellar and circumstellar clouds using the Long Wavelength spectrometer (LWS; Clegg et al. 1996, Swinyard et al. 1996) on board ISO1 (Kessler et al. 1996). Our main target for this search has been the C-rich evolved star IRC+10216.

The LWS grating data of IRC+10216 were analyzed and modeled by Cernicharo et al. (1996). The data have been reprocessed following the last pipeline product (TDT 19800158). No significant differences have been found with the already published analysis. After identifying all the HCN (including several vibrational states) and CO lines, a conspicuous feature composed by a blend of several lines (see Figure 1, top panel) was identified between 56 and 58 $`\mu `$m. These lines could not be assigned to any known species and no other similar patterns were found in the ISO far-infrared spectrum of IRC+10216. We have checked for all pure rotational lines of light species that could be potentially abundant without success. Hence, the only possible explanation is that the feature corresponds to a bending mode of an abundant heavy molecule, and most likely to a carbon chain. We searched for the same feature toward other objects selected from the same criteria of our previous searches of C3, C4H2, C6H2, and benzene (Cernicharo et al. 1996, 1997a, 2000, 2001a,b), i.e., IRC+10216, CRL618, CRL2688, NGC7027, and Sgr B2 (among other sources).

For Sgr B2 M we have used the LWS grating spectra obtained during revolutions 287 and 494 (TDTs 28701401, 287012130, 287012131, and 49400302). Some of these data have been already presented (Cernicharo et al. 1997b, Goicoechea & Cernicharo 2001). All together represent a total of 32 individual scans. The feature at 57.5 $`\mu `$m is clearly present. There are many other LWS grating observations of Sgr B2 with a reduced number of scans. We have inspected all of them and, again, the feature appears in each individual scan. Sgr B2 N was observed in the course of a fast mapping of the Sgr B2 star forming region with the LWS spectrometer and only four scans are available (see Cernicharo et al. 1997b).

For the Planetary Nebula NGC7027 we have used all the LWS grating data taken by ISO (see Herpin et al. 2002 for details). The data for the Proto-Planetary Nebula (PPN) CRL618 (TDTs 68800302 and 68800450) are those published by Herpin and Cernicharo (2000). For CRL2688, a young PPN, the TDT 02101504 data have been used. A previous analysis of the ISO/LWS spectrum of this source has been presented by Cox et al. (1996).

All the spectra have been analyzed using the ISO spectrometers data reduction package ISAP2. The data products correspond to version 10 of the pipeline. A polynomial baseline had to be removed from the final spectrum in each source. For some objects this baseline removal is critical to correctly extract the shape and intensity of the absorption features. Due to the large number of lines in Sgr B2 M and N spectra, the determination of the continuum level is not obvious. For sources with emission lines, although baseline removal is less critical, the global behavior of detector SW2 does not follow the expected shape for continuum emission and a polynomial of fourth degree has been removed after blanking the wavelength intervals where features above 5 sigma were present.

11footnotetext:  Based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with participation of ISAS and NASA.22footnotetext: The ISO Spectral Analysis Package (ISAP) is a joint development by the LWS and SWS Instruments Teams and Data Centers. Contributing institutes are CESR, IAS, IPAC,MPE, RAL, and SRON.

Figure 1 shows the 53-62 $`\mu `$m spectra of the selected sources. We have also analyzed the available observations for several other sources : O-rich stars, Orion, Sgr A, bright molecular sources, and even Mars to check for possible instrumental effects. The new infrared feature was not detected in any of these sources. We have also examined the spectral response of the LWS-SW2 filter in order to check for possible instrumental contributions to the blend of lines found around 57 $`\mu `$m but there is nothing at this wavelength (Garcia–Lario, private communication). Hence, we conclude that the blend of lines is real. The pure rotational lines of CO and HCN are indicated in the top panel of Figure 1. The model of Cernicharo et al. (1996) predicted a low flux for these high-J lines. Only in IRC+10216 some of the CO lines (J=43-42, 44-43 and 45-44) are marginally detected at a 3$`\sigma `$ level.

The observed band consists of a blend of at least 4 spectral lines or sub-bands (at 57.0, 57.6, 57.9 and 58.2 $`\mu `$m). There is an additional feature separated from the band at 56.15 $`\mu `$m that is detected in C-rich evolved stars but is less evident in Sgr B2 N and M. The J=47-46 line of CO could only make a modest contribution to this feature in C-rich evolved stars since other high-J CO lines are weak – even in IRC+10216. In Sgr B2 M this feature is blended with the absorption produced by NH3 (J=9-8). The 56.15 $`\mu `$m feature could be one of the fine structure lines of \[SI\] at 56.31 $`\mu `$m (3P03P1) but its companion \[SI\] (3P13P2) line at 25.249 $`\mu `$m is missing in all objects. The feature could also be the R(1) line of HD. However, it is unlikely to find such a high intensity for this line in evolved stars. We think, that if instrumental, the feature could arise from Uranus, the planet used as calibrator for the LWS/grating spectrometer. However, its intensity changes by more than a factor 10 between NGC7027 and IRC+10216 or CRL618 (the continuum flux at 58 $`\mu `$m is nearly the same for NGC7027 and CRL618, 1.15 and 0.7 10-16 W cm$`{}_{}{}^{2}\mu `$m-1 respectively, and 4 and 3 10-16 W cm$`{}_{}{}^{2}\mu `$m-1 for IRC+10216 and CRL2688 respectively). Therefore, the 56.2 $`\mu `$m feature in IRC+10216, CRL618, NGC7027 and CRL2688 cannot be explained alone by a possible instrumental contribution from the calibrator. Figure 1 shows that the 56.2 $`\mu `$m line follows the intensity of the 57-58 $`\mu `$m features. In the following we will consider that all them belongs to the same molecular band.

Two additional features at 53.7 and 54.4 $`\mu `$m in the spectrum of IRC+10216 (labelled U1 and U2 in the top panel of Figure 1), and another one at 52.89 $`\mu `$m (not shown), are also detected in NGC7027 but not in the other sources. These three lines remain, so far, unidentified. The contribution of OH (Goicoechea & Cernicharo 2002) and the maximum possible contribution from H2O (Cernicharo et al. 2002, in preparation) to the spectra of Sgr B2 N and M are indicated in the panel corresponding to Sgr B2 M. The contribution of H2O to the spectrum of CRL618 at 58$`\mu `$m is negligible (Herpin and Cernicharo 2000).

## 3 Discussion

The fact that the 58 $`\mu `$m feature is detected in all the selected C-rich post-AGB objects clearly points toward a carbon-rich molecule as the carrier of the band. Several molecules with four atoms have low frequency bending modes around 160-220 cm-1. HC3N has a low energy bending mode, $`\nu _7`$, at 222.413 cm-1 ($``$45 $`\mu `$m). Dyacetylene, C4H2, recently detected in CRL618 with a very large abundance (see Cernicharo et al. 2001a&b), has a low energy bending mode, $`\nu _9`$, also at $``$45 $`\mu `$m. Another example of molecule with bending modes in this wavelength range is C4H. Pure rotational lines from the $`\nu _7=1,2`$ bending states of C4H have been detected in IRC+10216 by Guélin et al. (1987) and Yamamoto et al. (1987). The frequency of this mode is expected from ab initio calculations (Graf et al. 2001) to be around 178 cm-1. Yamamoto et al. (1987) have estimated a frequency of 131 cm-1 from the l-doubling constant, q7, of the $`\nu _7`$=1,2 vibrational levels. A 20% of error could be expected on this value (Mikami et al. 1989). Photoelectron spectra of C4H indicated, however, a larger frequency for this mode (226 cm-1; Taylor et al. 1998). The low lying $`{}_{}{}^{2}\mathrm{\Pi }`$ electronic state of C4H is at $``$468 cm-1 above the ground (Taylor et al. 1998). Both transitions, the $`\nu _7=10`$ and the electronic $`A^2\mathrm{\Pi }X^2\mathrm{\Sigma }`$ are in the wavelength coverage of the Short Wavelength Spectrometer (SWS) on-board ISO. However, IRC+10216 is so bright in the mid-infrared that the expected range for these transitions is covered by strong ripples and not obvious features have been found so far (Cernicharo et al. 1999). The shape of the band in Figure 1 allows to rule out a $`{}_{}{}^{1}\mathrm{\Pi }^1\mathrm{\Sigma }`$ transition, i.e., the bending modes of HC3N, HC5N, CnH2 with n even, Cn with n odd, etc. These molecules have a $`{}_{}{}^{1}\mathrm{\Sigma }`$ ground electronic state and, at the limited resolution of the LWS grating spectrometer, a strong Q-branch with weak P and R-branches are expected for their bending modes. The HC3N $`\nu _7`$ and C4H2 $`\nu _9`$ modes are present in the LWS spectrum of CRL618 (to be published elsewhere), and correspond to these expectations.

Taking into account the large abundance of C3 found in Sgr B2 and IRC+10216, C4 could be a good candidate to explain the IR feature under investigation. Linear C4 was first observed in the laboratory using electron spin resonance spectroscopy of carbon clusters in matrices by Graham et al. (1976). They were able to derive the spin-coupling and spin-doubling constants, $`\lambda `$ and $`\gamma `$, for the ground vibrational state. C4 was first detected and characterized in the gas phase by Heath and Saykally (1991) by observing several ro-vibrational lines of the $`\nu _3`$ antisymmetric stretching mode. Moazzen-Ahmadi et al. (1994) observed the same band and the $`\nu _3+\nu _5\nu _5`$ hot band and determined the rotational constants, B and D, for both the ground and the $`\nu _5`$ states. The l-doubling constant, q5, for the latter vibrational state was also derived leading these authors to estimate a frequency of 160$`\pm `$4 cm-1 for the $`\nu _5`$ mode. The $`\nu _5`$ bending mode of C4 has been measured in matrix experiments at a wavelength of 58 $`\mu `$m (Withey et al. 1991), very close to that of the newly found IR feature. They observed two different lines separated by 2 cm-1 (170 and 172.4 cm-1) that they interpreted as two different trapping sites for C4 in the Ar matrix. Both values could correspond to the frequency of the $`\nu _5`$ mode of C4 but the authors adopted the value of 172.4 cm-1. This value agrees well with the series of peaks observed in the photoelectron spectroscopy of Arnold et al. (1991). The difference between the frequency of the new IR band and that derived from the matrix experiments is 2.6 cm-1, which could be consistent with a frequency shift introduced by the Ar matrix of 1.5%. Ab initio calculations predict wavenumbers in the 175-214 cm-1 range (47-57 $`\mu `$m) and an infrared intensity of 46.2 km mol-1 \[A($`\nu _5=10`$) $``$ 0.25 s-1; Martin et al. 1991; Watts et al. 1992; Taylor and Martin, 1996; see also Van Orden and Saykally, 1998\]. Due to the complexity of the C4 electronic energy levels and to the limited accuracy level of ab initio calculations, this intensity value has to be taken into account as a crude estimation.

The fact that the observed band (see Figure 1) contains several spectral features suggests a transition between terms of higher spin multiplicity, in particular toward a vibronic transition $`{}_{}{}^{2}\mathrm{\Pi }^2\mathrm{\Sigma }`$ or $`{}_{}{}^{3}\mathrm{\Pi }^3\mathrm{\Sigma }`$, i.e., just the type of bands we can expect for C4H and C4 which have $`{}_{}{}^{2}\mathrm{\Sigma }`$ and $`{}_{}{}^{3}\mathrm{\Sigma }_{g}^{}`$ ground electronic states respectively. Consequently the lowest bending modes of these species, $`\nu _7`$ for C4H and $`\nu _5`$ for C4, will have $`{}_{}{}^{2}\mathrm{\Pi }`$ and $`{}_{}{}^{3}\mathrm{\Pi }`$ vibronic character respectively. An argument against C4H is the following, while pure rotational lines of its bending mode have been detected towards IRC+10216 (Guélin et al. 1987; Yamamoto et al. 1987) these lines are missing in all other objects reported here (except in CRL618 where several weak lines from these states have been detected; Cernicharo et al. 2002, in preparation). Moreover, the abundance of C4H in Sgr B2 is rather low (Guélin, Mezaoui & Friberg, 1982), while C3 has a large abundance (Cernicharo, Goicoechea and Caux, 2000). Hence, it is tempting to assume that the observed features in post-AGB stars and in the ISM correspond to the $`\nu _5`$ mode of C4, although a definitive assignment will require better spectral resolution data and/or additional laboratory information on this species. The main problem in assigning the band to the $`\nu _5`$ mode of C4 is the large spin-orbit interaction constant, ASO, that should have the $`{}_{}{}^{3}\mathrm{\Pi }`$ vibronic state in order to reproduce the observations. As the ground electronic state of C4 is $`{}_{}{}^{3}\mathrm{\Sigma }`$, a very small value for ASO should be expected. However, if the molecule has a low lying excited electronic state with $`\mathrm{\Lambda }\mathrm{\#}0`$, then the ground and the excited states will be strongly coupled through vibronic mixing and a large spin-orbit constant could be expected. This is the case for CCH which has a $`{}_{}{}^{2}\mathrm{\Sigma }`$ ground state and an A$`{}_{}{}^{2}\mathrm{\Pi }`$ state 3600 cm-1 above. The $`\nu _2`$ bending mode of the ground electronic state has A$`{}_{SO}{}^{}`$0.3 cm-1. For C4H, ASO is very large, $``$3 cm-1 (Yamamoto et al. 1987), because its A$`{}_{}{}^{2}\mathrm{\Pi }`$ state is only $``$468 cm-1 above the $`{}_{}{}^{2}\mathrm{\Sigma }`$ ground state (see above). For C4 the situation could be very similar to that of C4H as its first excited electronic state is $`{}_{}{}^{1}\mathrm{\Delta }`$ and is only 2680 cm-1 above the ground electronic $`{}_{}{}^{3}\mathrm{\Sigma }`$ state (Arnold et al. 1991; Xu et al. 1997). Last panel in figure 1 shows the computed band shape for ASO=4 cm-1 at different temperatures. The similarity between the observed and computed shape is appealing. Nevertheless, we note that the adopted spin-orbit constant is particulary large. We have run several models with different ASO values and conclude that the shape of the band is only crudely reproduced for A$`{}_{SO}{}^{}>`$2 cm-1. Another possibility that may explain the complexity of the band shape and allows smaller ASO values is that 13CCCC and C13CCC isotopomers also make a remarkable contribution to the observed emission. This could also apply to C4H if finally this molecule is confirmed as the carrier. We note, however, that in this case we should also expect a very large column density of C4H, which looks incompatible with the upper limits found for this molecule in Sgr B2 (Guélin, Friberg and Mezaoui 1982). If the carrier is C4 then the column density will be particularly large in all observed sources. Assuming a vibrational excitation temperature of 50-100 K in AGB stars and PPNs, N(C4) has to be close to 1015 cm-2, i.e., a typical abundance of a few 10-7 \- 10-6. In Sgr B2 the computed column density (assuming a lower vibrational excitation temperature) is 1.2 1015 cm-2, i.e., similar to that derived by Cernicharo, Goicoechea and Caux (2000) for C3. The accuracy of the band intensity from ab initio calculations dominates the error in the estimated column densities.

Future space–born platforms, such as the Herschel Space Observatory, will be equipped with heterodyne instruments with much better spectral resolution which should permit the detection of longer Cn chains through their low lying bending modes in the far–IR. The detection of these molecules and the determination of their abundances will allow to stablish their role in the formation of the carriers of the unidentified infrared bands. Waiting for Herschel, important work has to be done in spectroscopy laboratories to fully characterize these molecular species.

We thank Spanish DGES and PNIE for funding support under grant ESP98–1351–E and PANAYA2000-1784. JRG acknowledges Universidad Autónoma de Madrid for a pre–doctoral fellowship. We would like to thank M. Guélin and J.R. Pardo for their useful comments and suggestions.