A new approach in an analytical method for diffusion dynamics for the presence of delocalized sink in a potential well: Application to different potential curves

TL;DR

This paper introduces a novel analytical approach to solve the one-dimensional Fokker-Planck equation with delocalized sinks, simplifying calculations of rate constants across various potential energy curves.

Contribution

The authors present a new method that avoids matrix equation solutions for rate constants, applicable to multiple potential energy profiles.

Findings

Method successfully applied to flat, linear, and parabolic potentials.

Simplifies calculation process compared to previous matrix-based methods.

Provides analytical solutions in the Laplace domain for diffusion with delocalized sinks.

Abstract

We provide a new approach to solve one dimension Fokker-Planck equation in the Laplace domain for the case where a particle is evolving in a potential energy curve in the presence of general delocalized sink. We also calculate rate constants in the presence of non-localized sink on different potential energy curves. In the previous method, we need to solve matrix equations to calculate rate constant but in our method, it is not required. We also calculate the rate constant by using the method to some known potential curve, viz, flat potential, linear potential, and parabolic potential.

Tables

Table 1. Table 1 : A comparison of rate constant with number of sink function n under the same value of rest of parameters. K I ​ S subscript 𝐾 𝐼 𝑆 K_{IS} is from Samanta-Ghosh’s method and K I ​ O subscript 𝐾 𝐼 𝑂 K_{IO} is from our method.

n	$K_{IS}$	$K_{IO}$
1	1.13977	1.13977
2	1.17566	1.17566
3	1.18468	1.18462
4	1.18693	1.18682
5	1.18749	1.18737
6	1.18763	1.18750

Equations

$$\frac{\partial P(x,t)}{\partial t}=[\mathcal{L}-S(x)-k_r]P(x,t).$$

$$\mathcal{L}=D(\frac{{\partial }^2}{\partial x^2}+\frac{1}{k_{B}T}\frac{\partial }{\partial x}\frac{dU(x)}{dx}).$$

$$[s+k_r-\mathcal{L}+S(x)]\mathcal{P}(x,s)=P(x,0),$$

$$\mathcal{G}(x,s+k_r|x_0)=<x|[s+k_r-\mathcal{L}+S(x){]}^{-1}|x_0>.$$

$$\mathcal{P}(x,s+k_r|x_0)={\int }_{-\infty }^{\infty }d{x}^{\prime }\mathcal{G}(x,s+k_r|x^{\prime })P(x^{\prime },0).$$

$$[s+k_r-\mathcal{L}+S(x){]}^{-1}=[s+k_r-\mathcal{L}{]}^{-1}-[s+k_r-\mathcal{L}{]}^{-1}S(x)[s+k_r-\mathcal{L}+S(x){]}^{-1}.$$

$${\mathcal{G}}_1(x,s+k_r|x_0)={\mathcal{G}}_0(x,s+k_r|x_0)-\frac{k_1{\mathcal{G}}_0(x,s+k_r|x_1){\mathcal{G}}_0(x_1,s+k_r|x_0)}{1+k_1{\mathcal{G}}_0(x_1,s+k_r|x_1)}.$$

$${\mathcal{G}}_0(x,s+k_r|x_0)=<x|[s+k_r-\mathcal{L}{]}^{-1}|x_0>.$$

$$[s+k_r+k_0[S_1(x)+S_2(x)]-\mathcal{L}{]}^{-1}=[s+k_r+k_0{S}_1(x)-\mathcal{L}{]}^{-1}$$

$$-[s+k_r+k_0{S}_1(x)-\mathcal{L}{]}^{-1}{k}_0{S}_2(x)[s+k_r+k_0[S_1(x)+S_2(x)]-\mathcal{L}{]}^{-1}.$$

$${\mathcal{G}}_2(x,s+k_r|x_0)={\mathcal{G}}_1(x,s+k_r|x_0)-\frac{k_0{k}_2{\mathcal{G}}_1(x,s+k_r|x_2){\mathcal{G}}_1(x_2,s+k_r|x_0)}{1+k_0{k}_2{\mathcal{G}}_1(x_2,s+k_r|x_2)}.$$

$${\mathcal{G}}_n(x,s+k_r|x_0)={\mathcal{G}}_{n-1}(x,s+k_r|x_0)-\frac{k_0{k}_i{\mathcal{G}}_{n-1}(x,s+k_r|x_n){\mathcal{G}}_{n-1}(x_n,s+k_r|x_0)}{1+k_0{k}_n{\mathcal{G}}_{n-1}(x_n,s+k_r|x_n)}.$$

$$K_I^{-1}={\mathcal{P}}_e(0)$$

$${\mathcal{P}}_0(s)={\mathcal{G}}_0(s+k_r|x_0)={\int }_{-\infty }^{\infty }{\mathcal{G}}_0(x,s+k_r|x_0)dx=\frac{1}{s+k_r}.$$

$${\mathcal{P}}_1(s)={\mathcal{G}}_1(s+k_r|x_0)={\mathcal{G}}_0(s+k_r|x_0)-\frac{k_1{\mathcal{G}}_0(s+k_r|x_1){\mathcal{G}}_0(x_1,s+k_r|x_0)}{1+k_1{\mathcal{G}}_0(x_1,s|x_1)}$$

$$K_{I1}^{-1}=\frac{1}{k_r}[1-\frac{k_1{\mathcal{G}}_0(x_1,k_r|x_0)}{1+k_1{\mathcal{G}}_0(x_1,k_r|x_1)}].$$

$$K_{I2}^{-1}=K_{I1}^{-1}-\frac{k_2{\mathcal{G}}_1(x_2,s+k_r|x_0)}{1+k_2{\mathcal{G}}_1(x_2,s+k_r|x_2)}{|}_{s=0}{K}_{I1}^{-1}{|}_{x_0=x_2}.$$

$$K_{In}^{-1}=K_{In-1}^{-1}-\frac{k_n{\mathcal{G}}_{n-1}(x_n,s+k_r|x_0)}{1+k_n{\mathcal{G}}_{n-1}(x_n,s+k_r|x_n)}{|}_{s=0}{K}_{In-1}^{-1}{|}_{x_0=x_n}.$$

$${\mathcal{G}}_0(x,s+k_r|x_0)=\frac{e^{-\sqrt{\frac{s+k_r}{D}}|x-x_0|}}{2\sqrt{D(s+k_r)}}.$$

$${\mathcal{P}}_1(s)=\frac{1}{s+k_r}[1-\frac{k_1{e}^{-\sqrt{\frac{s+k_r}{D}}(x_1-x_0)}}{2\sqrt{D(s+k_r)}+k_1}].$$

$${K_{I1}}^{-1}=\frac{1}{k_r}[1-\frac{k_1{e}^{-\sqrt{\frac{k_r}{D}}(x_1-x_0)}}{2\sqrt{D{k}_r}+k_1}].$$

$${K_{In}}^{-1}={K_{In-1}}^{-1}-\frac{k_0{k}_i{\mathcal{G}}_{n-1}(x_n,s+k_r|x_0)}{1+k_0{k}_n{\mathcal{G}}_{n-1}(x_n,s+k_r|x_n)}{|}_{s=0}{K}_{In-1}^{-1}{|}_{x_0=x_n}.$$

$$[s+k_r-\mathcal{L}+S(x)]\mathcal{P}(x,s)=\delta (x-x_0).$$

$${\mathcal{G}}_0(x,s+k_r|x_0)=\frac{e^{-(|x|-|x_0|)/2l}}{\Gamma \sqrt{1+\frac{4ls+k_r}{\Gamma }}}[e^{-\sqrt{1+\frac{4ls+k_r}{\Gamma }}|x-x_0|/2l}+\frac{e^{-\sqrt{1+\frac{4ls+k_r}{\Gamma }}(|x|+|x_0|)/2l}}{\sqrt{1+\frac{4ls+k_r}{\Gamma }}-1}].$$

$$K_{I1}^{-1}=\underset{s\to 0}{\lim }\frac{1}{s}[1-\frac{k_1{\mathcal{G}}_0(x_1,s|x_0)}{1+k_1{\mathcal{G}}_0(x_1,s|x_1)}].$$

$$K_{I1}^{-1}=\frac{1}{k_1\Gamma }(2l{k}_1{e}^{\frac{|x_1|}{l}}-2l{k}_1{e}^{\frac{|x_0|+|x_1|-|x_0-x_1|}{2l}}+2l\Gamma e^{\frac{|x_1|}{l}}+k_1|x_0|-k_1|x_1|).$$

$$\begin{array}{rl}\underset{s\to 0}{\lim }{\mathcal{G}}_1(x,s+k_r|x_j)& =\frac{1}{\Gamma }{e}^{-}\frac{2|x|+|x-x_j|+|x-x_1|+|x_j-x_1|}{2l}\\ & \times [e^{-}\frac{|x|+|x_j|+|x-x_1|+|x_j-x_1|}{2l}\\ & -e^{-}\frac{|x-x_j|+|x_j|+|x-x_1|+|x_1|}{2l}\\ & -e^{-}\frac{|x|+|x-x_j|+|x_j-x_1|+|x_1|}{2l}\\ & +(1+\frac{\Gamma }{k_1})e^{-}\frac{|x-x_j|+|x-x_1|+|x_j-x_1|+2|x_1|}{2l}].\end{array}$$

$$[s+k_r-D\frac{{\partial }^2}{\partial x^2}-\frac{D}{K_{B}T}\frac{\partial }{\partial x}(\frac{dU(x)}{dx})]{\mathcal{G}}_0=\delta (x-x_0)$$

$${\mathcal{G}}_0(x,s+k_r|x_0)=F(z,s+k_r|z_0)/(s+k_r),$$

$$F(z,s+k_r|z_0)=D_{\nu }(z)D_{\nu }(-z_0)exp[({z_0}^2-z^2)/4]\Gamma (1-\nu )[B/(2\pi D){]}^{1/2}.$$

Full text

††journal: Chemical Physics Letters

1 Introduction

The dynamics of diffusion motion of a particle in presence of sink on the potential energy surface (PES) provide a significant contribution to understanding various phenomena in chemical and biological process. It is modeled mathematically with the use of a Smoluchowski-like equation. This type of model has been used to calculate rates of diffusion-controlled reactions as well as cyclization of polymers chains in solutions[1, 2]. To describe the electron transfer reactions in polar solvent many authors have used the model of such a type[3, 4, 5]. A model of such a kind has been utilized by Sumi to explain the pressure influence on some isomerization reactions[6]. Even the one dimensional model of diffusion motion have been using to model several problems involving proton, an activated barrier crossing or barrierless reaction process such as relaxation from electronic excited state in solution [7, 8, 9, 10, 11]. In a biological process, the theoretical study of the phenomena including artificial photosynthetic machines, enzyme-ligand binding and home range formation in animal behavior [12, 13, 14] has often engaged the same diffusion equation. Due to absorption of light particle drifts over the excited state of PES and getting relaxed via radiative or nonradiative decay from anywhere from the surface. In most of the model the sink term often consider to be a single Dirac delta function on the excited PES. In 1984 Szabo, Lamm and Weiss [15] provided a generalized method when the sink is a sum of $\delta$ functions and the solution can be written in terms Green’s function when there is no sink. Later Samanta and Ghosh provided an extension to the previous method and explored it for a general delocalized sink [16]. But In their technique, it’s essentially needed to solve a matrix equation and moreover, the dimension of the matrix is proportional to a number of $\delta$ function sink present on the PES. So whenever we generalize the form of sink term a large number of $\delta$ functions is needed and thus the complexity of the rate calculation is therefore reflected.

The purpose of this paper is to construct a new method to find the solution of Smoluchowski like equation for general delocalized sink in the Laplace domain. Our method doesn’t include any matrix calculation. It is quite simple and easy to calculate the rate constant. We also know that it is not easy to calculate the non-radiative rate constant in the Laplace domain and special treatment has to make for different PES in calculating the same. The beauty of our method is if we have the non-radiative decay rate expression for a single $\delta$ function sink then we can calculate the same decay rate for the rest of $\delta$ functions sink smoothly by the recursion formula, which will be shown below. Further we apply this method for different potential energy curve (PEC) and have calculated rate constant numerically which is shown graphically.

2 Methodology

It is quite common to assume the motion on the PES to be one-dimensional and diffusive, the relevant coordinate being denoted by x. It’s additionally usual to assume that the motion on the potential energy curve (PEC) is overdamped. The probability distribution $P(x,t)$ corresponding to diffusion in a potential well $U(x)$ in the presence of sink function $S(x)$ is governed by the Fokker-Planck equation in the Smoluchowski form which is

[TABLE]

In the above,

[TABLE]

$k_{r}$ is position independent radiative decay rate and D is diffusion coefficient at temperature $T$ and $k_{B}$ is the Boltzmann constant. The above equation will often be called as a Smoluchawski equation when $U(x)$ is considered to be quadratic in $x$. Our method starts with taking the Laplace transform of Eq. (1), which is

[TABLE]

where $s$ is Laplace variable and $\mathscr{P}(x,s)$ denotes Laplace transform of $P(x,t)$. The initial distribution is assumed to be localized at $x_{0}$ and thus $P(x,0)=\delta(x-x_{0})$. The solution of the Eq. (2) in terms of Green’s function is given by

[TABLE]

The solution in Laplace domain will be

[TABLE]

By using the operator identity of quantum mechanics we can get

[TABLE]

The forms of the potential and also the sink function rely upon the physical problem of interest. In several cases, the potential is chosen to be parabolic and also the simplest sink function is the localized sink at an acceptable location $x_{i}$ and it is given by $S(x)=k_{i}\delta(x-x_{i})$. Although in majority problem of interest the consideration of non localized sink function is required for a proper description of dynamics. We consider the sink function as a linear combination of $\delta$ function. This form of the sink function was used by Szabo, Lamm, and Weiss [15]. Samanta and Ghosh [16] also used the same form of the sink function. This function is quite general because an arbitrary sink function can be written as $\int_{-\infty}^{\infty}dx^{\prime}S(x^{\prime})\delta(x-x^{\prime})$. Farther we can discretize the integral as $S=\sum_{i}S_{i}$, where $S_{i}=k_{i}\delta(x-x_{i})$ and $k_{i}$ is strength of the $i$th $\delta$ function sink. By putting the operator identity in Eq. (5) into Eq. (3) we get the Green’s function for a single sink function, $S_{1}=k_{1}\delta(x-x_{1})$, which is

[TABLE]

In the above $\mathscr{G}_{0}(x,s+k_{r}|x_{0})$, the Green’s function corresponding to propagation of the particle placed initially at $x_{0}$ when there is no sink on the potential energy curve, is defined as

[TABLE]

Now we add another sink function, $S_{2}=k_{2}\delta(x-x_{2})$, along with the first one and the corresponding operator identity in Eq. (5) will be

[TABLE]

[TABLE]

Like the case of a single sink function, we obtain the Green’s function for the presence of two sink functions and it is given by

[TABLE]

Obviously $\mathscr{G}_{1}(x,s+k_{r}|x_{0})$ will be calculated from Eq. (6). So the above derivations suggest that by knowing $\mathscr{G}_{0}(x,s+k_{r}|x_{0})$ one can find the solution for any number of $\delta$ function sink presence on the potential energy curve analytically. The general form of the Green’s function solution for $n$ number of $\delta$ function sink will be

[TABLE]

Here we see that if we are able to find $\mathscr{G}_{0}$ for a potential energy curve then we can find the solution when there are $n$ numbers of $\delta$ sink function on PEC analytically.

3 Calculation of Rate Constant

In the following we provide a general procedure for calculating the average rate constant when multiple sink is presented on PEC. The average rate constant $K_{I}$ in Laplace domain has form ${\mathscr{P}_{e}}^{-1}(s=0)$ and it can be found from $\mathscr{P}_{e}(x,s)$ through the expression $\mathscr{P}_{e}(s)=\int_{-\infty}^{\infty}\mathscr{P}(x,s)dx$. Thus

[TABLE]

For a potential energy curve $U(x)$ with initial distribution $P(x,0)=\delta(x-x_{0})$ we may find

[TABLE]

Which correspondingly refers that the total probability in the time domain is unity when there is no sink present on PEC. The Eq. (11) is also true for any value of $x_{0}$. Therefor for one $\delta$ function sink at $x_{1}$ we get

[TABLE]

Hence the inverse of average rate constant is

[TABLE]

$K_{I1}$ depends upon initial distribution $P(x,0)$. By using Eq. (9) we can calculate

[TABLE]

Here $K_{I1}^{-1}|_{x_{0}=x_{2}}=\int_{-\infty}^{\infty}\mathscr{G}_{1}(x,k_{r}|x_{2})dx$. So in the above equation we are able to express $K_{I2}$ in terms of $K_{I1}$ and $\mathscr{G}_{1}(x,s+k_{r}|x_{0})$. The general expression for the rate constant in case of $n$ number of $\delta$ sink function can be written as

[TABLE]

However for nonzero of $k_{r}$ there is no problem as such in the evaluation of numerical value for the decay rate. But in the absence of radiative decay, i.e., $k_{r}=0$, divergences come in the evaluation of ${K_{I1}}^{-1}$ in Eq. (13) in the limit $s\rightarrow 0$ if the stationary distribution is nonzero. So in order to calculate the non-radiative decay rate, a special technique is required to overcome the divergence. For that we have to find $K_{I1}^{-1}$ in the limit $k_{r}\rightarrow 0$.

3.1 Flat Potential

As a first example, we consider a simple case where a particle is influenced by a constant potential, for which the Green’s function

[TABLE]

Therefor the explicit form of $\mathscr{P}_{1}(s)$ is

[TABLE]

Hence average rate constant $K_{I1}$ for single sink is

[TABLE]

By using the above equation we can calculate $K_{I}^{-1}$ for more than one sink on the PEC. Numerical value of average rate constant have been calculated for different strength of sink function by using the below expression

[TABLE]

In Figure 1, we plot $K_{I}$ against a number of the sink, $n$, on PEC for three values of $k_{i}$, viz. 1, 2 and 3.

This figure is showing $K_{I}$ gets saturated for the presence of a higher number of sink function. Also overall $K_{I}$ is increasing with increase of strength of each $\delta$ function sink. In Table 1, a small comparison between our method and Samanta-Ghosh’s method have been made by calculating numerical values of $K_{I}$.

3.2 Linear Potential

We can get the Green’s function for the linear potential $U(x)=|x|$ from the following equation

[TABLE]

with $\mathscr{L}=D\frac{d^{2}}{dx^{2}}+\frac{D}{k_{B}T}\frac{d}{dx}\frac{|x|}{x}$. By using proper boundary condition, the Green’s function for the linear or V shaped potential will be [17]

[TABLE]

Where $\Gamma=D/k_{B}T$ and $l=D/\Gamma$. The solution for generalized sink can be determined by using Eq. (10). For non radiative decay rate calculation the expression for, in case of single $\delta$ function sink, average rate constant from Eq. (13)

[TABLE]

By using the Eq. (21), with $k_{r}=0$, and after calculating the limiting value of right hand site of of the above equation we find

[TABLE]

We also calculate $\mathscr{G}_{1}(x,s+k_{r}|x_{j})$ in the limit $s\rightarrow 0$, which will be needed to calculate rate constant for a higher number of sink. The relevant expression is

[TABLE]

The above expression is derived by employing the Green’s function of the linear potential. So the form of inverse of nonradiative decay rate for two $\delta$ function sink at positions $x_{1}$ and $x_{2}$ is nothing but the Eq. (14) with $k_{r}\rightarrow 0$. In that equation we use Eq. (24) to calculate the numerical value of $K_{I2}$. So we can easily calculate the rate constant for a higher number of $\delta$ function sink by using our recursion formula in Eq. (15) and the Eq. (24). In Figure (2), we have plotted $K_{In}$ against a number of sinks $n$ on PEC for the cases, one for when initial distribution is at the left side of all sinks and another for when it is at the right side of all sinks.

3.3 Parabolic Potential

In the following, we get the propagator for the parabolic potential $U(x)=\frac{1}{2}\mu{\omega}^{2}x^{2}$, $\mu$ is mass and $\omega$ is frequency, by solving the following equation:

[TABLE]

to get

[TABLE]

with

[TABLE]

In the above we have introduced new variable $z$ such that $z=x\sqrt{\frac{B}{D}}$ and $z_{i}=x_{i}\sqrt{\frac{B}{D}}$ [9]. We also have consider $B=\mu\omega^{2}D/(k_{B}T)$ and $\nu=-s/B$. $D_{\nu}$ are parabolic cylinder functions and $\Gamma(\nu)$ is the gamma function.

Now we consider Eq. (13) for this PEC, which is

[TABLE]

In the limit $\nu\rightarrow 0$ we know $D_{\nu}(z)=0$. So in this limit of $\nu$, as $k_{r}\rightarrow 0$, we find $F(z_{1},k_{r}|z_{0})$ and $F(z_{1},k_{r}|z_{1})\rightarrow exp(-{z_{1}}^{2}/2)[B/(2\pi D)^{1/2}]$ so that $k_{1}F(z_{1},k_{r}|z_{0})/k_{r}+k_{1}F(z_{1},k_{r}|z_{1})$ $\rightarrow 1$. Therefor Eq. (28) can be written as

[TABLE]

We consider that $x_{0}<x_{1}$, without loss of generality, so that the initial distribution is placed to the left of the sink and the above equation can farther be evaluated as

[TABLE]

After using proper expression for $D_{\nu}(-z)$ we can obtain

[TABLE]

Where $erf(z)$ is a error function and has a form, $erf(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}dt$. We can also find

[TABLE]

Using the above general expression and the Eq. (31) we may calculate $K_{I2}$ and the rate constant corresponding to a higher number of $\delta$ function sink when $k_{r}=0$. We have calculated nonradiative decay rate numerically, graphically shown in Figure (3), in presence of more than one sink on the PEC when $P(x,0)$ is sharply picked at the center of the PEC and sinks are considered on the right side of $P(x,0)$.

4 Summary and Conclusions

Our Eq. (10) is quit general for generalized sinks present on PEC in Laplace domain. Since the form of sink is most likely to be a Gaussian function and we can approximate this form by choosing strength of each $\delta$ function sink properly. Eq. (15) represents the rate constant for a generalized sink. If we are able to find the nonradiative decay rate for a localized sink, then It is very easy to calculate a nonradiative decay rate for a generalized sink by using our method than the Samanta and Ghosh’s method.

5 Acknowledgments

The author (C.S.) would like to thank Indian Institute of Technology Mandi for Half-Time Research Assistantship fellowship. The author is grateful to Aniruddha Chakraborty for suggestion and valuable discussion.

References