text:
# Optimally Shaped Terahertz Pulses for Phase Retrieval in a Rydberg Atom Data Register

## I INTRODUCTION

The interaction of terahertz frequency electromagnetic pulses with Rydberg atoms has produced many insights into the dynamical properties of atomic systems . The comparable time scales of terahertz pulses with those of Rydberg state lifetimes make it possible to envision schemes of quantum control. In this paper, we propose a method of controlling Rydberg wave packets using shaped terahertz pulses, and theoretically show how these pulses can be designed to execute a quantum algorithm on a Rydberg atom data register.

It has been shown that information can be stored in the phases of the constituent orbitals of a Rydberg wave packet . Recently, a terahertz half-cycle pulse was used to decode the information stored in a Rydberg atom data register . This half-cycle pulse decodes the phase structure by retaining the population only in the orbital that was initially $`180^{}`$ out of phase with respect to the other orbitals, i. e. , the marked bit. However, this guess (unshaped) pulse does not decode all marked bits of the register with the same efficiency. In this paper, we aim to find the shaped THz pulse that will optimally transfer most of the population to any marked bit of the quantum data register.

## II OPTIMAL CONTROL THEORY

To design the terahertz pulse, we use a method that has been used extensively in mathematical and engineering applications — Optimal Control Theory (OCT) . This theory has also been applied to the control of quantum systems  with some success in experimental implementation . We use OCT to design a terahertz frequency pulse that can be used to achieve a desired target state. We then modify the OCT target state to make it possible to discover not a single target, but an optimized quantum algorithm. The wave function of the Rydberg electron is the state variable, and the electric field of the terahertz pulse is the control parameter. A functional $`J`$ is defined, whose extremum must be calculated. The functional consists of two parts, representing the desired target and the cost. Our aim is to maximize the fraction of the electron probability density in a target orbital $`|a_k`$ at a time T (after the end of the terahertz pulse). That is, the target functional $`P_k(T)=\psi (T)|a_ka_k|\psi (T)`$ must be a maximum. The cost functional represents the constraint on the control parameter, the terahertz field $`E(t)`$. The integrated energy of the pulse must be kept low, therefore the cost functional is defined as $`Y(T)=_0^T𝑑t\mathrm{}(t)|E(t)|^2`$. Here $`\mathrm{}`$ is a penalty parameter, in general time dependent, that controls the cost functional, and hence the peak terahertz field. The functional $`J`$ written as

$`J`$ $`=`$ $`\psi (T)|P_k|\psi (T)Y(T),`$ (1)
must be maximized.

The aim is to find an optimal control function $`E(t)`$ that maximizes $`J`$. The wave packet evolution is governed by the Schrödinger equation. In atomic units, $`e=m_e=\mathrm{}=1`$,

$`|\dot{\psi }(t)`$ $`=`$ $`iH(t)|\psi (t),`$ (2)
where $`H(t)=H_0+E(t)z`$. The equation of motion acts as a constraint on the evolution of the state, and in a manner similar to that used in variational calculus, we introduce a Lagrange multiplier $`|\lambda (t)`$ . The unconstrained functional that must be optimized is written as

$`\overline{J}`$ $`=`$ $`J{\displaystyle _0^T}dt[\lambda (t)|\dot{\psi }t+\lambda (t)|iH|\psi (t)+\psi (t)|\lambda (t)`$ (4)
$`i\dot{\psi }t|H|\lambda (t)].`$
$`=`$ $`\psi (T)|P_k|\psi (T){\displaystyle _0^T}𝑑t\mathrm{}(t)|E(t)|^22Re\lambda (t)|\psi (t)|_0^T`$ (6)
$`+{\displaystyle _0^T}𝑑t2Re\dot{\lambda (t)}|\psi (t){\displaystyle _0^T}𝑑t2Re\lambda (t)|iH|\psi (t).`$
Several iterative techniques for determining the optimal solution have been developed . Following the scheme for an iterative solution proposed in Ref. , the functional $`\overline{J}`$ is written as the sum of a terminal part and an integral

$`\overline{J}`$ $`=`$ $`G+{\displaystyle _0^T}𝑑tR,\mathrm{where}`$ (7)
$`G`$ $`=`$ $`\psi (T)|P_k|\psi (T)2Re\lambda (t)|\psi (t)|_0^T,`$ (8)
$`R`$ $`=`$ $`\mathrm{}(t)|E(t)|^2+2Re\left[\dot{\lambda (t)}|\psi (t)i\lambda (t)|H|\psi (t)\right].`$ (9)
The maximum of both G and R is sufficient to ensure the maximum of $`\overline{J}`$. Our objective is to iteratively determine the optimal function $`E(t)`$ that maximizes $`\overline{J}`$. The functional $`\overline{J}`$ at the $`k^{th}`$ and $`(k+1)^{st}`$ iteration is defined by the changes in the wave function and the field. That is,

$`|\psi (t)^{k+1}`$ $``$ $`|\psi (t)^k+|\mathrm{\Delta }\psi (t),`$ (10)
$`E^{k+1}`$ $``$ $`E^k+\mathrm{\Delta }E.`$ (11)
The difference in $`\overline{J}`$ between two successive iterations is written as

$`\overline{J}^{k+1}\overline{J}^k`$ $`=`$ $`\mathrm{\Delta }_1+\mathrm{\Delta }_2+\mathrm{\Delta }_3,\mathrm{where}`$ (12)
$`\mathrm{\Delta }_1`$ $``$ $`G(\psi ^{k+1}(T))G(\psi ^k(T))`$ (13)
$`=`$ $`2Re\left[\psi ^k(T)|P_k|\mathrm{\Delta }\psi (T)\lambda (T)|\mathrm{\Delta }\psi (T)\right]`$ (15)
$`+\mathrm{\Delta }\psi (T)|P_k|\mathrm{\Delta }\psi (T),`$
$`\mathrm{\Delta }_2`$ $``$ $`{\displaystyle _0^T}𝑑t[R(t,\psi ^{k+1},E^{k+1})R(t,\psi ^{k+1},E^k)]`$ (16)
$`=`$ $`2Re{\displaystyle _0^T}𝑑t\mathrm{}(t)E^k(t)\mathrm{\Delta }E(t)`$ (20)
$`{\displaystyle _0^T}𝑑t\mathrm{}(t)|\mathrm{\Delta }E(t)|^2`$
$`2Re{\displaystyle _0^T}dt[i\lambda ^k(t)|z\mathrm{\Delta }E(t)|\psi ^k(t)`$
$`+i\lambda ^k(t)|z\mathrm{\Delta }E(t)|\mathrm{\Delta }\psi (t)],`$
$`\mathrm{\Delta }_3`$ $``$ $`{\displaystyle _0^T}𝑑t[R(t,\psi ^{k+1},E^k)R(t,\psi ^k,E^k)]`$ (21)
$`=`$ $`2Re[{\displaystyle _0^T}dt\dot{\lambda (t)}|\mathrm{\Delta }\psi (t)`$ (23)
$`i{\displaystyle _0^T}dt\lambda (t)|H^k|\mathrm{\Delta }\psi (t)].`$
Choosing

$`|\lambda (T)`$ $`=`$ $`P_k|\psi ^k(T),\mathrm{and}`$ (24)
$`|\dot{\lambda (t)}`$ $`=`$ $`iH^k|\lambda (t),\mathrm{we}\mathrm{find}`$ (25)
$`\mathrm{\Delta }_1`$ $`=`$ $`\mathrm{\Delta }\psi (T)|P_k|\mathrm{\Delta }\psi (T).`$ (26)
$`\mathrm{\Delta }_2`$ $`=`$ $`2Re{\displaystyle _0^T}dt[\mathrm{}(t)E^k(t)\mathrm{\Delta }E(t)+{\displaystyle \frac{1}{2}}\mathrm{}(t)|\mathrm{\Delta }E(t)|^2`$ (29)
$`+i\lambda ^k(t)|z\mathrm{\Delta }E(t)|\psi ^k(t)+`$
$`i\lambda ^k(t)|z\mathrm{\Delta }E(t)|\mathrm{\Delta }\psi (t)].`$
$`\mathrm{\Delta }_3`$ $`=`$ $`0.`$ (30)
The solution improves when $`\overline{J}`$ increases or stays the same at each iteration. $`P_k`$ is a positive semidefinite operator, therefore $`\mathrm{\Delta }_1`$ is greater than or equal to zero. The change in the control parameter $`\mathrm{\Delta }E`$ is chosen to maximize $`\mathrm{\Delta }_2`$. The expression for $`\mathrm{\Delta }_2`$ suggests the following change in the field at the (k+1)th iteration :

$`\mathrm{\Delta }E(t)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}(t)}}\lambda ^k(t)|z|\psi ^{k+1}(t).`$ (31)
The appearance of $`|\psi ^{k+1}(t)`$ in the above expression implies that the overlap of $`|\psi (t)`$ and $`|\lambda (t)`$ is fed back immediately to find the field at the next time step.

The optimal control algorithm consists of the following steps:

1. Starting from the initial wave packet $`|\psi ^{(0)}(0)=|\psi (0)`$ and a first guess for the terahertz field $`E^{(0)}(t)`$, the wave packet is propagated according to Eq. 2 to find $`|\psi ^{(0)}(T)`$.
2. Using Eq. 24 to find $`|\lambda (T)`$, Eq. 25 is iterated backward to time $`t=0`$, and $`|\lambda (t)`$ is found at every time step.
3. With $`|\psi ^1(0)=|\psi (0)`$, equation 31 is then used to find a new value of the control field $`E^1(t)`$ and equation 2 is used to propagate the wavepacket forward in time.

The second two steps are repeated until the target yield converges to within the desired accuracy.

## III OPTIMAL PULSE FOR A SINGLE TARGET STATE

The first step in our approach is to find the optimal THz field required to decode a single flipped state. The initial state of the Rydberg data register is a wavepacket made of the $`24p`$ through $`29p`$ orbitals of equal amplitudes and the phase of the $`26p`$ (the marked bit) orbital opposite to that of the others. The initial guess THz pulse is a half-cycle pulse  of pulse width $`1ps`$. The desired initial phase structure occurs at the peak of the half-cycle pulse (i. e. at $`0.5`$ps). We find the THz pulse that will optimally transfer most of the population to the marked bit. The best value of the penalty parameter, $`\mathrm{}`$, that controls the peak field of the THz pulse to a reasonable value, and at the same time produces the desired final state, was found to be roughly $`10^{10}`$. One feature of several optimal fields obtained theoretically  is that the fields do not go to zero smoothly at the times $`t=0`$ and $`t=T`$. In the present system, the correct evolution of the wave packet depends very sensitively on the fields at the end points. The addition of a smooth switch-on and switch-off of the calculated optimal field drastically changes the evolution of the wave packet. Therefore, the condition that the THz field continuously goes to zero before and after the time interval of choice must be built into the algorithm. To ensure that the THz field goes smoothly to zero at times $`t=0`$ and $`t=T`$, the penalty parameter $`\mathrm{}(t)`$ is made a smoothly varying time-dependent function. The penalty on the pulse fluence is a thousand times more at the end points than at the rest of the pulse duration. The smoothness of the penalty function ensures a smooth switching on and switching off of the THz pulse.

This OCT implementation is very successful in describing the terahertz control of a Rydberg wave packet. The method takes macrosteps in the control field at every iteration and convergence is swift. The computational complexity is of the same order as that of the wave packet propagation. Therefore, we use a split-operator method in a restricted basis of essential states . The energy eigenstates of cesium are calculated using a pseudopotential method on a nonlinear radial grid . A time step of 10fs ensures the accuracy of the propagator, which is correct through the second-order in the time step. The numerical implementation of the local iterative algorithm is extremely sensitive to numerical error, and $`\mathrm{\Delta }_3`$ must be maintained equal to zero to very high precision . The unitary nature of the symmetrized-product propagator maintains this condition. The restricted basis consists of 195 energy eigenstates with principal quantum number $`n`$ ranging between 21 and 31, and angular momentum quantum number $`\mathrm{}<17`$. Absorbers ensure that population does not get reflected from the $`n=21`$, $`n=31`$, and $`\mathrm{}=16`$ ‘boundaries’. Using the selected state basis also imposes the condition that the spectrum of the THz pulse is decided by the energy range of the selected state basis. In this implementation of optimal control theory, we have chosen a fixed pulse length, $`T`$, of roughly 8ps. In other formalisms, this time $`T`$ may also be varied as a parameter.

The THz field that optimizes the population in the marked $`26p`$ state is shown in Fig. LABEL:fig1(a). The initial population in the $`26p`$ state is $`16.7\%`$. The optimal pulse will decode the information stored as phase by transferring most of the population into the $`26p`$ state. With the initial guess pulse, the population is $`29.5\%`$. After 50 iterations, the target yield is increased to $`52.8\%`$. The spectrum and Husimi distribution  of this optimal pulse are shown in Fig. LABEL:fig1(b) and Fig. LABEL:fig1(c) respectively. Notably, the strong peaks in the spectrum and the Husimi distribution do not correspond to any resonance between the energy levels of the selected state basis. The optimal terahertz pulse does not drive the system to any particular resonant condition. Instead, it alters the phases of the constituent orbitals of the wave packet so that they interfere to produce the desired probability distribution.

Figure LABEL:fig2 shows the evolution of the wavepacket as a function of time while the optimal pulse is on. During the THz pulse, probability can leak into other states not in the register (the other states in the essential basis). At the end of the pulse, a large fraction of the electron probability density lies in the flipped orbital (marked bit) of the data register. This can be thought of as using the other states of the data register as working qubits, which are used during the computation, but are not measured for any useful retrieval of information.

One interesting feature of this optimal pulse is that the peak field of roughly $`1\mathrm{K}\mathrm{V}/\mathrm{cm}`$ lasts for roughly $`0.5`$ps. For a $`\overline{n}=26`$ wavepacket, this field which is beyond the field ionization limit lasts for more than half the Kepler period ($`2\pi n^3`$). Yet, $`99\%`$ of the population remains in the selected state basis. This feature is an example of interferometric stabilization , seen in other atomic systems.

## IV OPTIMAL PULSE FOR A QUANTUM ALGORITHM

This THz pulse is optimal only for decoding the flipped $`26p`$ orbital. That is, if the phase of a different state were flipped, this pulse will not decode it. We wish to design a universal THz pulse that will optimally decode any flipped orbital of the wave packet register. Therefore we redefine the optimal control problem by considering an initial state that is a product state of independent wave packets with singly flipped orbitals.

$`|\mathrm{\Psi }(0)`$ $`=`$ $`|\psi _{25p}^{(1)}(0)\psi _{26p}^{(2)}(0)\psi _{27p}^{(3)}(0)\mathrm{}`$ (32)
The terahertz pulse acts simultaneously, but independently on all these wave packets. The desired final state is also a product state of independent wave packets with the flipped bit correctly decoded.

$`|\mathrm{\Psi }(T)`$ $`=`$ $`|25p^{(1)}|26p^{(2)}|27p^{(3)}\mathrm{}`$ (33)
The counterparts of Eq. 24 and Eq. 25 are straightforward. At every time step, the updated THz field is found by using a modified version of Eq. 31, with the matrix element of $`z`$ replaced by a sum of matrix elements of $`z`$, one from each independent ‘subspace’.

$`\mathrm{\Delta }E(t)`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}(t)}}{\displaystyle \underset{i=1}{\overset{N}{}}}\lambda _{(i)}^k(t)|iz|\psi _{(i)}^{k+1}(t).`$ (34)
Using this method, we find the terahertz pulse that detects any flipped orbital of the $`N`$-bit data register. The advantage of this refinement is that the computational resources needed increase only by a factor of the number of constituent states in the wave packet register.

We now find the optimal terahertz pulse that will decode any flipped state in a six state Rydberg data register. The register consists of $`np`$ states of cesium, with $`n`$ from $`2429`$. Population in the flipped orbital is amplified by the diffusion of probability density from the adjacent states. This is an example of the implementation of Grover’s search algorithm, where information is stored in states with differing phases, and a marked bit is amplified by “quantum diffusion” . The outer states $`n=24`$ and $`n=29`$ are therefore not included in the optimization. The universal decoding pulse and its effect on a wave packet with different marked bits is shown in Fig. LABEL:fig3. After the pulse, the wave packet population is distributed so that the flipped state is clearly amplified.

## V CONCLUSIONS

In conclusion, we have designed a terahertz pulse to implement a search algorithm on a quantum data register. Phase information stored in a Rydberg wave packet was optimally retrieved through the interaction with the pulse. Careful attention was paid to the smooth switch on and switch off of the THz pulse. We also show that it is possible to design an optimal pulse that can achieve not only a desired target state of an atom, but also implement a desired algorithm.

To our knowledge, this is the first time that optimal control theory has been applied to the terahertz control of a quantum system. This theoretical study motivates the experimental design and control of terahertz frequency pulses. Beyond quantum control, these results point to the possibilities of using Rydberg atoms as quantum computers, and terahertz pulses to implement quantum algorithms.

## ACKNOWLEDGEMENTS

It is a pleasure to acknowledge discussions and help from Professors D. J. Tannor, A. M. Bloch, V. F. Krotov and V. Malinovsky. C. R. acknowledges support from the Fellows and Visitor’s program of the Center for Ultrafast Optical Science at the University of Michigan. This research was supported by the National Science Foundation grant 9987916.


format: text

text:
# An Analytical and numerical investigation of escape rate for a noise driven bath

## I  Introduction

More than half a century ago Kramers proposed a diffusion model for chemical reactions in terms of the theory of Brownian motion in phase space . Since then the model and several of its variants have been ubiquitous in many areas of physics, chemistry and biology for understanding the nature of activated processes in classical , quantum and semiclassical  systems, in general. These have become the subject of several reviews  and monograph  in the recent past.

In the majority of these treatments one is essentially concerned with an equilibrium thermal bath at a finite temperature which stimulates the reaction co-ordinate to cross the activation energy barrier. The inherent noise of the medium is of internal origin. This implies that the dissipative force which the system experiences in course of its motion in the medium and the stochastic force acting on the system as a result of random impact from the constituents of the medium arise from a common mechanism. From a microscopic point of view the system-reservoir Hamiltonian description  developed over the decades suggests that the coupling of the system and the reservoir co-ordinates determines both the noise and the dissipative terms in the Langevin equation describing the motion of the system. It is therefore not difficult to anticipate that these two entities get related through a fluctuation-dissipation relation  ( these systems are sometimes classified as thermodynamically closed system in contrast to the systems driven by external noise in nonequilibrium statistical mechanics  ). However, when the reservoir is modulated by an external noise it is likely that this relation gets affected in a significant way. Since the modulation of the reservoir crucially depend on its response function, one can further envisage a connection between the dissipation of the system and the response function of the reservoir due to the external noise from a microscopic point of view.

In the present paper we explore this connection in the context of activated rate processes when the reservoir is modulated by an external noise. Specifically our object here is twofold: (i) to explore the role of reservoir response as a function of external noise on the system dynamics in contrast to direct driving of the system by external noise, (ii) to calculate the generalized Kramers’ rate for the steady state of this nonequilibrium open system taking full care of thermodynamic consistency condition. Both the internal and the external noises are Gaussian, stationary and are characterized by arbitrary decaying correlation functions. While the internal noise of the reservoir is thermal, the external noise may be of thermal or non-thermal type. We consider the stochastic motion to be spatial diffusion limited and calculate the rate of escape in the intermediate to strong damping regime. It is worth to mention that the externally generated nonequilibrium fluctuations can bias the Brownian motion of a particle in an anisotropic medium and may be used for designing molecular motors and pumps . We further mention that nonequilibrium, non-thermal systems have also been investigated phenomenologically by a number of workers in several other contexts, e.g., for examining the role of color noise in stationary probabilities , properties of nonlinear systems , nature of crossover , rate of diffusion limited coagulation processes , effect of monochromatic noise , etc. While these treatments concern direct driving of the system by an external noise the present consideration is based on modulation of the bath. A number of different situations depicting the modulation of the bath by an external noise may be physically relevant. As, for example, we consider a simple unimolecular conversion (say, an isomerization reaction) from $`AB`$. The reaction is carried out in a photochemically active solvent under the influence of external fluctuating light intensity. Since the fluctuations in the light intensity result in the fluctuations in the polarization of the solvent molecules, the effective reaction field around the reactant system gets modified. Provided the required stationarity of this nonequilibrium open system is maintained (which is not difficult in view of the experiments performed in the studies of external noise-induced transitions in photochemical systems ) the dynamics of barrier crossing becomes amenable to the present theoretical analysis that follows.

The remaining part of this paper is organized as follows: In Sec. II we discuss a system-reservoir model where the later is modulated by an external noise and establish an interesting connection between the dissipation of the system and the response function of the reservoir due to external noise. The stochastic motion in a linearized potential field is described in terms of a Fokker-Planck equation in Sec. III. Based on the traditional flux over population method  we derive in Sec. IV the generalized expression for the Kramers’ rate of escape from a metastable well. In Sec. V we numerically analyze the model and the bath modulated dynamics for the full potential and verify the theoretical rate with numerical simulation. The paper is concluded in Sec. VI.

## II  The system-reservoir model : The reservoir modulated by an external noise

We consider a classical particle of mass $`M`$ is linearly coupled to a heat bath of $`N`$ harmonic oscillators driven by an external noise. The total Hamiltonian is given by

$$H=\frac{p^2}{2M}+V(x)+\frac{1}{2}\underset{i=1}{\overset{N}{}}\left\{\frac{p_i^2}{m_i}+m_i\omega _i^2(q_ig_ix)^2\right\}+H_{int}.$$
(1)
In Eq.(1), $`x`$ and $`p`$ are the co-ordinate and momentum of the system particle; $`(q_i,p_i)`$ are the variables associated with the $`i`$-th oscillator and $`\omega _i`$ and $`m_i`$ are the corresponding frequency and mass, respectively. $`g_ix`$ measures the interaction between the particle and the bath. $`V(x)`$ is the potential energy of the particle. $`H_{int}`$ is assumed to be of the form

$$H_{int}=\frac{1}{2}\underset{i=1}{\overset{N}{}}\kappa _iq_iϵ(t).$$
(2)
The coupling function $`\kappa _i`$ measures the strength of interaction and $`ϵ(t)`$ is the external noise which we assume to be stationary and Gaussian with zero mean, i.e., $`ϵ(t)_e=0`$ and is characterized by an arbitrary correlation function as follows;

$`ϵ(t)ϵ(t^{})_e=2D\mathrm{\Psi }(tt^{}).`$
Here $`\mathrm{}_e`$ implies the averaging over the external noise.

We then eliminate the bath degrees of freedom in the usual way  to obtain the following generalized Langevin equation

$`\dot{x}`$ $`=`$ $`v,`$ (3)
$`\dot{v}`$ $`=`$ $`{\displaystyle \frac{dV}{dx}}{\displaystyle _0^t}𝑑t^{}\gamma (tt^{})v(t^{})+f(t)+\pi (t)`$ (4)
\[ while constructing Eq.(3) we have set $`M`$ and $`m_i`$ equal to unity \] where

$$\gamma (t)=\underset{i=1}{\overset{N}{}}g_i^2\omega _i^2\mathrm{cos}\omega _it.$$
(5)
$`f(t)`$ is the internal forcing function generated through the coupling between the system and the heat bath and is given by

$$f(t)=\underset{i=1}{\overset{N}{}}g_i\left\{\left[q_i(0)g_ix(0)\right]\omega _i^2\mathrm{cos}\omega _it+v_i(0)\omega _i\mathrm{sin}\omega _it\right\}.$$
(6)
In Eq.(3), $`\pi (t)`$ is a fluctuating force term due to the external noise $`ϵ(t)`$ and is given by

$$\pi (t)=_0^t\phi (tt^{})ϵ(t^{})𝑑t^{},$$
(7)
where

$$\phi (t)=\underset{i=1}{\overset{N}{}}g_i\kappa _i\omega _i\mathrm{sin}\omega _it.$$
(8)
The form of Eq.(3) therefore suggests that the system is driven by two forcing functions $`f(t)`$ and $`\pi (t)`$. $`f(t)`$ depends on the initial conditions of the bath oscillators for a fixed choice of the initial condition of the system degrees of freedom. To define the statistical properties of $`f(t)`$, we assume that the initial distribution is one in which the bath is equilibrated at $`t=0`$ in the presence of the system but in the absence of the external noise agency such that $`f(t)=0`$ and $`f(t)f(t^{})=k_BT\gamma (tt^{})`$.

Now, at $`t=0_+`$, the external noise agency is switched on and the bath is modulated by $`ϵ(t)`$. The system is governed by Eq.(3), where, apart from the internal noise $`f(t)`$, another fluctuating force $`\pi (t)`$ appears, that depends on the external noise $`ϵ(t)`$. Therefore, one can define an effective noise $`\xi (t)(=f(t)+\pi (t))`$ whose correlation is given by

$$\xi (t)\xi (t^{})=k_BT\gamma (tt^{})+2D_0^t𝑑t^{}_0^t^{}𝑑t^{\prime \prime }\phi (tt^{})\phi (t^{}t^{\prime \prime })\mathrm{\Psi }(t^{}t^{\prime \prime }),$$
(9)
along with $`\xi (t)=0`$, where $`\mathrm{}`$ means we have taken two averages independently. It should be emphasized that the above relation (9) is not a fluctuation-dissipation relation due to the appearance of the external noise intensity. Rather it serves as a thermodynamic consistency condition.

Let us now digress a little bit about $`\pi (t)`$. The statistical properties of $`\pi (t)`$ are determined by the normal mode density of the bath frequencies, the coupling of the system with the bath, the coupling of the bath with the external noise and the external noise itself. Eq.(7) is reminiscent of the familiar linear relation between the polarization and the external field where $`\pi `$ and $`ϵ`$ play the role of the former and the later, respectively. $`\phi (t)`$ can then be interpreted as a response function of the reservoir due to external noise $`ϵ(t)`$. The very structure of $`\pi (t)`$ suggests that this forcing function although originating from an external force is different from a direct driving force acting on the system. The distinction lies at the very nature of the bath characteristics (rather than system characteristics) as reflected in the relation (7) and (8).

With the coupling coefficients $`g(\omega )=g_0/\sqrt{\tau _c}\omega `$ and $`\kappa (\omega )=\sqrt{\tau _c}\omega \kappa _0`$, in the continuum limit  $`\gamma (t)`$ and $`\phi (t)`$ reduce to the following forms

$$\gamma (t)=\frac{g_0^2}{\tau _c}𝑑\omega 𝒟(\omega )\mathrm{cos}\omega t$$
(10)
and

$$\phi (t)=g_0\kappa _0𝑑\omega 𝒟(\omega )\omega \mathrm{sin}\omega t.$$
(11)
where $`g_0`$ and $`\kappa _0`$ are constants and $`\tau _c^1`$ is the cutoff frequency of the oscillator. $`𝒟(\omega )`$ is the density of modes of the heat bath.

From the above two relations, we obtain

$$\frac{d\gamma }{dt}=\frac{g_0}{\kappa _0}\frac{1}{\tau _c}\phi (t).$$
(12)
Eq.(12) is an important content of the present model. This expresses how the dissipative kernel $`\gamma (t)`$ depends on the response function $`\phi (t)`$ of the medium due to external noise $`ϵ(t)`$ \[ see Eq.(7) \]. Such a relation for the open system can be anticipated in view of the fact that both the dissipation and the response function crucially depend on the properties of the reservoir especially on its density of modes and its coupling to the system and the external noise source. In what follows we shall be concerned with the consequences of this relation in terms of the Langevin description in the next section ( Eq.(13) ) and numerical analysis of the full model potential in Sec. V.

## III  The Fokker-Planck description of the linearized motion : Asymptotic analysis of the Fokker-Planck coefficients

We now consider the system to be a harmonically bound particle of unit mass and of frequency $`\omega _0`$. Then because of Eq.(7) the Langevin equation (3) becomes

$`\dot{x}`$ $`=`$ $`v,`$ (13)
$`\dot{v}`$ $`=`$ $`\omega _0^2x{\displaystyle _0^t}𝑑t^{}\gamma (tt^{})v(t^{})+f(t){\displaystyle _0^t}𝑑t^{}\phi (tt^{})ϵ(t^{})`$ (14)
The Laplace transform of Eq. (13) allows us to write a formal solution for the displacement of the form

$`x(t)`$ $`=`$ $`x(t)+{\displaystyle _0^t}𝑑t^{}h(tt^{})f(t^{}){\displaystyle \frac{\kappa _0}{g_0}}\tau _c\omega _0^2{\displaystyle _0^t}𝑑t^{}h(tt^{})ϵ(t^{})`$ (16)
$`{\displaystyle \frac{\kappa _0}{g_0}}\tau _c{\displaystyle _0^t}𝑑t^{}h_2(tt^{})ϵ(t^{}),`$
where we have made use of the relation (12) explicitly.

Here

$$x(t)=\chi _x(t)x(0)+h(t)v(0)$$
(17)
with $`x(0)`$ and $`v(0)`$ being the initial position and initial velocity of the oscillator, respectively, which are nonrandom and

$$\chi _x(t)=\left[1\omega _0^2_0^th(\tau )𝑑\tau \right].$$
(18)
The kernel $`h(t)`$ is the Laplace inversion of

$$\stackrel{~}{h}(s)=\frac{1}{s^2+\stackrel{~}{\gamma }(s)s+\omega _0^2}$$
(19)
where, $`\stackrel{~}{\gamma }(s)=_0^{\mathrm{}}e^{st}\gamma (t)𝑑t`$, is the Laplace transform of the friction kernel $`\gamma (t)`$, and

$$h_2(t)=\frac{d^2h(t)}{dt^2}.$$
(20)
The time derivative of Eq.(16) yields

$`v(t)`$ $`=`$ $`v(t)+{\displaystyle _0^t}𝑑t^{}h_1(tt^{})f(t^{}){\displaystyle \frac{\kappa _0}{g_0}}\tau _c\omega _0^2{\displaystyle _0^t}𝑑t^{}h_1(tt^{})ϵ(t^{})`$ (22)
$`{\displaystyle \frac{\kappa _0}{g_0}}\tau _c{\displaystyle _0^t}𝑑t^{}h_3(tt^{})ϵ(t^{})`$
where

$$v(t)=\omega _0^2h(t)+v(0)h_1(t),$$
(23)
$$h_1(t)=\frac{dh(t)}{dt}\mathrm{and}h_3(t)=\frac{d^3h(t)}{dt^3}.$$
(24)
Next we calculate the variances. From the formal solution of $`x(t)`$ and $`v(t)`$, the explicit expressions for the variances are obtained which are given below;

$`\sigma _{xx}^2(t)`$ $`=`$ $`[x(t)x(t)]^2`$ (25)
$`=`$ $`2{\displaystyle _0^t}𝑑t_1h(t_1){\displaystyle _0^{t_1}}𝑑t_2h(t_2)f(t_1)f(t_2)`$ (29)
$`+2\left({\displaystyle \frac{\kappa _0}{g_0}}\tau _c\omega _0^2\right)^2{\displaystyle _0^t}𝑑t_1h(t_1){\displaystyle _0^{t_1}}𝑑t_2h(t_2)ϵ(t_1)ϵ(t_2)_e`$
$`+2\left({\displaystyle \frac{\kappa _0}{g_0}}\tau _c\right)^2{\displaystyle _0^t}𝑑t_1h_2(t_1){\displaystyle _0^{t_1}}𝑑t_2h_2(t_2)ϵ(t_1)ϵ(t_2)_e`$
$`+2\left({\displaystyle \frac{\kappa _0}{g_0}}\tau _c\right)^2\omega _0^2{\displaystyle _0^t}𝑑t_1h(t_1){\displaystyle _0^{t_1}}𝑑t_2h_2(t_2)ϵ(t_1)ϵ(t_2)_e,`$
$`\sigma _{vv}^2(t)`$ $`=`$ $`[v(t)v(t)]^2`$ (30)
$`=`$ $`2{\displaystyle _0^t}𝑑t_1h_1(t_1){\displaystyle _0^{t_1}}𝑑t_2h_1(t_2)f(t_1)f(t_2)`$ (34)
$`+2\left({\displaystyle \frac{\kappa _0}{g_0}}\tau _c\omega _0^2\right)^2{\displaystyle _0^t}𝑑t_1h_1(t_1){\displaystyle _0^{t_1}}𝑑t_2h_1(t_2)ϵ(t_1)ϵ(t_2)_e`$
$`+2\left({\displaystyle \frac{\kappa _0}{g_0}}\tau _c\right)^2{\displaystyle _0^t}𝑑t_1h_3(t_1){\displaystyle _0^{t_1}}𝑑t_2h_3(t_2)ϵ(t_1)ϵ(t_2)_e`$
$`+2\left({\displaystyle \frac{\kappa _0}{g_0}}\tau _c\right)^2\omega _0^2{\displaystyle _0^t}𝑑t_1h_1(t_1){\displaystyle _0^{t_1}}𝑑t_2h_3(t_2)ϵ(t_1)ϵ(t_2)_e,`$
and

$`\sigma _{xv}^2(t)`$ $`=`$ $`[x(t)x(t)][v(t)v(t)]`$ (35)
$`=`$ $`{\displaystyle \frac{1}{2}}\dot{\sigma }_{xx}^2(t)`$ (36)
where we have assumed that the noises $`f(t)`$ and $`ϵ(t)`$ are symmetric with respect to the time argument and have made use the fact that $`f(t)`$ and $`ϵ(t)`$ are uncorrelated.

Due to the Gaussian property of the noises $`f(t)`$ and $`ϵ(t)`$ and the linearity of the Langevin equation (13), we see that the joint probability density $`p(x,v,t)`$ of the oscillator must be Gaussian. The joint characteristic function associated with the density is

$$\stackrel{~}{p}(\mu ,\rho ,t)=\mathrm{exp}\left\{ix(t)\mu +iv(t)\rho \frac{1}{2}\left[\sigma _{xx}^2(t)\mu ^2+2\sigma _{xv}^2(t)\rho \mu +\sigma _{vv}^2(t)\rho ^2\right]\right\}.$$
(37)
Using the method of characteristic function  and the above expression (37) we find the general Fokker-Planck equation associated with the probability density function $`p(x,v,t)`$ for the process (13);

$$\frac{p}{t}=v\frac{p}{x}+\overline{\omega }_0^2(t)x\frac{p}{v}+\overline{\gamma }(t)\frac{}{v}(vp)+\varphi (t)\frac{^2p}{v^2}+\psi (t)\frac{^2p}{vx},$$
(38)
where

$`\overline{\gamma }(t)`$ $`=`$ $`{\displaystyle \frac{d}{dt}}\mathrm{ln}\mathrm{{\rm Y}}(t),\overline{\omega }_0^2(t)={\displaystyle \frac{h(t)h_1(t)+h_1^2(t)}{\mathrm{{\rm Y}}(t)}}\mathrm{and}`$
$`\mathrm{{\rm Y}}(t)`$ $`=`$ $`{\displaystyle \frac{h_1(t)}{\omega _0^2}}\left[1\omega _0^2{\displaystyle _0^\tau }𝑑\tau h(\tau )\right]+h^2(t).`$
The functions $`\varphi (t)`$ and $`\psi (t)`$ are defined by

$$\varphi (t)=\overline{\omega }_0^2(t)\sigma _{xv}^2+\overline{\gamma }\sigma _{vv}^2+\frac{1}{2}\dot{\sigma }_{xv}^2\mathrm{and}\psi (t)=\dot{\sigma }_{xv}^2+\overline{\gamma }(t)\sigma _{xv}^2+\overline{\omega }_0^2\sigma _{xx}^2\sigma _{vv}^2$$
(39)
where the covariances are to be calculated for a particular given noise process.

For the internal noise processes it had been shown earlier that for several models the various time dependent parameters $`\overline{\omega }_0^2(t)`$, $`\overline{\gamma }(t)`$, etc. do exist asymptotically as $`t\mathrm{}`$. The above consideration shows that $`h(t)`$, $`h_1(t)`$, etc. do not depend on the nature of the noise but depend only on the relaxation $`\overline{\gamma }(t)`$.

We now discuss the asymptotic properties of $`\varphi (t)`$ and $`\psi (t)`$, which in turn are dependent on the variances $`\sigma _{xx}^2(t)`$ and $`\sigma _{vv}^2(t)`$, as $`t\mathrm{}`$ since they play a significant role in our further analysis that follows.

From Eqs. (25) and (30), we may write

$`\sigma _{xx}^2(t)=\sigma _{xx}^{2(i)}(t)+\sigma _{xx}^{2(e)}(t)\mathrm{and}\sigma _{vv}^2(t)=\sigma _{vv}^{2(i)}(t)+\sigma _{vv}^{2(e)}(t).`$
where ‘$`i`$’ denotes the part corresponding to internal noise $`f(t)`$ and ‘$`e`$’ corresponds to the external noise $`ϵ(t)`$. Since the average velocity of the oscillator is zero as $`t\mathrm{}`$ we see from Eq.(23) that $`h(t)`$ and $`h_1(t)`$ must be zero as $`t\mathrm{}`$. Also from Eq.(17) we observe that the function $`\chi _x(t)`$ must decay to zero for long times. Hence, from Eq.(18) we see that the stationary value of the integral of $`h(t)`$ is $`1/\omega _0^2`$, i.e.,

$$_0^{\mathrm{}}h(t)𝑑t=\frac{1}{\omega _0^2}.$$
(40)
Now, $`\sigma _{xx}^{2(i)}(t)`$ and $`\sigma _{vv}^{2(i)}(t)`$ of Eqs.(25) and (30) can be written in the form

$`\sigma _{xx}^{2(i)}(t)`$ $`=`$ $`2{\displaystyle _0^t}𝑑t_1h(t_1){\displaystyle _0^{t_1}}𝑑t_2h(t_2)f(t_1)f(t_2)`$ (41)
$`=`$ $`k_BT\left[2{\displaystyle _0^t}𝑑\tau h(\tau )h^2(t)\omega _0^2\left\{{\displaystyle _0^t}𝑑\tau h(\tau )\right\}^2\right]`$ (42)
and

$`\sigma _{vv}^{2(i)}(t)`$ $`=`$ $`2{\displaystyle _0^t}𝑑t_1h_1(t_1){\displaystyle _0^{t_1}}𝑑t_2h_1(t_2)f(t_1)f(t_2)`$ (43)
$`=`$ $`k_BT\left[1h_1^2(t)\omega _0^2h^2(t)\right].`$ (44)
From the above two expressions \[ Eqs.(42) and (44) \] we see that

$$\sigma _{xx}^{2(i)}(\mathrm{})=\frac{k_BT}{\omega _0^2}\mathrm{and}\sigma _{vv}^{2(i)}(\mathrm{})=k_BT.$$
(45)
It is important to note that these stationary values are not related to the intensity and correlation time of the internal noise.

We next consider the parts, $`\sigma _{xx}^{2(e)}(t)`$ and $`\sigma _{vv}^{2(e)}(t)`$, due to the presence of the external noise. The Laplace transform of Eq.(16) yields the expression

$$\stackrel{~}{x}(s)\stackrel{~}{x}(s)=\stackrel{~}{h}(s)\stackrel{~}{f}(s)\frac{\kappa _0}{g_0}\tau _c\omega _0^2\stackrel{~}{h}(s)\stackrel{~}{ϵ}(s)\frac{\kappa _0}{g_0}\tau _cs^2\stackrel{~}{h}(s)\stackrel{~}{ϵ}(s)$$
(46)
where

$`\stackrel{~}{x}(s)`$ $`=`$ $`\left\{{\displaystyle \frac{1}{s}}{\displaystyle \frac{\omega _0^2}{s[s^2+s\stackrel{~}{\gamma }(s)+\omega _0^2]}}\right\}x(0)+{\displaystyle \frac{1}{s^2+s\stackrel{~}{\gamma }(s)+\omega _0^2}}v(0)`$ (47)
$`=`$ $`\left\{{\displaystyle \frac{1}{s}}\omega _0^2{\displaystyle \frac{\stackrel{~}{h}(s)}{s}}\right\}x(0)+\stackrel{~}{h}(s)v(0).`$ (48)
From the above equation (46) we can calculate the variance $`\sigma _{xx}^2`$ in the Laplace-transformed space which can be identified as the Laplace transform of Eq.(25). Thus, for the part $`\sigma _{xx}^{2(e)}(t)`$ we observe that, $`\stackrel{~}{\sigma }_{xx}^{2(e)}(s)`$ contains terms like $`\left(\frac{\kappa _0}{g_0}\tau _c\omega _0^2\stackrel{~}{h}(s)\right)^2\stackrel{~}{ϵ}^2(s)_e`$. Since, we have assumed the stationarity of the noise $`ϵ(t)`$, we conclude that if $`\stackrel{~}{C}(0)`$ exists \[ where $`C(tt^{})=ϵ(t)ϵ(t^{})_e`$ \], then the stationary value of $`\sigma _{xx}^{2(e)}(t)`$ exists and becomes a constant that depends on the correlation time and the strength of the noise. Similar argument is also valid for $`\sigma _{vv}^{2(e)}(t)`$.

Summarizing the above discussions we note that,

(i) the internal noise-driven parts of $`\sigma _{xx}^2(t)`$ and $`\sigma _{vv}^2(t)`$, i.e., $`\sigma _{xx}^{2(i)}`$ and $`\sigma _{vv}^{2(i)}`$ approach the fixed values which are independent of the noise correlation and the intensity as $`t\mathrm{}`$,

(ii) the external noise driven parts of variances also approach the constant values at the stationary ($`t\mathrm{}`$) limit which are dependent on the strength and the correlation time of the noise.

Hence we conclude, following the Ref.(33) and our preceding discussions that even in presence of an external noise the above terms do exist asymptotically and we write the steady state Fokker-Planck equation for the asymptotic values of the parameters as,

$$v\frac{p}{x}+\overline{\omega }_0^2x\frac{p}{v}+\overline{\gamma }\frac{}{v}(vp)+\varphi (\mathrm{})\frac{^2p}{v^2}+\psi (\mathrm{})\frac{^2p}{vx}=0,$$
(49)
where, $`\overline{\omega }_0^2`$, $`\overline{\gamma }`$, $`\varphi (\mathrm{})`$, $`\psi (\mathrm{})`$, etc. are to be calculated from the general definition (39) for the steady state. As an explicit example we consider the case of a $`\delta `$-correlated external noise and Ornstein-Uhlenbeck internal noise for which we provide the expressions for variances $`\sigma _{xx}^2(t)`$, $`\sigma _{vv}^2(t)`$ and $`\sigma _{xv}^2(t)`$ and the relaxation function $`h(t)`$ given in the Appendix-A.

The general steady state solution of the above equation (49) is

$$p_{st}(x,v)=\frac{1}{Z}\mathrm{exp}\left[\left\{\frac{v^2}{2D_0}+\frac{\overline{\omega }_0^2x}{2(D_0+\psi (\mathrm{}))}\right\}\right]$$
(50)
where

$$D_0=\frac{\varphi (\mathrm{})}{\overline{\gamma }}$$
(51)
and $`Z`$ is the normalization constant. The solution (50) can be verified by direct substitution. The distribution (50) is not an equilibrium distribution. This stationary distribution for the nonequilibrium open system plays the role of an equilibrium distribution of the closed system which may, however, be recovered in the absence of external noise term.

## IV  Kramers’ escape rate

We now turn to the problem of decay of a metastable state. In Kramers approach , the particle coordinate $`x`$ corresponds to the reaction coordinate and its values at the minima of the potential $`V(x)`$ denotes the reactant and product states.

Linearizing the motion around barrier top at $`x=x_b`$ the Langevin equation (3) can be written down as

$`\dot{y}`$ $`=`$ $`v,`$ (52)
$`\dot{v}`$ $`=`$ $`\omega _b^2y{\displaystyle _0^t}𝑑t^{}\gamma (tt^{})v(t^{})+f(t)+\pi (t),`$ (53)
where, $`y=xx_b`$ and the barrier frequency $`\omega _b^2`$ is defined by

$$V(y)=V_b\frac{1}{2}\omega _b^2y^2;\omega _b^2>0.$$
(54)
Correspondingly the motion of the particle is governed by the Fokker-Planck equation (38)

$$\frac{p}{t}=v\frac{p}{y}\overline{\omega }_b^2(t)y\frac{p}{v}+\overline{\gamma }_b(t)\frac{}{v}(vp)+\varphi _b(t)\frac{^2p}{v^2}+\psi _b(t)\frac{^2p}{vy},$$
(55)
where, the suffix ‘$`b`$’ indicates that all the coefficients are to be calculated using the general definition (39) for the barrier top region.

It is apparent from Eqs.(49) and (55) that since the dynamics is non-Markovian and the system is thermodynamically open one has to deal with the renormalized frequencies $`\overline{\omega }_0`$ and $`\overline{\omega }_b`$ near the bottom or top of the well, respectively. We make the ansatz that the nonequilibrium, steady state probability $`p_b`$, generating a nonvanishing diffusion current $`j`$, across the barrier is given by

$$p_b(x,v)=\mathrm{exp}\left[\left\{\frac{v^2}{2D_b}+\frac{\stackrel{~}{V}(x)}{D_b+\psi _b(\mathrm{})}\right\}\right]\xi (x,v)$$
(56)
where

$$D_b=\frac{\varphi _b(\mathrm{})}{\overline{\gamma }_b}.$$
(57)
$`\stackrel{~}{V}(x)`$ is the renormalized linear potential as

$`\stackrel{~}{V}(x)`$ $`=`$ $`V(x_0)+{\displaystyle \frac{1}{2}}\overline{\omega }_0^2(xx_0)^2,\mathrm{near}\mathrm{the}\mathrm{bottom}`$ (58)
$`\stackrel{~}{V}(x)`$ $`=`$ $`V(x_b){\displaystyle \frac{1}{2}}\overline{\omega }_b^2(xx_b)^2,\mathrm{near}\mathrm{the}\mathrm{top}`$ (59)
with $`\overline{\omega }_0^2`$, $`\overline{\omega }_b^2>0`$. The unknown function $`\xi (x,v)`$ obeys the natural boundary condition that for $`x\mathrm{}`$, $`\xi (x,v)`$ vanishes.

The ansatz of the form (56) denoting the steady state distribution is motivated by the local analysis near the bottom and the top of the barrier in the Kramers’ sense . For a stationary nonequilibrium system, on the other hand, the relative population of the two regions, in general, depends on the global properties of the potential leading to an additional factor in the rate expression. Although because of the Kramers’ type ansatz which is valid for the local analysis, such a consideration is outside the scope of the present treatment, we point out a distinctive feature in the ansatz (56) compared to Kramers’ ansatz. While in the latter case one considers a complete factorization of the equilibrium part (Boltzmann) and the dynamical part, the ansatz (56) incorporates the additional dynamical contribution through dissipation and strength of the noise into the exponential part. This modification of Kramers’ ansatz (by dynamics) is due to nonequilibrium nature of the system. Thus unlike Kramers’, the exponential factors in (56) and in the stationary distribution (50) which serves as a boundary condition are markedly different. Before carrying out global analysis in the present section our aim here is to understand the modification of the rate due to modulation of the bath driven by an external noise, within the perview of Kramers’ type ansatz. The internal consistency of the treatment, however, can be checked by recovering the Kramers’ result when the external noise is switched off.

From equation (55), using (56) we obtain the equation for $`\xi (y,v)`$ in the steady state in the neighborhood of $`x_b`$, the equation

$$\left(1+\frac{\psi _b(\mathrm{})}{D_b}\right)v\frac{\xi }{y}\left[\frac{D_b}{D_b+\psi _b(\mathrm{})}\overline{\omega }_b^2y+\overline{\gamma }_bv\right]\frac{\xi }{v}+\varphi _b(\mathrm{})\frac{^2\xi }{v^2}+\psi _b(\mathrm{})\frac{^2\xi }{vy}=0.$$
(60)
After making use of the appropriate transformations and boundary conditions for reduced distribution functions we obtain the barrier crossing rate $`k`$ given by

$$k=\frac{\overline{\omega }_0}{2\pi }\frac{D_b}{\{D_0+\psi (\mathrm{})\}^{1/2}}\left(\frac{\mathrm{\Lambda }}{1+\mathrm{\Lambda }D_b}\right)^{1/2}\mathrm{exp}\left[\frac{E_0}{D_b+\psi _b(\mathrm{})}\right]$$
(61)
where

$`\mathrm{\Lambda }={\displaystyle \frac{\lambda }{\varphi _b(\mathrm{})+a\psi _b(\mathrm{})}}`$
with

$`a={\displaystyle \frac{D_b}{2(D_b+\psi _b(\mathrm{}))}}\left\{\overline{\gamma }_b\sqrt{\overline{\gamma }_b^2+4\overline{\omega }_b^2}\right\}\mathrm{and}\lambda =\overline{\gamma }_ba\left(1+{\displaystyle \frac{\psi _b(\mathrm{})}{D_b}}\right).`$
Here $`E_0`$ is the activation energy, $`E_0=V(x_b)V(x_0)`$. Since the temperature due to internal thermal noise, the strength of the external noise and the damping constant are buried in the parameters $`D_0`$, $`D_b`$, $`\psi _0`$, $`\psi _0`$ and $`\mathrm{\Lambda }`$ the generalized expression look somewhat cumbersome. We point out that the subscripts ‘$`0`$’ and ‘$`b`$’ in $`D`$ and $`\psi `$ refer to the well or barrier top region, respectively. Eq.(61) is one of the key results of this paper. We note here that $`(D_b+\psi _b(\mathrm{}))/k_B`$ in the exponential factor defines a new effective temperature characteristic of the steady state of the nonequilibrium open system and an effective transmission factor is contained in the prefactor controlling the barrier crossing dynamics. As expected both are the functions of the external noise strength $`D`$ and coupling of noise to the bath modes. The dependence of the rate on the parameters can be exposed explicitly once we consider the typical cases.

## V  Numerical analysis of the model and the barrier crossing rate

### A  The model : bath modulation vs. direct driving

We first explore the distinctive aspects of the external noise-driven-reservoir-modulated dynamics of the system ( in contact with the external noise driven bath ) in contrast to direct driving of the system by the external noise. This will help us to elucidate the special role of the reservoir response function in controlling the rate. To this end we first return to our basic equation of the model, Eq.(3) where $`\phi (tt^{})`$ as given by Eq.(8) characterizes the response function. We solve this equation numerically using second order stochastic algorithm of Fox  for the full potential $`V(x)=\frac{1}{4}x^4\frac{1}{2}x^2`$ and plot the results of computation of the inverse of the mean first passage time as a function of the damping constant $`g_0^2`$ ( see Eq.(10) ) in Fig.(1). The density of modes of the bath has been assumed to be of the standard Lorentzian form, $`𝒟(\omega )=\frac{2}{\pi }\frac{1}{\tau _c}\left(\frac{1}{\omega ^2+\tau _c^2}\right)`$ with vanishing internal correlation time $`\tau _c`$. Typically for the other scaled parameters we have used $`k_BT`$ = 0.1, the strength of the $`\delta `$-correlated external noise, $`D`$ = 0.1. The mean first passage times have been calculated by averaging over 10000 trajectories. We first show the continuous curve in Fig.(1) which represents the case of direct driving by the external $`\delta `$-correlated noise. What is immediately apparent is the loss of turnover of the rate as one allows the variation of the dissipation constant $`g_0^2`$ from low to high friction regime. This is in sharp contrast to what is shown by the dotted curve in Fig.(1) which depicts the situation of bath modulation with the response function $`\phi (tt^{})`$ determined by Eq.(11). Both the curves ( continuous and dotted ) representing the open system are compared to that for the standard Kramers’ turnover, i.e., when the dynamics is calculated in absence of the external noise ( dashed curve ). As expected the rate in the driven system ( be it directly or through bath modulation ) is always higher that in the undriven one. Fig.(1) also shows that although in the low damping region the direct driving causes a much higher rate, it is, however, the bath modulation which becomes more effective in inducing activated barrier crossing in the high friction regime. It is thus interesting to note that the barrier crossing dynamics of the system in contact with an external noise modulated bath captures the essential turnover features of the Kramers’ dynamics of the closed system. We therefore realize that although open, the nature of the response function of the reservoir as well as the thermodynamic consistency condition (9), make the open system feel like a closed system.

In Fig.(2) we show the variation of inverse of the mean first passage time as a function of the strength of the external noise $`D`$, keeping all other parameter same as before but for $`g_0^2`$ = 2.0. It is apparent that the barrier crossing is more facilitated by modulating the bath than driving the system directly for higher values of external noise strength.

### B The rate : Theoretical results vs. numerical simulation

So far we have considered the full potential $`V(x)=\frac{1}{4}x^4\frac{1}{2}x^2`$ and vanishingly small correlation times for external and internal noises for numerical computation. We now turn to our basic theoretical result Eq.(61) which is a generalization of Kramers’ rate for bath modulated dynamics for intermediate to strong damping regime. To examine its validity we calculate the rate (61) as a function of the damping constant $`g_0^2`$ ( $`g_0`$ is related to both $`\gamma (t)`$ and $`\phi (t)`$ in Eqs.(10) and (11) respectively ) for several values of external noise strength $`D`$. The scaled barrier height $`E_0`$ and $`k_BT`$ have been set to 0.25 and 0.1, respectively. The results are plotted in Fig.(3) ( continuous curve ) and compared to the rate, inverse of the mean first passage time ( dotted curve ) calculated numerically using full potential for the set of parameter values ( as given in the earlier subsection ) with Eq.(3). In Figs.(4) and (5) we further compare our theoretical results ( continuous curve ) with numerical simulation ( dotted curve ) for the variation of rate as a function of external noise strength and its correlation time, respectively for several values of damping constant. It is apparent from Figs.(3-5) that the theoretical and numerical results are in good agreement.

## VI  Conclusions

Based on a system-reservoir microscopic model where the reservoir is modulated by an external, stationary and Gaussian noise with arbitrary decaying correlation function, we have numerically analyzed the model and generalized the Kramers’ theory to calculate the steady state rate of escape from a metastable well. The main conclusions of this study are as follows;

(i) We have shown that since the reservoir is driven by the external noise and the dissipative properties of the system depend on the reservoir, a simple connection between the dissipation and the response function of the medium due to the external noise can be established. This connection is important for realizing an effective temperature of the reservoir characterizing the stationary state of the thermodynamically open system, as well as an effective transmission factor controlling the rate. Both of these quantities depend on the strength and correlation time of the external noise.

(ii) Many of the earlier treatments of the rate concern direct (phenomenological) driving of the system and did not emphasize the question of thermodynamic consistency in the context of open systems. The present theory being microscopic the fluctuation-dissipation-like relation (9) remains an inbuilt characteristic of the model itself as an essential thermodynamic consistency condition.

(iii) Based on numerical simulation of the full model potential we show that one can recover the turnover features of the Kramers’ dynamics when the external noise modulates the reservoir rather than the system directly. This recovery, we believe, is an offshoot of the derived thermodynamic consistency condition (9).

(iv) Provided the long time limit of the moments for the stochastic processes pertaining to the external and internal noises characterized by arbitrary decaying correlation functions exist, the expression for generalized Kramers’ rate of barrier crossing for the open system we derive here is fairly general. We have shown that it agrees reasonably well with numerical simulation using full potential for this model.

The creation of a typical nonequilibrium open situation by modulating a bath with the help of an external noise is not an uncommon phenomenon in applications and industrial processing. The external agency generating noise does work on the bath by stirring, pumping, agitating, etc., to which the system dissipates internally. In the present treatment we are concerned with a nonequilibrium steady state which signifies a constant throughput of energy in contrast to thermal equilibrium defined by an constant temperature. We believe that these considerations are likely to be important in other related issues in nonequilibrium open systems and may serve as a basis for studying processes occurring within irreversibly driven environments  and for thermal ratchet problems .

###### Acknowledgements.

We thank E. Pollak for his kind interest in this work and for pointing out the Ref. (35) to us. SKB in indebted to Council of Scientific and Industrial Research (CSIR), Government of India for financial support.

## A Calculation of variances

We consider a particular case as an example where the external noise $`ϵ(t)`$ is $`\delta `$-correlated and the internal noise is an Ornstein-Uhlenbeck process, i.e.,

$`ϵ(t)ϵ(t^{})_e=2D\delta (tt^{})`$
and

$`f(t)f(t^{})={\displaystyle \frac{g_0^2k_BT}{\tau _c}}e^{|tt^{}|/\tau _c}.`$
Consequently, from the fluctuation-dissipation relation we derive the dissipative kernel as,

$`\gamma (tt^{})={\displaystyle \frac{g_0^2}{\tau _c}}e^{|tt^{}|/\tau _c}.`$
It should be noted that for $`\tau _c0`$, the above noise process become $`\delta `$-correlated.

The Laplace transform of $`\gamma (t)`$ as given above can be written as

$`\stackrel{~}{\gamma }(s)={\displaystyle \frac{g_0^2}{s\tau _c+1}},`$
and subsequently, we have for $`\tau _c0`$,

$`\stackrel{~}{h}(s)={\displaystyle \frac{s+a}{s^3+as^2+bs+c_0}}`$
where

$`a={\displaystyle \frac{1}{\tau _c}},b=\omega _0^2+{\displaystyle \frac{g_0^2}{\tau _c}}\mathrm{and}c_0={\displaystyle \frac{\omega _0^2}{\tau _c}}.`$
We find that the inverse Laplace transform of $`\stackrel{~}{h}(s)`$ reads

$$h(t)=c_1e^{\mathrm{\Delta }_1t}+c_2e^{\mathrm{\Delta }_2t}\mathrm{sin}(\beta t+\alpha )$$
(A1)
where the coefficients $`c_1`$, $`c_2`$, $`\mathrm{\Delta }_1`$, $`\mathrm{\Delta }_2`$, $`\beta `$ and $`\alpha `$ are given by

$`\mathrm{\Delta }_1`$ $`=`$ $`𝒜+{\displaystyle \frac{a}{3}},`$ (A3)
$`\mathrm{\Delta }_2`$ $`=`$ $`{\displaystyle \frac{1}{2}}(𝒜+)+{\displaystyle \frac{a}{3}},`$ (A4)
$`\beta `$ $`=`$ $`{\displaystyle \frac{\sqrt{3}}{2}}(𝒜),`$ (A5)
$`c_1`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{\Delta }_2\mathrm{\Delta }_1d}},`$ (A6)
$`d`$ $`=`$ $`{\displaystyle \frac{a(2\mathrm{\Delta }_2\mathrm{\Delta }_1)\mathrm{\Delta }_2^2\beta ^2}{a\mathrm{\Delta }_1}},`$ (A7)
$`𝒜`$ $`=`$ $`\left({\displaystyle \frac{a^3}{27}}+{\displaystyle \frac{ab}{6}}{\displaystyle \frac{c_0}{2}}+\sqrt{𝒬}\right)^{1/3},`$ (A8)
$``$ $`=`$ $`\left({\displaystyle \frac{a^3}{27}}+{\displaystyle \frac{ab}{6}}{\displaystyle \frac{c_0}{2}}\sqrt{𝒬}\right)^{1/3},`$ (A9)
$`c_2`$ $`=`$ $`{\displaystyle \frac{c_1}{\beta }}[(d\mathrm{\Delta }_2)^2+\beta ^2]^{1/2},`$ (A10)
$`\alpha `$ $`=`$ $`\mathrm{tan}^1\left({\displaystyle \frac{\beta }{d\mathrm{\Delta }_2}}\right)\mathrm{and}`$ (A11)
$`𝒬`$ $`=`$ $`{\displaystyle \frac{a^2b^2}{108}}+{\displaystyle \frac{b^3}{27}}+{\displaystyle \frac{a^3c_0}{27}}{\displaystyle \frac{abc_0}{6}}+{\displaystyle \frac{c_0^2}{4}}.`$ (A12)
Here we note that for a physically allowed solution $`\mathrm{\Delta }_1`$, $`\mathrm{\Delta }_2`$ must be positive. Since by Eq.(19) $`h(t)`$ depends on the memory kernel $`\gamma (t)`$ which is of decaying type and all the moments, in general, reach asymptotic constancy as shown in Sec. III, these quantities are positive (which depends on the correlation time $`\tau _c`$, the strength of the noise and other potential parameters) which may be checked (after some algebra) by considering the limiting cases such as $`\tau _c0`$ and $`\tau _c`$ large.

Substituting Eq.(A1) into the expressions for variances \[external noise is $`\delta `$-correlated\], namely into (25) and (30) we have after some lengthy algebra

$`\sigma _{xx}^2(t)=\sigma _{xx}^{2(i)}(t)+\sigma _{xx}^{2(e)}(t)`$
where

$`\sigma _{xx}^{2(i)}(t)`$ $`=`$ $`k_BT\left(c_2R+{\displaystyle \frac{c_1}{\mathrm{\Delta }_1}}\right)\left[2\omega _0^2\left(c_2R+{\displaystyle \frac{c_1}{\mathrm{\Delta }_1}}\right)\right]`$ (A18)
$`+k_BT\{{\displaystyle \frac{c_1}{\mathrm{\Delta }_1}}e^{\mathrm{\Delta }_1t}[22\omega _0^2c_2R{\displaystyle \frac{2\omega _0^2c_1}{\mathrm{\Delta }_1}}+e^{\mathrm{\Delta }_1t}(\mathrm{\Delta }_1c_1+{\displaystyle \frac{\omega _0^2c_1}{\mathrm{\Delta }_1}})]`$
$`{\displaystyle \frac{2c_2e^{\mathrm{\Delta }_2t}}{\mathrm{\Delta }_2^2+\beta ^2}}\left[1\omega _0^2c_2R+{\displaystyle \frac{\omega _0^2c_1}{\mathrm{\Delta }_1}}(e^{\mathrm{\Delta }_1t}1)\right]\left[\mathrm{\Delta }_2\mathrm{sin}(\beta t+\alpha )+\beta \mathrm{cos}(\beta t+\alpha )\right]`$
$`2c_1c_2e^{(\mathrm{\Delta }_1+\mathrm{\Delta }_2)t}\mathrm{sin}(\beta t+\alpha )`$
$`{\displaystyle \frac{\mathrm{\Delta }_2\beta \omega _0^2c_2^2e^{2\mathrm{\Delta }_2t}}{(\mathrm{\Delta }_2^2+\beta ^2)^2}}\mathrm{sin}2(\beta t+\alpha ){\displaystyle \frac{\beta ^2\omega _0^2c_2^2e^{2\mathrm{\Delta }_2t}}{(\mathrm{\Delta }_2^2+\beta ^2)^2}}`$
$`+[{\displaystyle \frac{\omega _0^2(2\beta ^2\mathrm{\Delta }_2^2)}{(\mathrm{\Delta }_2^2+\beta ^2)^2}}1]c_2^2e^{2\mathrm{\Delta }_2t}\mathrm{sin}^2(\beta t+\alpha )\}`$
with

$$R=\frac{1}{\mathrm{\Delta }_2^2+\beta ^2}(\mathrm{\Delta }_2\mathrm{sin}\alpha +\beta \mathrm{cos}\alpha )$$
(A19)
and

$`\sigma _{xx}^{2(e)}(t)`$ $`=`$ $`2D\left({\displaystyle \frac{\kappa _0}{g_0}}\tau _c\right)^2[c_1^2(\omega _0^4+\mathrm{\Delta }_1^4+2\omega _0^2\mathrm{\Delta }_1^2)I_A(t)`$ (A23)
$`+c_2^2\{\omega _0^4+(\mathrm{\Delta }_2^2\beta ^2)^24\beta ^2\mathrm{\Delta }_2^2+2\omega _0^2(\mathrm{\Delta }_2^2\beta ^2)\}I_B(t)`$
$`+2c_1c_2\{\omega _0^4+\mathrm{\Delta }_1^2(\mathrm{\Delta }_2^2\beta ^2)+\omega _0^2(\mathrm{\Delta }_1^2+\mathrm{\Delta }_2^2\beta ^2)\}I_C(t)`$
$`2c_2^2\beta \mathrm{\Delta }_2(\mathrm{\Delta }_2^2\beta ^2+\omega _0^2)I_D(t)+4c_2^2\beta ^2\mathrm{\Delta }_2^2I_E(t)4c_1c_2\beta \mathrm{\Delta }_2(\mathrm{\Delta }_1^2+\omega _0^2)I_F(t)].`$
Here the $`I`$’s are defined by

$`I_A(t)`$ $`=`$ $`{\displaystyle _0^t}e^{2\mathrm{\Delta }_1t}𝑑t,`$ (A25)
$`I_B(t)`$ $`=`$ $`{\displaystyle _0^t}e^{2\mathrm{\Delta }_2t}\mathrm{sin}^2(\beta t+\alpha )𝑑t,`$ (A26)
$`I_C(t)`$ $`=`$ $`{\displaystyle _0^t}e^{(\mathrm{\Delta }_1+\mathrm{\Delta }_2)t}\mathrm{sin}(\beta t+\alpha )𝑑t,`$ (A27)
$`I_D(t)`$ $`=`$ $`{\displaystyle _0^t}e^{2\mathrm{\Delta }_2t}\mathrm{sin}2(\beta t+\alpha )𝑑t,`$ (A28)
$`I_E(t)`$ $`=`$ $`{\displaystyle _0^t}e^{2\mathrm{\Delta }_2t}𝑑t\mathrm{and}`$ (A29)
$`I_F(t)`$ $`=`$ $`{\displaystyle _0^t}e^{(\mathrm{\Delta }_1+\mathrm{\Delta }_2)t}\mathrm{cos}(\beta t+\alpha )𝑑t.`$ (A30)
Similarly

$`\sigma _{vv}^2(t)=\sigma _{vv}^{2(i)}(t)+\sigma _{vv}^{2(e)}(t)`$
where

$`\sigma _{vv}^{2(i)}(t)`$ $`=`$ $`k_BT[(\mathrm{\Delta }_1^2+\omega _0^2)c_1^2e^{2\mathrm{\Delta }_1t}+\beta ^2c_2^2e^{2\mathrm{\Delta }_2t}`$ (A33)
$`\beta \mathrm{\Delta }_2c_2^2e^{2\mathrm{\Delta }_2t}\mathrm{sin}2(\beta t+\alpha )+(\mathrm{\Delta }_2^2+\omega _0^2\beta ^2)c_2^2e^{2\mathrm{\Delta }_2t}\mathrm{sin}^2(\beta t+\alpha )`$
$`+e^{(\mathrm{\Delta }_1+\mathrm{\Delta }_2)t}\{2c_1c_2(\omega _0^2+\mathrm{\Delta }_1\mathrm{\Delta }_2)\mathrm{sin}(\beta t+\alpha )2\mathrm{\Delta }_1\beta c_1c_2\mathrm{cos}(\beta t+\alpha )\}]`$
and

$`\sigma _{vv}^{2(e)}(t)`$ $`=`$ $`2D\left({\displaystyle \frac{\kappa _0}{g_0}}\tau _c\right)^2[c_1^2\mathrm{\Delta }_1^2(\omega _0^2+\mathrm{\Delta }_1^2)^2I_A(t)`$ (A39)
$`+c_2^2\{(\omega _0^2+\mathrm{\Delta }_2^23\beta ^2)^2\mathrm{\Delta }_2^2(\omega _0^2+3\mathrm{\Delta }_2^2\beta ^2)^2\beta ^2\}I_B(t)`$
$`+2c_1c_2\mathrm{\Delta }_1\mathrm{\Delta }_2(\omega _0^2+\mathrm{\Delta }_1^2)(\omega _0^23\beta ^2+\mathrm{\Delta }_2^2)I_C(t)`$
$`c_2^2\beta \mathrm{\Delta }_2(3\mathrm{\Delta }_2^2\beta ^2+\omega _0^2)(\omega _0^23\beta ^2+\mathrm{\Delta }_2^2)I_D(t)`$
$`+c_2^2\beta ^2(\omega _0^2\beta ^2+3\mathrm{\Delta }_2^2)I_E(t)`$
$`2c_1c_2\beta \mathrm{\Delta }_1(\mathrm{\Delta }_1^2+\omega _0^2)(\omega _0^2+3\mathrm{\Delta }_2^2\beta ^2)I_F(t)]`$
where, $`I`$’s are defined in Eq.(A25-A30). The explicit expression for $`\sigma _{xv}^2(t)`$ can be derived from Eq.(35). In the limit $`t\mathrm{}`$ we calculate the stationary values of the variances. The variances $`\sigma _{xx}^2(\mathrm{})`$, $`\sigma _{vv}^2(\mathrm{})`$ and $`\sigma _{xv}^2(\mathrm{})`$ yield $`\varphi (\mathrm{})`$ and $`\psi (\mathrm{})`$ and other relevant quantities.


text:
# Remark on Periodic Solutions of Non Linear Oscillators

## 1 Introduction

The non linear nonharmonic motion of an oscillator may be given by the following differential equation

$$u^{\prime \prime }+\omega ^2u=\beta u^2$$
(1)
$`\beta `$ and $`\omega `$ being constants, with initial conditions

$$u(0)=a_0,u^{}(0)=0$$
(2)
To solve this problem A. Shidfar and A. Sadeghi , have given two series solutions. They describe a general approach in which the differential equation, rather than the solutions series , is majorized.
Notice that if we write $`a_0=\frac{\omega ^2}{\beta }`$ then

$$u(t)\frac{\omega ^2}{\beta }$$
(3)
is a trivial solution of (1) and (2).
They gave a series solutions of (1) and (2), which includes (3) as a special case.

By writing

$$u=v\frac{\omega ^2}{2\beta }$$

the problem becomes

$$v^{\prime \prime }+\beta v^2=\frac{\omega ^4}{4\beta }$$
(4)
under the initial conditions

$`\{\begin{array}{cc}v(0)=a_0+\frac{\omega ^2}{2\beta }\hfill & \\ v^{}(0)=0\hfill & \end{array}`$ (5)
The method of  consists to solve equations (4) and (5) in the form

$$v(t)=c_0+c_1\mathrm{sin}\omega t+c_2\mathrm{sin}^2\omega t+c_3\mathrm{sin}^3\omega t+\mathrm{}..$$
(6)
where $`c_i,i=0,1,2,\mathrm{}.`$ are coefficients to be determined by the substitution of (6) in (4).
In fact, $`\omega =\frac{\pi }{T}`$ where $`T`$ is the period of the solution, which can be expressed in terms of the Weierstrass function $`\mathrm{}(z,2T,2T^{})`$.
So, we find that

$$2\omega ^2c_2+\beta c_0^2=\frac{\omega ^2}{4\beta }.$$

For $`n1`$, the recursion formula for these coefficients are

$$(n+1)(n+2)c_{n+2}=n^2c_n\frac{\beta }{\omega ^2}\underset{r=0}{\overset{n}{}}c_rc_{nr}.$$
(7)
Equations (5) and (6) imply that $`c_1=0`$. Relations (7) yields

$$c_3=0,c_5=0,\mathrm{}..$$

The even order coefficients simply are

$$c_0=a_0+\frac{\omega ^2}{2\beta },$$

$$c_2=\frac{a_0}{2\omega ^2}(\omega ^2+a_0\beta ),$$

$$c_4=\frac{\beta }{6\omega ^2}a_0(\omega ^2+a_0\beta )(\frac{3}{4}\frac{a_0\beta }{2\omega ^2}),$$

$$c_6=\frac{\beta }{180\omega ^2}a_0(\omega ^2+a_0\beta )(\frac{3}{4}\frac{a_0\beta }{2\omega ^2})(15\frac{2a_0\beta }{\omega ^2})\frac{\beta a_0^2}{120\omega ^6}(\omega ^2+a_0\beta )^2,$$

etc.
The coefficient $`c_0`$ follows from the condition (5). The solution for the equations (4) and (5) can now be written as

$$u(t)=a_0\frac{a_0}{2\omega ^2}(\omega ^2+a_0\beta )\mathrm{sin}^2\omega t+\mathrm{}.$$
(8)
Relations of the coefficients and further induction show that $`c_{2i},i=1,2,\mathrm{}`$ all vanish for $`a_0=\frac{\omega ^2}{\beta }.`$ So, the trivial solution (3) is included in (6) as a special case.

## 2 Convergence of the solutions

We now show the convergence of these series. In  one proved the following

Lemma 1 The serie (6) solution of Equation (4)-(5) is absolutely convergent for all $`t`$.

Proof We firstly note that if $`c_0>0,c_2>0`$ and $`\beta <0,`$ then all coefficients $`c_n`$ in the serie expansion (6) are positive. Indeed, we may write

$$\underset{n0}{}((n+1)(n+2)c_{n+2}=\underset{n0}{}n^2c_n\frac{\beta }{\omega ^2}(\underset{n0}{}c_n)^2+\frac{\beta }{\omega ^2}c_0^2,$$

or

$$\beta (\underset{n0}{}c_n)^2+\omega ^2\underset{n0}{}nc_n=\beta c_0^22\omega ^2c_2.$$
(9)
Since the right hand side of (9) is finite and $`c_i`$ are positive, the series $`_{n0}c_n`$ converges. Following , if we put

$$c_0^{}=c_0,c_1^{}=0,c_2^{}=c_2$$

and for $`n2`$

$$c_{n+2}^{}=\frac{n^2}{(n+1)(n+2)}c_n^{}+\frac{\beta }{\omega ^2(n+1)(n+2)}\underset{r=0}{\overset{n}{}}c_r^{}c_{nr}^{},$$

then the series $`_{n0}c_n^{}`$ converges. Since $`c_nc_n^{}`$, it follows that the solution series (6) is absolutely convergent, and hence the series expansion solution of (1)-(2) converges for all $`t.`$

We notice that we may deduce Lemma 1 a previous result concerning Equation (4).
We have shown that the coefficients verify a more general properties. Indeed, we have

Lemma 2 For any positive number $`ϵ`$ small enough (but $`ϵ0`$), there exists a positive constant $`k`$ verifying

$$k<\frac{\beta }{\omega ^2}\frac{3}{4}ϵ4^{ϵ\frac{1}{2}}$$

such that the coefficients $`c_n`$ of the series expansion (6) solution of the differential equation (4)-(5) satisfy the inequality

$$c_n<\frac{k}{n^{\frac{3}{2}ϵ}}.$$
(10)
Proof We first notice that Lemma 2 gives an optimal result, because our method do not run for $`ϵ=0`$.
The coefficients $`c_n`$ of the power series solution, satisfy the recursion formula (7). We shall prove there exist two positive constants $`k>0,`$ and $`\alpha >1`$, such that the following inequality holds

$$c_n<\frac{k}{n^\alpha }$$

for any integer $`n1`$. Suppose for any $`np`$ , we get $`c_n<\frac{k}{n^\alpha }`$. In particular, it implies that

$$\underset{0<r<p}{}c_rc_{pr}<\underset{0<r<p}{}\frac{k^2}{r^\alpha (rp)^\alpha }\frac{k^2}{(p1)^{\alpha 1}}.$$

Equality (7) gives

$$c_{p+2}=\frac{p^22\frac{\beta }{\omega ^2}c_0}{(p+1)(p+2)}c_p\frac{\beta }{\omega ^2(p+1)(p+2)}\underset{r=1}{\overset{p1}{}}c_rc_{pr}.$$

Thus, if we prove the following inequality

$$\frac{p^22\frac{\beta }{\omega ^2}c_0}{(p+1)(p+2)}\frac{k}{p^\alpha }+\frac{\beta }{\omega ^2(p+1)(p+2)}\frac{k^2}{(p1)^{\alpha 1}}\frac{k}{(p+2)^\alpha }$$
(11)
so

$$c_{p+2}<\frac{k}{(p+2)^\alpha }$$
(12)
Notice that (11) implies that

$$\frac{p^22\frac{\beta }{\omega ^2}c_0}{(p+1)(p+2)}\frac{k}{p^\alpha }\frac{k}{(p+2)^\alpha }.$$

Thus, a necessary condition to (11) holds is : $`\alpha \frac{3}{2}.`$ Inequality (11) is equivalent to

$$k<\frac{\beta }{\omega ^2}pf(p)g(p)$$
(13)
where

$$f(p)=\frac{p+1}{p}(\frac{p1}{p+2})^{\alpha 1}$$

$$g(p)=1\frac{(p^2\frac{\beta }{\omega ^2}c_0)(p+2)^{\alpha 1}}{(p+1)p^\alpha }$$

By using MAPLE, we are able to prove that $`f(p)`$ is an increasing positive function in $`p`$ . Moreover, for any $`p1`$, $`f(p)`$ is minorated

$$f(p)(\frac{3}{2})4^{1\alpha }.$$

The function $`g(p)`$ is such that

$$pg(p)=p\frac{(p^2\frac{\beta }{\omega ^2}c_0)(\frac{p+2}{p})^{\alpha 1}}{(p+1)}$$

is a strictly decreasing and bounded function .
More exactly, we may calculate the lower bound

$$g(p)>\frac{(32\alpha )}{p}.$$

Thus, if $`(32\alpha )=ϵ>0,`$ it suffices to choice

$$k(\frac{3}{2})4^{1\alpha }(32\alpha )$$

to inequality (13) holds.

Remark for the case $`ϵ=0`$ :
Notice that the choice of $`k`$ depends on $`\alpha `$ value.
For $`\alpha =\frac{3}{2},`$ we then prove by MAPLE that the function

$$pg(p)=p\frac{(p^2\frac{3}{2\omega ^2}c_0)(\frac{p+2}{p})^{\frac{1}{2}}}{(p+1)}$$

is positive and strictly decreasing to $`0`$. While $`p^2g(p)`$ is a bounded function.
Moreover, it appears that $`pf(p)g(p)`$ is a decreasing function which tends to $`0`$ when $`p`$ tends to infinity. Thus, our method falls since it do not permit to determine a non negative constant $`k`$.

Following , it is interesting to write the power series solution for the system (4)-(5),

$$v(x)=\underset{n=0}{\overset{\mathrm{}}{}}b_nx^n.$$
(14)
We find again that

$$b_{2p+1}=0p=0,1,2,\mathrm{},$$

while

$$b_0=a_0+\frac{\omega ^2}{2\beta },$$

$$2b_2+\beta b_0^2=\frac{\omega ^4}{4\beta },$$

$$(n+2)(n+1)b_{n+2}=\beta \underset{r=0}{\overset{r=n}{}}b_rb_{nr},$$

where $`n`$ is even and non zero.
The coefficients $`b_{2p},p=0,1,2,\mathrm{},`$ again vanish for $`a_0=\frac{\omega ^2}{\beta }.`$ We may verify that the solutions (6) and (14) are identical.
The latter method permits to compare approximate solutions of the anharmonic motion of the oscillator.

## 3 Another differential equation

We now examine the following differential equation

$$u^{\prime \prime }+\omega ^2u=\beta u^3$$
(15)
$`\beta `$ and $`\omega `$ being constants, with initial conditions

$$u(0)=a_0,u^{}(0)=0.$$
(16)
We put $`v=\frac{u}{a_0}`$ and $`t=\omega x.`$ We then obtain from (15) and (16)

$$\frac{d^2v}{dt^2}+v+\beta v^3=0,v(0)=1,v^{}(0)=0$$
(17)
where $`\beta =\frac{\beta a_0^2}{\omega }.`$
A. Shidfar and A. Sadeghi  solved (17) by series method in Sinus power

$$v(t)=c_0+c_1\mathrm{sin}\omega t+c_2\mathrm{sin}^2\omega t+c_3\mathrm{sin}^3\omega t+\mathrm{}..$$
(18)
Here, $`\omega =\frac{\pi }{T}`$ where $`T`$ is the period of the solution, which can be expressed in terms of the Jacobi function $`sn(z,2T,2T^{})`$.
So,

$$c_0=a_0.$$

For $`n1`$, we get the recursion formula

$$(n+1)(n+2)c_{n+2}=n^2c_n\frac{\beta }{\omega ^2}\underset{r=0}{\overset{n}{}}\underset{m=0}{\overset{nr}{}}c_mc_rc_{nmr}.$$
(19)
Under some conditions, they proved estimates of the coefficients

$$c_nR^n$$

where $`\frac{1}{R}`$ is a radius of convergence.
In fact, we may prove an analog of Lemma 2 for this equation. Indeed, we have

$$\underset{r=0}{\overset{n}{}}\underset{m=0}{\overset{nr}{}}c_mc_rc_{nmr}=2c_0^2c_n+2c_0c_1c_{n1}+c_0\underset{m=0}{\overset{n}{}}c_nc_{nm}+\underset{r=2}{\overset{n2}{}}\underset{m=0}{\overset{nr}{}}c_mc_rc_{nmr}.$$

Then,

$$\underset{r=2}{\overset{n2}{}}\underset{m=0}{\overset{nr}{}}c_mc_rc_{nmr}=\underset{r=2}{\overset{n2}{}}c_r[2c_0c_{nr}+\underset{m=1}{\overset{nr1}{}}c_mc_{nmr}]$$

$$<\underset{r=2}{\overset{n2}{}}c_r[2c_0c_{nr}+\frac{k^2}{(nr1)^{\alpha 1}}]$$

$$<2c_0k^2\underset{r=2}{\overset{n2}{}}\frac{1}{r^\alpha (nr)^\alpha }+k^3\underset{r=2}{\overset{n2}{}}\frac{1}{r^\alpha (nr1)^{\alpha 1}}<\frac{2c_0k^2}{(n1)^{\alpha 1}}+\frac{k^3}{(n2)^{\alpha 1}}$$

Finally,

$$c_{p+2}<\frac{k}{(p+2)^\alpha }$$

as soon as the non negative constant $`k`$ satisfies the inequality

$$\frac{k}{n^{\alpha 2}}+\frac{2c_0k^2}{n^\alpha }+\frac{2c_0c_1k}{(n1)^\alpha }+\frac{3c_0k^2}{(n2)^{\alpha 1}}+\frac{k^3}{(n2)^{\alpha 1}}<\frac{(n+1)k}{(n+2)^{\alpha 1}},$$

So,

$$\frac{1}{n^{\alpha 2}}+\frac{k^2+3c_0k}{(n2)^{\alpha 1}}+\frac{2c_0+2c_0c_1}{(n1)^\alpha }<\frac{(n+1)}{(n+2)^{\alpha 1}}.$$

By using MAPLE, we verify it is possible to find a such constant. Moreover, we find again the necessary condition : $`\alpha \frac{3}{2},`$ since we get $`\frac{1}{n^{\alpha 2}}<\frac{(n+1)}{(n+2)^{\alpha 1}}.`$

General remarks : It is wellknown from the theory of elliptic functions that solutions of equations (4) and (17) are related. This allows to express the series expansion of the solution of (17) from a series expansion of a solution of (4) and conversely. Indeed, one has

$$\mathrm{}(z)=C\frac{\delta ^2}{sn^2(\delta z)},$$

where $`\mathrm{}(z)`$ is the elliptic Weierstrass function and $`sn(u)`$ is the Jacobi fonction , $`\delta `$ is a constant , only dependent on the initial parameters.
Notice that series expansion of the $`sn(u)`$ in sinus power was given in a previous paper (see Proposition (2.1) in ).

REFERENCES

 A. Shidfar and A. Sadeghi Some Series Solutions of the Anharmonic Motion Equation, J. Math. Anal. Appl. 120, p. 488-493 (1986).

 A. Shidfar and A. Sadeghi The Periodic Solutions of Certain Non-linear Oscillators, Appl. Math. Lett., vol 3, n 4, p. 21-24, (1990).

 R. Chouikha Sur des developpements de fonctions elliptiques, Publ. Math. Fac des Sci de Besancon, Fasc. Th des Nombres, p.1-9, (1989).

 R. Chouikha Fonctions elliptiques et bifurcations d’equations differentielles, Canad. Math. Bull., vol. 40 (3), p. 276-284, (1997)