Communication: Propagator for diffusive dynamics of an interacting molecular pair

Lee, Sangyoub; Son, Chang Yun; Sung, Jaeyoung; Chong, Song-Ho

doi:10.1063/1.3565476

We introduce a new method of solution for the Fredholm integral equations of the second kind. The method would be useful when the direct iterative approach leads to a divergent perturbation series solution. By using the method, we obtain an accurate expression of the propagator for diffusive dynamics of a pair of particles interacting via an arbitrary central potential and hydrodynamic interaction. We test the accuracy of the propagator expression by calculating the diffusion-controlled geminate and bimolecular reaction rates. It is shown that our propagator expression provides very accurate results for the whole time region.

The propagator G(r, t|r₀), which represents a probability density of finding a pair of particles at a separation r at time t, given that their initial separation was r₀, is one of the most important dynamic quantities in the liquid state physics.¹ Its knowledge enables the calculation of various dynamic quantities such as time correlation functions of the physical variables that depend on the particle separation, the first encounter time distribution of a pair of particles, the time-dependent survival probability of reactive molecules, and so on.^2,3

When the medium surrounding the particles is very viscous, the time evolution of G(r, t|r₀) is often described by the Smoluchowski equation²:

\begin{equation}\frac{{\partial G(r,t| {r_0 } )}}{{\partial t}} = D\frac{1}{{r^2 }}\frac{\partial }{{\partial r}}r^2 h(r)e^{ - U(r)} \frac{\partial }{{\partial r}}e^{U(r)} G(r,t| {r_0 } ).\end{equation}\vspace*{-5pt}

\frac{\partial G (r, t | r_{0})}{\partial t} = D \frac{1}{r^{2}} \frac{\partial}{\partial r} r^{2} h (r) e^{- U (r)} \frac{\partial}{\partial r} e^{U (r)} G (r, t | r_{0}) .

(1)

Here, U(r) denotes an arbitrary central interaction potential in units of k_BT, with k_B and T denoting the Boltzmann constant and the absolute temperature. The relative diffusion coefficient at large separation is D, but its value decreases with decreasing r due to the hydrodynamic interaction, which is described by the function h(r). When the hydrodynamic interaction can be described by the Oseen tensor and the hydrodynamic radii of the particles is given by σ/2 with σ denoting the contact distance, we have h(r) = 1 − 3σ/4r.² However, other forms of h(r) may also be used. The initial and boundary conditions associated with Eq. (1) are |$G(r,t = 0| {r_0 } ) = \delta (r - r_0)/(4\pi r_0^2)$| $G (r, t = 0 | r_{0}) = δ (r - r_{0}) / (4 π r_{0}^{2})$ ⁠, |$\mathop {\lim }\limits_{r \to \infty } G(r,t | {r_0 } ) = 0$| $lim_{r \to \infty} G (r, t | r_{0}) = 0$ ⁠, and |$ {{{\partial [e^{U(r)} G(r,t| {r_0 } )]} \mathord{/ {\vphantom {{\partial [e^{U(r)} G(r,t| {r_0 } )]} {\partial r}}} \kern-\nulldelimiterspace} {\partial r}}} |_{r = \sigma } = 0$| $\partial [e^{U (r)} G (r, t | r_{0})] / \partial r |_{r = σ} = 0$ ⁠.

An exact general expression for G(r, t|r₀) is still absent. When the diffusion coefficient is constant [h(r) = 1] and U(r) is given by the Coulomb potential, the Laplace transform expression |$\hat G(r,s| {r_0 } )$| $\hat{G} (r, s | r_{0})$ is known, but its usability is limited because it is given as an infinite series of complicated functions.⁴

In the present work, we will present a simple and accurate expression for the propagator. The derivation is based on a new method of solution for the linear integral equations and the transformation of variables given by

\begin{equation}\tilde r(r) = \left\{ {\int_r^\infty {dr_1 } \,\big[r_1 ^2 h(r_1)e^{ - U(r_1)}\big]^{ - 1} } \right\}^{ - 1},\end{equation}\vspace*{-5pt}

\tilde{r} (r) = {\{\int_{r}^{\infty} d r_{1} {[r_{1}^{2}, h, (r_{1}), e^{- U (r_{1})}]}^{- 1}\}}^{- 1},

(2)

\begin{equation}\tilde G(\tilde r,t\left| {\tilde r_0 } \right.) = e^{U(r)} \varphi (r_0)^{ - 1} G(r,t\left| {r_0 } \right.).\end{equation}

\tilde{G} (\tilde{r}, t |{\tilde{r}}_{0})) = e^{U (r)} ϕ {(r_{0})}^{- 1} G (r, t |r_{0})) .

(3)

The transformation in Eq. (2) is the generalization of the Flannery's transformation⁵ to the case with h(r) ≠ 1. In Eq. (3), φ(r) is a function defined by

\begin{equation}\varphi (r) = h(r)\left( {{{d\tilde r}}/{{dr}}} \right)^2 = {{({{{\tilde r^4 }}/{{r^4}}})e^{2U(r)} }}/{{h(r)}}.\end{equation}

ϕ (r) = h (r) {(d \tilde{r} / d r)}^{2} = ({\tilde{r}}^{4} / r^{4}) e^{2 U (r)} / h (r) .

(4)

We have constructed |$\tilde G(\tilde r,t| {\tilde r_0 } )$| $\tilde{G} (\tilde{r}, t | {\tilde{r}}_{0})$ such that its initial value is given by |$\tilde G(\tilde r,t = 0 | {\tilde r_0 } ) = \delta (\tilde r - \tilde r_0)/(4\pi \tilde r_0^2)$| $\tilde{G} (\tilde{r}, t = 0 | {\tilde{r}}_{0}) = δ (\tilde{r} - {\tilde{r}}_{0}) / (4 π {\tilde{r}}_{0}^{2})$ ⁠, where |$\tilde r_0\break = \tilde r(r_0)$| ${\tilde{r}}_{0} = \tilde{r} (r_{0})$ ⁠. From Eqs. (1) – (3), we obtain

\begin{eqnarray}\frac{{\partial \tilde G(\tilde r,t\left| {\tilde r_0 } \right.)}}{{\partial t}} = L(\tilde r)\tilde G(\tilde r,t\left| {\tilde r_0 } \right.);\quad L(\tilde r) = D\varphi (\tilde r)\frac{1}{{\tilde r^2 }}\frac{\partial }{{\partial \tilde r}}\tilde r^2 \frac{\partial }{{\partial \tilde r}}.\hskip-5pt\nonumber\\\end{eqnarray}\vspace*{-15pt}

\begin{matrix} \frac{\partial \tilde{G} (\tilde{r}, t |{\tilde{r}}_{0}))}{\partial t} = L (\tilde{r}) \tilde{G} (\tilde{r}, t |{\tilde{r}}_{0})); L (\tilde{r}) = D ϕ (\tilde{r}) \frac{1}{{\tilde{r}}^{2}} \frac{\partial}{\partial \tilde{r}} {\tilde{r}}^{2} \frac{\partial}{\partial \tilde{r}} . \end{matrix}

(5)

The boundary conditions are given by |$\mathop {\lim }_{\tilde r \to \infty } \tilde G(\tilde r,t | {\tilde r_0 } )\break = 0$| $\lim_{\tilde{r} \to \infty} \tilde{G} (\tilde{r}, t | {\tilde{r}}_{0}) = 0$ and |$ {{{\partial \tilde G(\tilde r,t | {\tilde r_0 } )} \mathord{/ {\vphantom {{\partial \tilde G(\tilde r,t | {\tilde r_0 } )} {\partial \tilde r}}} \kern-\nulldelimiterspace} {\partial \tilde r}}} |_{\tilde r = \tilde \sigma } = 0$| $\partial \tilde{G} (\tilde{r}, t | {\tilde{r}}_{0}) / \partial \tilde{r} |_{\tilde{r} = \tilde{σ}} = 0$ ⁠, where |$\tilde \sigma = \tilde r(\sigma)$| $\tilde{σ} = \tilde{r} (σ)$ ⁠. Then the Laplace transformation of Eq. (5) gives

\begin{eqnarray}s{\hat {\tilde G}}(\tilde r,s|\tilde r_0) - {{\delta (\tilde r - \tilde r_0)}}/{{4\pi \tilde r_0^2 }} = L(\tilde r){\hat {\tilde G}}(\tilde r,s|\tilde r_0).\end{eqnarray}\vspace*{-5pt}

\begin{matrix} s \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) - δ (\tilde{r} - {\tilde{r}}_{0}) / 4 π {\tilde{r}}_{0}^{2} = L (\tilde{r}) \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) . \end{matrix}

(6)

If φ(r) = 1, Eq. (6) reduces to that for a seemingly interaction-free problem with the solution given by²

\begin{eqnarray}{\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0) = \frac{1}{{8\pi D\zeta \tilde r\tilde r_0 }}\left[ {e^{ - \zeta |\tilde r - \tilde r_0 |} + \frac{{\tilde \sigma \zeta - 1}}{{\tilde \sigma \zeta + 1}}e^{ - \zeta (\tilde r + \tilde r_0 - 2\tilde \sigma)} } \right],\hskip-6pt\nonumber\\\end{eqnarray}

\begin{matrix} {\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{0}) = \frac{1}{8 π D ζ \tilde{r} {\tilde{r}}_{0}} [e^{- ζ | \tilde{r} - {\tilde{r}}_{0} |} + \frac{\tilde{σ} ζ - 1}{\tilde{σ} ζ + 1} e^{- ζ (\tilde{r} + {\tilde{r}}_{0} - 2 \tilde{σ})}], \end{matrix}

(7)

where ζ = (s/D)^1/2. For the general case with φ(r) ≠ 1, we can express |${\hat {\tilde G}}(\tilde r,s|\tilde r_0)$| $\hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0})$ in terms of |${\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0)$| ${\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{0})$ as

\begin{eqnarray}{\hat {\tilde G}}(\tilde r,s|\tilde r_0) &=& {\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0)\nonumber\\&& - \int {d{\bf \tilde r}_1 } \;{\hat {\tilde G}}(\tilde r,s|\tilde r_1)d(\tilde r_1)L_0 (\tilde r_1){\hat {\tilde G}}_0 (\tilde r_1,s|\tilde r_0),\nonumber\\[-4pt]\end{eqnarray}

\begin{matrix} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) & = & {\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{0}) \\ - \int d {\tilde{r}}_{1} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{1}) d ({\tilde{r}}_{1}) L_{0} ({\tilde{r}}_{1}) {\hat{\tilde{G}}}_{0} ({\tilde{r}}_{1}, s | {\tilde{r}}_{0}), \end{matrix}

(8)

where |$\int \!\!{d{\bf \tilde r}} = \int_{\tilde \sigma }^\infty {d\tilde r} \,4\pi \tilde r^2$| $\int d \tilde{r} = \int_{\tilde{σ}}^{\infty} d \tilde{r} 4 π {\tilde{r}}^{2}$ ⁠, |$d(\tilde r) = 1 - \varphi (\tilde r)$| $d (\tilde{r}) = 1 - ϕ (\tilde{r})$ ⁠, and |$L_0 (\tilde r)\break = D(1/\tilde r^2) \,(\partial /\partial \tilde r)\tilde r^2 (\partial /\partial \tilde r)$| $L_{0} (\tilde{r}) = D (1 / {\tilde{r}}^{2}) (\partial / \partial \tilde{r}) {\tilde{r}}^{2} (\partial / \partial \tilde{r})$ ⁠. For notational convenience, we will not distinguish the dimensionless function of |$\tilde r$| $\tilde{r}$ from that of r representing the same quantity; that is, |$\varphi (\tilde r(r)) = \varphi (r)$| $ϕ (\tilde{r} (r)) = ϕ (r)$ and so on. Noting that |${\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0)$| ${\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{0})$ satisfies Eq. (6) with |$L(\tilde r)$| $L (\tilde{r})$ replaced by |$L_0 (\tilde r)$| $L_{0} (\tilde{r})$ ⁠, from Eq. (8) we obtain the integral equation for |${\hat {\tilde G}}(\tilde r,s|\tilde r_0)$| $\hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0})$ in the following form:

\begin{eqnarray}{\hat {\tilde G}}(\tilde r,s|\tilde r_0) &=& \varphi (\tilde r_0)^{ - 1} \left[ {\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0)\vphantom{\int}\right.\nonumber\\&&\left. -\, s\!\!\int {d{\bf \tilde r}_1 } \;{\hat {\tilde G}}(\tilde r,s|\tilde r_1)d(\tilde r_1){\hat {\tilde G}}_0 (\tilde r_1,s|\tilde r_0) \right]\!.\end{eqnarray}

\begin{matrix} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) & = & ϕ {({\tilde{r}}_{0})}^{- 1} [{\hat{\tilde{G}}}_{0}, (\tilde{r}, s | {\tilde{r}}_{0})) \\ (- s \int d {\tilde{r}}_{1} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{1}) d ({\tilde{r}}_{1}) {\hat{\tilde{G}}}_{0} ({\tilde{r}}_{1}, s | {\tilde{r}}_{0})] . \end{matrix}

(9)

The usual iterative method for solving the Fredholm integral equation of the second kind gives the following series solution⁶:

\begin{eqnarray}{\hat {\tilde G}}(\tilde r,s|\tilde r_0) &=& \varphi (\tilde r_0)^{ - 1} \left[\vphantom{\left.\times \int {d{\bf \tilde r}_1 } \;{\hat {\tilde G}}_0 (\tilde r_2,s|\tilde r_1)\psi (\tilde r_1){\hat {\tilde G}}_0 (\tilde r_1,s|\tilde r_0) + \cdot \cdot \cdot \right.} {\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0)\right.\nonumber\\ [-2pt]&& \ -\, s\int {d{\bf \tilde r}_1 } \;{\hat {\tilde G}}_0 (\tilde r,s|\tilde r_1)\psi (\tilde r_1){\hat {\tilde G}}_0 (\tilde r_1,s|\tilde r_0)\nonumber \\[-2pt]&& \ +\, s^2 \int {d{\bf \tilde r}_2 } \;{\hat {\tilde G}}_0 (\tilde r,s|\tilde r_2)\psi (\tilde r_2)\\[-2pt]&&\left.\times \int {d{\bf \tilde r}_1 } \;{\hat {\tilde G}}_0 (\tilde r_2,s|\tilde r_1)\psi (\tilde r_1){\hat {\tilde G}}_0 (\tilde r_1,s|\tilde r_0) + \cdot \cdot \cdot \right],\nonumber\end{eqnarray}

\begin{matrix} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) & = & ϕ {({\tilde{r}}_{0})}^{- 1} [{\hat{\tilde{G}}}_{0}, (\tilde{r}, s | {\tilde{r}}_{0})) \\ - s \int d {\tilde{r}}_{1} {\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{1}) ψ ({\tilde{r}}_{1}) {\hat{\tilde{G}}}_{0} ({\tilde{r}}_{1}, s | {\tilde{r}}_{0}) \\ + s^{2} \int d {\tilde{r}}_{2} {\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{2}) ψ ({\tilde{r}}_{2}) \\ (\times \int d {\tilde{r}}_{1} {\hat{\tilde{G}}}_{0} ({\tilde{r}}_{2}, s | {\tilde{r}}_{1}) ψ ({\tilde{r}}_{1}) {\hat{\tilde{G}}}_{0} ({\tilde{r}}_{1}, s | {\tilde{r}}_{0}) + \cdot \cdot \cdot], \end{matrix}

(10)

where |$\psi (\tilde r) = d(\tilde r)/\varphi (\tilde r)$| $ψ (\tilde{r}) = d (\tilde{r}) / ϕ (\tilde{r})$ ⁠. However, this solution would be useful only in the small-s limit. For large s, the series will diverge. We can circumvent this difficulty by rewriting Eq. (9) as

\begin{eqnarray}\hskip-3pt{\hat {\tilde G}}(\tilde r,s|\tilde r_0) &=& {{\varphi (\tilde r_0)}}^{-1}\left[ {\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0) - s {\hat {\tilde G}}(\tilde r,s|\tilde r_0) \int \!\! {d{\bf \tilde r}_1 }\right.\nonumber\\&&\times\,\left.{{{\hat {\tilde G}}(\tilde r,s|\tilde r_0)^{-1}}}{{{\hat {\tilde G}}(\tilde r,s|\tilde r_1)}}d(\tilde r_1){\hat {\tilde G}}_0 (\tilde r_1,s|\tilde r_0)\! \right] \!,\nonumber\hskip-6pt\\\end{eqnarray}

\begin{matrix} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) & = & {ϕ ({\tilde{r}}_{0})}^{- 1} [{\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{0}) - s \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) \int d {\tilde{r}}_{1}) \\ \times (\hat{\tilde{G}} {(\tilde{r}, s | {\tilde{r}}_{0})}^{- 1}, \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{1}), d, ({\tilde{r}}_{1}), {\hat{\tilde{G}}}_{0}, ({\tilde{r}}_{1}, s | {\tilde{r}}_{0})], \end{matrix}

(11)

which yields the following renormalized equation,

\begin{eqnarray}{\hat {\tilde G}}(\tilde r,s|\tilde r_0) &=& {{\varphi (\tilde r_0)^{-1}}}{{{\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0)}}\left[ 1 + s{{\varphi (\tilde r_0)^{-1}}}\vphantom{\frac{{{\hat {\tilde G}}(\tilde r,s|\tilde r_1)}}{{{\hat {\tilde G}}(\tilde r,s|\tilde r_0)}}}\!\! \int \!\! {d{\bf \tilde r}_1 } \;\right.\nonumber\\&&\left.\times\,{{{\hat {\tilde G}}(\tilde r,s|\tilde r_0)^{-1}}}{{{\hat {\tilde G}}(\tilde r,s|\tilde r_1)}}d(\tilde r_1){\hat {\tilde G}}_0 (\tilde r_1,s|\tilde r_0) \right]^{ - 1}\!\!\!\!\!\! .\nonumber\\\end{eqnarray}

\begin{matrix} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) & = & ϕ {({\tilde{r}}_{0})}^{- 1} {\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{0}) [1 + s ϕ {({\tilde{r}}_{0})}^{- 1} \int d {\tilde{r}}_{1}) \\ {(\times \hat{\tilde{G}} {(\tilde{r}, s | {\tilde{r}}_{0})}^{- 1} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{1}) d ({\tilde{r}}_{1}) {\hat{\tilde{G}}}_{0} ({\tilde{r}}_{1}, s | {\tilde{r}}_{0})]}^{- 1} . \end{matrix}

(12)

Note that this is still an exact equation. We can then make an approximation given by |${\hat {\tilde G}}(\tilde r,s|\tilde r_1)/{\hat {\tilde G}}(\tilde r,s|\tilde r_0)\break \cong {\hat {\tilde G}}_a (\tilde r,s|\tilde r_1)/{\hat {\tilde G}}_a (\tilde r,s|\tilde r_0)$| $\hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{1}) / \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) ≅ {\hat{\tilde{G}}}_{a} (\tilde{r}, s | {\tilde{r}}_{1}) / {\hat{\tilde{G}}}_{a} (\tilde{r}, s | {\tilde{r}}_{0})$ ⁠, where |${\hat {\tilde G}}_a (\tilde r,s|\tilde r_0)$| ${\hat{\tilde{G}}}_{a} (\tilde{r}, s | {\tilde{r}}_{0})$ is a suitable approximation for |${\hat {\tilde G}}(\tilde r,s|\tilde r_0)$| $\hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0})$ ⁠. This procedure provides a convergent solution even in the s → ∞ limit. The error caused by using the approximate function |${\hat {\tilde G}}_a (\tilde r,s|\tilde r_0)$| ${\hat{\tilde{G}}}_{a} (\tilde{r}, s | {\tilde{r}}_{0})$ would be partially offset because it enters as a ratio.

A simple choice for |${\hat {\tilde G}}_a (\tilde r,s|\tilde r_0)$| ${\hat{\tilde{G}}}_{a} (\tilde{r}, s | {\tilde{r}}_{0})$ is the leading term approximation to |${\hat {\tilde G}}(\tilde r,s|\tilde r_0)$| $\hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0})$ in Eq. (10); that is, |${\hat {\tilde G}}_a (\tilde r,s|\tilde r_0)\break = \varphi (\tilde r_0)^{ - 1} {\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0)$| ${\hat{\tilde{G}}}_{a} (\tilde{r}, s | {\tilde{r}}_{0}) = ϕ {({\tilde{r}}_{0})}^{- 1} {\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{0})$ ⁠. We then have

\begin{eqnarray}{\hat {\tilde G}}(\tilde r,s|\tilde r_0) &=& {{\varphi (\tilde r_0)^{-1}}}{\hat {\tilde G}}_0 (\tilde r,s|\tilde r_0)\bigg[ 1 + s \! \int \! {d{\bf \tilde r}_1 } \;\psi (\tilde r_1)\nonumber\\&& \times{\,} {\hat{\tilde G}_0}(\tilde r,s|\tilde r_0)^{-1}{\hat{\tilde G}_0} (\tilde r,s|\tilde r_1) {\hat{\tilde G}_0} (\tilde r_1,s|\tilde r_0)\bigg]^{-1}\!\!\!.\nonumber\\ \end{eqnarray}

\begin{matrix} \hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0}) & = & ϕ {({\tilde{r}}_{0})}^{- 1} {\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{0}) [1 + s \int d {\tilde{r}}_{1} ψ ({\tilde{r}}_{1}) \\ \times {\hat{\tilde{G}}}_{0} {(\tilde{r}, s | {\tilde{r}}_{0})}^{- 1} {\hat{\tilde{G}}}_{0} (\tilde{r}, s | {\tilde{r}}_{1}) {\hat{\tilde{G}}}_{0} ({\tilde{r}}_{1}, s | {\tilde{r}}_{0})]^{- 1} . \end{matrix}

(13)

This expression for |${\hat {\tilde G}}(\tilde r,s|\tilde r_0)$| $\hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0})$ is exact in the small-s limit and remains quite accurate for larger s values as will be shown below. For better accuracy, one may use the higher order expressions in Eq. (10) for |${\hat {\tilde G}}_a (\tilde r,s|\tilde r_0)$| ${\hat{\tilde{G}}}_{a} (\tilde{r}, s | {\tilde{r}}_{0})$ ⁠.

As an application of the propagator expression given by Eq. (13), we now consider the recombination rate of a geminate pair of reactant molecules. We assume that the initial distance between the reactants is r₀; for the general situation, we may average the final r₀-dependent rate expression over the initial distance distribution. When the initial concentrations of the reactants are very small, we may consider the reaction dynamics of an isolated pair of reactants from the same parent molecule. The effect of the recombination reaction event can be incorporated into the Smoluchowski equation in Eq. (1) in terms of a reaction sink function.⁷ We will assume a delta-function reaction sink, λδ(r − σ)/(4πσ²), which is equivalent to imposing the radiation boundary condition.⁷

Then, the recombination rate of the reactant pair at time t is given by

\begin{equation}R(t\left| {r_0 } \right.) = \int {d{\bf r}} \,\lambda \frac{{\delta (r - \sigma)}}{{4\pi \sigma ^2 }}G_R (r,t\left| {r_0 } \right.) = \lambda G_R (\sigma,t\left| {r_0 } \right.),\end{equation}

R (t |r_{0})) = \int d r λ \frac{δ (r - σ)}{4 π σ^{2}} G_{R} (r, t |r_{0})) = λ G_{R} (σ, t |r_{0})),

(14)

where G_R(r, t|r₀) is the propagator in the presence of the reaction sink. It can be shown that G_R(r, t|r₀) can be related to the reaction-free propagator G(r, t|r₀) in the Laplace transform domain as |$\hat G_R (r,s | {r_0 } ) = \hat G(r,s | {r_0 } )\break - \lambda [ {1 + \lambda \hat G(\sigma,s| \sigma )} ]^{ - 1} \hat G(r,s| \sigma )\hat G(\sigma,s| {r_0 } )$| ${\hat{G}}_{R} (r, s | r_{0}) = \hat{G} (r, s | r_{0}) - λ {[1 + λ \hat{G} (σ, s | σ)]}^{- 1} \hat{G} (r, s | σ) \hat{G} (σ, s | r_{0})$ ⁠. We thus obtain

\begin{equation}\hat R(s\left| {r_0 } \right.) = \frac{{\lambda e^{ - U(\sigma)} \varphi (\tilde r_0){\hat {\tilde G}}(\tilde \sigma,s\left| {\tilde r_0 } \right.)}}{{1 + \lambda e^{ - U(\sigma)} \varphi (\tilde \sigma){\hat {\tilde G}}(\tilde \sigma,s\left| {\tilde \sigma } \right.)}}.\end{equation}\vspace*{-6pt}

\hat{R} (s |r_{0})) = \frac{λ e^{- U (σ)} ϕ ({\tilde{r}}_{0}) \hat{\tilde{G}} (\tilde{σ}, s |{\tilde{r}}_{0}))}{1 + λ e^{- U (σ)} ϕ (\tilde{σ}) \hat{\tilde{G}} (\tilde{σ}, s |\tilde{σ}))} .

(15)

For a perfect absorbing sink with λ → ∞, corresponding to the Smoluchowski absorbing boundary condition, Eq. (15) gives the first encounter time distribution for a pair of particles separated initially by r₀. Using the expression for |${\hat {\tilde G}}(\tilde r,s|\tilde r_0)$| $\hat{\tilde{G}} (\tilde{r}, s | {\tilde{r}}_{0})$ given by Eqs. (7) and (13), we have

\begin{equation}\hat R(s\left| {r_0 } \right.) = \frac{{\tilde \sigma }}{{\tilde r_0 }}e^{ - \zeta (\tilde r_0 - \tilde \sigma)} \left( {\frac{{1 + (\zeta /2)I_d (s|\tilde \sigma)}}{{1 + (\zeta /2)I_d (s|\tilde r_0)}}} \right),\end{equation}

\hat{R} (s |r_{0})) = \frac{\tilde{σ}}{{\tilde{r}}_{0}} e^{- ζ ({\tilde{r}}_{0} - \tilde{σ})} (\frac{1 + (ζ / 2) I_{d} (s | \tilde{σ})}{1 + (ζ / 2) I_{d} (s | {\tilde{r}}_{0})}),

(16)

where

\begin{equation*}I_d (s|\tilde r_0) \equiv \! \int_{\tilde \sigma }^\infty \! {d\tilde r} \,\psi (\tilde r) \! \left[ {e^{ - \zeta (\tilde r - \tilde r_0 + |\tilde r - \tilde r_0 |)} + \frac{{\tilde \sigma \zeta - 1}}{{\tilde \sigma \zeta + 1}}e^{ - 2\zeta (\tilde r - \tilde \sigma)} } \right].\end{equation*}

I_{d} (s | {\tilde{r}}_{0}) \equiv \int_{\tilde{σ}}^{\infty} d \tilde{r} ψ (\tilde{r}) [e^{- ζ (\tilde{r} - {\tilde{r}}_{0} + | \tilde{r} - {\tilde{r}}_{0} |)} + \frac{\tilde{σ} ζ - 1}{\tilde{σ} ζ + 1} e^{- 2 ζ (\tilde{r} - \tilde{σ})}] .

To obtain an explicit expression in the time domain, we now approximate |$\psi (\tilde r)$| $ψ (\tilde{r})$ as a biexponential function:

\begin{equation}y(\tilde r) \equiv {{\psi (\tilde r)}}/{{\psi (\tilde \sigma)}} \cong c_1 e^{ - \alpha _1 (\tilde r - \tilde \sigma)} + c_2 e^{ - \alpha _2 (\tilde r - \tilde \sigma)}.\end{equation}

y (\tilde{r}) \equiv ψ (\tilde{r}) / ψ (\tilde{σ}) ≅ c_{1} e^{- α_{1} (\tilde{r} - \tilde{σ})} + c_{2} e^{- α_{2} (\tilde{r} - \tilde{σ})} .

(17)

The four parameters {c₁, c₂, α₁, α₂} can be determined by solving the following equations:

\begin{equation*}c_1 + c_2 = 1,\end{equation*}\vspace*{-10pt}

c_{1} + c_{2} = 1,

\begin{eqnarray*}\alpha &\equiv& \left. { - {({dy}}/{{d\tilde r}})} \right|_{\tilde r = \tilde \sigma }\nonumber\\&=& \frac{{2h(\sigma)U^\prime(\sigma) - h^\prime (\sigma) - 4\sigma ^{ - 1} h(\sigma) + 4\sigma ^{ - 2} \tilde \sigma e^{U(\sigma)} }}{{\psi (\tilde \sigma) ({{\tilde \sigma }}/{\sigma })^6 h(\sigma)^{ - 1} e^{3U(\sigma)} }}\nonumber\\&=& c_1 \alpha _1 + c_2 \alpha _2,\end{eqnarray*}\vspace*{-12pt}

\begin{matrix} α & \equiv & {(- (d y / d \tilde{r})|}_{\tilde{r} = \tilde{σ}} \\ = & \frac{2 h (σ) U^{'} (σ) - h^{'} (σ) - 4 σ^{- 1} h (σ) + 4 σ^{- 2} \tilde{σ} e^{U (σ)}}{ψ (\tilde{σ}) {(\tilde{σ} / σ)}^{6} h {(σ)}^{- 1} e^{3 U (σ)}} \\ = & c_{1} α_{1} + c_{2} α_{2}, \end{matrix}

\begin{eqnarray*}\beta &\equiv& \int_{\tilde \sigma }^\infty {d\tilde r} \,y(\tilde r) = \frac{1}{{\psi (\tilde \sigma)}}\int_\sigma ^\infty {dr} \,\frac{{r^4 h(r) - \tilde r^4 e^{2U(r)} }}{{r^2 \tilde r^2 h(r)e^{U(r)} }}\nonumber\\&=& ({{c_1 }}/{{\alpha _1 }}) + ({{c_2 }}/{{\alpha _2 }}),\end{eqnarray*}

\begin{matrix} β & \equiv & \int_{\tilde{σ}}^{\infty} d \tilde{r} y (\tilde{r}) = \frac{1}{ψ (\tilde{σ})} \int_{σ}^{\infty} d r \frac{r^{4} h (r) - {\tilde{r}}^{4} e^{2 U (r)}}{r^{2} {\tilde{r}}^{2} h (r) e^{U (r)}} \\ = & (c_{1} / α_{1}) + (c_{2} / α_{2}), \end{matrix}

\begin{eqnarray}\delta &\equiv& \int_{\tilde \sigma }^\infty {d\tilde r} \,y(\tilde r)^2 = \frac{1}{{\psi (\tilde \sigma)^2 }}\int_\sigma ^\infty {dr} \,\frac{{[r^4 h(r) - \tilde r^4 e^{2U(r)}]^2 }}{{r^2 \tilde r^6 h(r)e^{3U(r)} }}\nonumber\\&=& \big({{c_1 ^2 }}/{{2\alpha _1 }}\big) + \big({{c_2 ^2 }}/{{2\alpha _2 }}\big) + [{{2c_1 c_2 }}/({{\alpha _1 + \alpha _2 }})],\end{eqnarray}

\begin{matrix} δ & \equiv & \int_{\tilde{σ}}^{\infty} d \tilde{r} y {(\tilde{r})}^{2} = \frac{1}{ψ {(\tilde{σ})}^{2}} \int_{σ}^{\infty} d r \frac{{[r^{4} h (r) - {\tilde{r}}^{4} e^{2 U (r)}]}^{2}}{r^{2} {\tilde{r}}^{6} h (r) e^{3 U (r)}} \\ = & (c_{1}^{2} / 2 α_{1}) + (c_{2}^{2} / 2 α_{2}) + [2 c_{1} c_{2} / (α_{1} + α_{2})], \end{matrix}

(18)

where |$\psi (\tilde \sigma) = {{[ {\sigma ^4 h(\sigma) - \tilde \sigma ^4 e^{2U(\sigma)} } ]} \mathord{/ {\vphantom {{[ {\sigma ^4 h(\sigma) - \tilde \sigma ^4 e^{2U(\sigma)} } ]} {[ {\tilde \sigma ^4 e^{2U(\sigma)} } ]}}} \kern-\nulldelimiterspace} {[ {\tilde \sigma ^4 e^{2U(\sigma)} } ]}}$| $ψ (\tilde{σ}) = [σ^{4} h (σ) - {\tilde{σ}}^{4} e^{2 U (σ)}] / [{\tilde{σ}}^{4} e^{2 U (σ)}]$ ⁠, h′(σ) = dh(r)/dr|_{r = σ}, and U′(σ) = .dU(r)/dr|_{r = σ}. We have

\begin{eqnarray*}\alpha _i \hspace*{-1pt} &=& \hspace*{-1pt} \frac{{\alpha \beta \hspace*{-1pt} - \hspace*{-1pt} 1 \hspace*{-1pt} + \hspace*{-1pt} ( - 1)^i \sqrt {\left( {\alpha \beta \hspace*{-1pt} - \hspace*{-1pt} 1} \right)^2 \hspace*{-1pt} + \hspace*{-1pt} 4\left( {2\alpha \delta - 1} \right)\left( {2\beta ^{ - 1} \delta - 1} \right)} }}{{2(\beta - 2\delta)}}\nonumber\\&&{\rm (}i = 1,2),\end{eqnarray*}\vspace*{-8pt}

\begin{matrix} α_{i} & = & \frac{α β - 1 + {(- 1)}^{i} \sqrt{{(α β - 1)}^{2} + 4 (2 α δ - 1) (2 β^{- 1} δ - 1)}}{2 (β - 2 δ)} \\ (i = 1, 2), \end{matrix}

\begin{equation}c_1 = ({{\alpha - \alpha _2 }})/({{\alpha _1 - \alpha _2 }}),\quad c_2 = ({{\alpha - \alpha _1 }})/({{\alpha _2 - \alpha _1 }}).\end{equation}

c_{1} = (α - α_{2}) / (α_{1} - α_{2}), c_{2} = (α - α_{1}) / (α_{2} - α_{1}) .

(19)

Then, Eq. (16) becomes

\begin{eqnarray}\hat R(s\left| {r_0 } \right.) \cong \frac{{\tilde \sigma }}{{\tilde r_0 }}e^{ - \zeta (\tilde r_0 - \tilde \sigma)} \left( {\frac{{b_0 + b_1 \zeta + b_2 \zeta ^2 + b_3 \zeta ^3 }}{{a_0 + a_1 \zeta + a_2 \zeta ^2 + a_3 \zeta ^3 + a_4 \zeta ^4 }}} \right),\hspace*{-9pt}\nonumber\\\end{eqnarray}

\begin{matrix} \hat{R} (s |r_{0})) ≅ \frac{\tilde{σ}}{{\tilde{r}}_{0}} e^{- ζ ({\tilde{r}}_{0} - \tilde{σ})} (\frac{b_{0} + b_{1} ζ + b_{2} ζ^{2} + b_{3} ζ^{3}}{a_{0} + a_{1} ζ + a_{2} ζ^{2} + a_{3} ζ^{3} + a_{4} ζ^{4}}), \end{matrix}

(20)

where the various coefficients are given by

\begin{equation*}b_0 \equiv \alpha _1 \alpha _2,\quad b_1 \equiv 2\beta \alpha _1 \alpha _2 + 2\alpha + \alpha _1 \alpha _2 \tilde \sigma,\end{equation*}\vspace*{-14pt}

b_{0} \equiv α_{1} α_{2}, b_{1} \equiv 2 β α_{1} α_{2} + 2 α + α_{1} α_{2} \tilde{σ},

\begin{equation*}b_2 \equiv 4 \hspace*{-1pt} + \hspace*{-1pt} 2\alpha \tilde \sigma \hspace*{-1pt} + \hspace*{-1pt} [2 \hspace*{-1pt} + \hspace*{-1pt} \psi (\tilde \sigma)]\beta \alpha _1 \alpha _2 \tilde \sigma,\quad b_3 \hspace*{-1pt} \equiv \hspace*{-1pt} 2[2 \hspace*{-1pt} + \hspace*{-1pt} \psi (\tilde \sigma)]\tilde \sigma,\end{equation*}\vspace*{-12pt}

b_{2} \equiv 4 + 2 α \tilde{σ} + [2 + ψ (\tilde{σ})] β α_{1} α_{2} \tilde{σ}, b_{3} \equiv 2 [2 + ψ (\tilde{σ})] \tilde{σ},

\begin{eqnarray*}&&\hspace*{48pt}a_0 \equiv b_0,\quad a_1 \equiv b_1,\nonumber\\a_2 \hspace*{-1.5pt} &\equiv& \hspace*{-1.5pt} b_2 \hspace*{-1.5pt} + \hspace*{-1.5pt} \left[ \hspace*{-1.5pt} {\frac{{c_1 \alpha _2 }}{{\alpha _1 }} \hspace*{-1.5pt} \left( \hspace*{-1.5pt} {1 \hspace*{-1.5pt} - \hspace*{-1pt} e^{ - \alpha _1 (\tilde r_0 - \tilde \sigma)} } \right) \hspace*{-1.5pt} + \hspace*{-1.5pt} \frac{{c_2 \alpha _1 }}{{\alpha _2 }} \hspace*{-1pt} \left( \hspace*{-1pt} {1 - e^{ - \alpha _2 (\tilde r_0 - \tilde \sigma)} } \right)} \hspace*{-2pt} \right] \hspace*{-1pt} \psi (\tilde \sigma),\end{eqnarray*}

\begin{matrix} a_{0} \equiv b_{0}, a_{1} \equiv b_{1}, \\ a_{2} & \equiv & b_{2} + [\frac{c_{1} α_{2}}{α_{1}} (1 - e^{- α_{1} ({\tilde{r}}_{0} - \tilde{σ})}) + \frac{c_{2} α_{1}}{α_{2}} (1 - e^{- α_{2} ({\tilde{r}}_{0} - \tilde{σ})})] ψ (\tilde{σ}), \end{matrix}

\begin{eqnarray*}a_3 &\equiv& [4 + \psi (\tilde \sigma)]\tilde \sigma + \beta \psi (\tilde \sigma)(2 + \alpha \tilde \sigma + \beta \alpha _1 \alpha _2 \tilde \sigma)\nonumber\\&&-\,\big[ ({{c_1 }}/{{\alpha _1 }})\left( {2 \hspace*{-1pt} + \hspace*{-1pt} \alpha _2 \tilde \sigma } \right)e^{ - \alpha _1 (\tilde r_0 - \tilde \sigma)} \hspace*{-1pt} \nonumber\\&&+\, ({{c_2 }}/{{\alpha _2 }})\left( {2 \hspace*{-1pt} + \hspace*{-1pt} \alpha _1 \tilde \sigma } \right)e^{ - \alpha _2 (\tilde r_0 - \tilde \sigma)} \big]\psi (\tilde \sigma),\end{eqnarray*}

\begin{matrix} a_{3} & \equiv & [4 + ψ (\tilde{σ})] \tilde{σ} + β ψ (\tilde{σ}) (2 + α \tilde{σ} + β α_{1} α_{2} \tilde{σ}) \\ - [(c_{1} / α_{1}) (2 + α_{2} \tilde{σ}) e^{- α_{1} ({\tilde{r}}_{0} - \tilde{σ})} \\ + (c_{2} / α_{2}) (2 + α_{1} \tilde{σ}) e^{- α_{2} ({\tilde{r}}_{0} - \tilde{σ})}] ψ (\tilde{σ}), \end{matrix}

\begin{equation}a_4 \equiv 2\psi (\tilde \sigma)\tilde \sigma \left( {\beta - \frac{{c_1 }}{{\alpha _1 }}e^{ - \alpha _1 (\tilde r_0 - \tilde \sigma)} - \frac{{c_2 }}{{\alpha _2 }}e^{ - \alpha _2 (\tilde r_0 - \tilde \sigma)} } \right).\end{equation}

a_{4} \equiv 2 ψ (\tilde{σ}) \tilde{σ} (β - \frac{c_{1}}{α_{1}} e^{- α_{1} ({\tilde{r}}_{0} - \tilde{σ})} - \frac{c_{2}}{α_{2}} e^{- α_{2} ({\tilde{r}}_{0} - \tilde{σ})}) .

(21)

By first making a small-s approximation to the ζ³ and ζ⁴ terms in the denominator of Eq. (20) and then making an error-compensating large-t approximation to the resulting time-domain expression, we obtain

\begin{equation}\hspace*{-18pt}R_{DC} (t\left| {r_0 } \right.) \cong \frac{{\tilde \sigma }}{{\tilde r_0 }}\left( {\frac{{\tau _g }}{{\pi t^3 }}} \right)^{1/2} \exp \left( { - \frac{{\tau _g }}{t}} \right)\left( {\frac{{1 + I_1 (t)}}{{1 + I_2 (t)}}} \right).\end{equation}

R_{D C} (t |r_{0})) ≅ \frac{\tilde{σ}}{{\tilde{r}}_{0}} {(\frac{τ_{g}}{π t^{3}})}^{1 / 2} \exp (- \frac{τ_{g}}{t}) (\frac{1 + I_{1} (t)}{1 + I_{2} (t)}) .

(22)

Here, |$\tau _g [ \equiv (\tilde r_0 - \tilde \sigma)^2 /(4D)]$| $τ_{g} [\equiv {({\tilde{r}}_{0} - \tilde{σ})}^{2} / (4 D)]$ roughly measures the encounter time scale for the reactant pair and

\begin{eqnarray*}I_1 (t) &\equiv& \frac{2}{{\tilde r_0 - \tilde \sigma }}\sum^2_{i=1}(-1)^i \frac{{b_2 - a_2 - b_3 \gamma _ i }}{{a_2 \left( {\gamma _ 1 - \gamma _ 2 } \right)}}M_1 (t|\gamma _ i),\end{eqnarray*}

\begin{matrix} I_{1} (t) & \equiv & \frac{2}{{\tilde{r}}_{0} - \tilde{σ}} \sum_{i = 1}^{2} {(- 1)}^{i} \frac{b_{2} - a_{2} - b_{3} γ_{i}}{a_{2} (γ_{1} - γ_{2})} M_{1} (t | γ_{i}), \end{matrix}

\begin{eqnarray*}M_1 (t|\gamma) &\equiv& \gamma ^2 Dt - ({1}/{2}) - \gamma \sqrt{D\tau _g} + (\tau _g /t)\nonumber\\&& - \gamma ^3 \sqrt {\pi D^3 t^3 } \Omega \left( {\sqrt {{{\tau _g }}/{t}} + \gamma \sqrt {Dt} } \right),\end{eqnarray*}

\begin{matrix} M_{1} (t | γ) & \equiv & γ^{2} D t - (1 / 2) - γ \sqrt{D τ_{g}} + (τ_{g} / t) \\ - γ^{3} \sqrt{π D^{3} t^{3}} Ω (\sqrt{τ_{g} / t} + γ \sqrt{D t}), \end{matrix}

\begin{eqnarray*}I_2 (t) &\equiv& \frac{2}{{\tilde r_0 - \tilde \sigma}} \sum^2_{i=1}(-1)^i \frac{{a_3 - a_4 \gamma _ i }}{{a_2 \left( {\gamma _ 1 - \gamma _ 2 } \right)}}M_2 (t|\gamma _ i),\end{eqnarray*}

\begin{matrix} I_{2} (t) & \equiv & \frac{2}{{\tilde{r}}_{0} - \tilde{σ}} \sum_{i = 1}^{2} {(- 1)}^{i} \frac{a_{3} - a_{4} γ_{i}}{a_{2} (γ_{1} - γ_{2})} M_{2} (t | γ_{i}), \end{matrix}

\begin{eqnarray*}M_2 (t|\gamma) &\equiv& \sqrt {{1}/{{Dt}}} \left\{ \big[ ( {{{\tau _g }}/{t}) + \gamma ^2 Dt} \big] \big( {\sqrt {{{\tau _g }}/{t}}} \right.\nonumber\\&&\left.-\, \gamma \sqrt {Dt} \big) - ({3}/{2})\sqrt {{{\tau _g }}/{t}} + {{(\gamma/2) \sqrt {Dt} }} \right \}\nonumber\\&&\,+\, \gamma ^4 \sqrt {\pi D^3 t^3 } \Omega \big( {\sqrt {{{\tau _g }}/{t}} + \gamma \sqrt {Dt} } \big),\end{eqnarray*}

\begin{matrix} M_{2} (t | γ) & \equiv & \sqrt{1 / D t} \{[(τ_{g} / t) + γ^{2} D t] (\sqrt{τ_{g} / t}) \\ (- γ \sqrt{D t}) - (3 / 2) \sqrt{τ_{g} / t} + (γ / 2) \sqrt{D t}\} \\ + γ^{4} \sqrt{π D^{3} t^{3}} Ω (\sqrt{τ_{g} / t} + γ \sqrt{D t}), \end{matrix}

where Ω(x) ≡ exp (x²)erfc(x), |$\gamma _ i = [a_1\break + (-1)^i\sqrt {a{\smash{_1^2}} - 4a{\smash{_2}} a{\smash{_0}} }]/ {(2a_2)}$| $γ_{i} = [a_{1} + {(- 1)}^{i} \sqrt{a_{1}^{2} - 4 a_{2} a_{0}}] / (2 a_{2})$ ⁠.

Another important quantity in the chemical kinetics is the time-dependent rate coefficient for diffusion-influenced bimolecular reactions.² In the low reactant concentration limit, it can be expressed in terms of the reaction-free propagator as

\begin{equation}s\hat k_f (s) = [ {1 + k_f (0)\varphi (\tilde \sigma){\hat {\tilde G}}(\tilde \sigma,s\left| {\tilde \sigma } \right.)} ]^{ - 1} k_f (0).\end{equation}

s {\hat{k}}_{f} (s) = {[1 + k_{f} (0) ϕ (\tilde{σ}) \hat{\tilde{G}} (\tilde{σ}, s |\tilde{σ}))]}^{- 1} k_{f} (0) .

(23)

Here, k_f(0) [ = λe^−U(σ)] is the rate constant at t = 0 when the reactants are distributed in equilibrium. In the diffusion-controlled limit (λ → ∞), we have from Eqs. (7), (13), and (23)

\begin{equation}\hat k_f (s) = k_f (\infty)\left[ {\frac{{1 + \tilde \sigma \zeta }}{s} + \frac{{\tilde \sigma }}{D}\int_{\tilde \sigma }^\infty {d\tilde r} \,\psi (\tilde r)e^{ - 2\zeta (\tilde r - \tilde \sigma)} } \right],\end{equation}

{\hat{k}}_{f} (s) = k_{f} (\infty) [\frac{1 + \tilde{σ} ζ}{s} + \frac{\tilde{σ}}{D} \int_{\tilde{σ}}^{\infty} d \tilde{r} ψ (\tilde{r}) e^{- 2 ζ (\tilde{r} - \tilde{σ})}],

(24)

whose inverse Laplace transform gives

\begin{eqnarray}k_f (t) &=& k_f (\infty)\left[ 1 + \frac{{\tilde \sigma }}{{(\pi Dt)^{1/2} }} + \frac{{\pi \tilde \sigma }}{{(\pi Dt)^{3/2} }}\int_\sigma ^\infty{dr} \right.\nonumber\\&&\times\,\left.\frac{{1 - \varphi (r)}}{{[h(r)\varphi (r)]^{1/2} }}(\tilde r - \tilde \sigma){\rm exp}\left[ { - \frac{{(\tilde r - \tilde \sigma)^2 }}{{Dt}}} \right] \right],\nonumber\\\end{eqnarray}

\begin{matrix} k_{f} (t) & = & k_{f} (\infty) [1 + \frac{\tilde{σ}}{{(π D t)}^{1 / 2}} + \frac{π \tilde{σ}}{{(π D t)}^{3 / 2}} \int_{σ}^{\infty} d r) \\ \times (\frac{1 - ϕ (r)}{{[h (r) ϕ (r)]}^{1 / 2}}, (\tilde{r} - \tilde{σ}), \exp, [- \frac{{(\tilde{r} - \tilde{σ})}^{2}}{D t}]], \end{matrix}

(25)

With the approximation to |$\psi (\tilde r)$| $ψ (\tilde{r})$ given by Eq. (17) the radial integral can be evaluated approximately and thus we have

\begin{eqnarray}k_f (t) &=& k_f (\infty)\Bigg\{\vphantom{\left.\left. - \frac{{c_1 \alpha _1 }}{4}\Omega \left( {\frac{{\alpha _1 \sqrt {Dt} }}{2}} \right) - \frac{{c_2 \alpha _2 }}{4}\Omega \left( {\frac{{\alpha _2 \sqrt {Dt} }}{2}} \right) \right] \right.} 1 + \frac{{\tilde \sigma }}{{(\pi Dt)^{1/2} }} + \tilde \sigma \psi (\tilde \sigma)\Bigg[\vphantom{\left.\left. - \frac{{c_1 \alpha _1 }}{4}\Omega \left( {\frac{{\alpha _1 \sqrt {Dt} }}{2}} \right) - \frac{{c_2 \alpha _2 }}{4}\Omega \left( {\frac{{\alpha _2 \sqrt {Dt} }}{2}} \right) \right] \right.} \frac{1}{{2(\pi Dt)^{1/2} }}\qquad\nonumber\\&& -\, \frac{{c_1 \alpha _1 }}{4}\Omega \Bigg( \!{\frac{{\alpha _1 \sqrt {Dt} }}{2}} \Bigg) {-\,} \frac{{c_2 \alpha _2 }}{4}\Omega \Bigg( {\frac{{\alpha _2 \sqrt {Dt} }}{2}} \Bigg)\! \Bigg] \Bigg\}.\end{eqnarray}

\begin{matrix} k_{f} (t) & = & k_{f} (\infty) {1 + \frac{\tilde{σ}}{{(π D t)}^{1 / 2}} + \tilde{σ} ψ (\tilde{σ}) [\frac{1}{2 {(π D t)}^{1 / 2}} \\ - \frac{c_{1} α_{1}}{4} Ω (\frac{α_{1} \sqrt{D t}}{2}) - \frac{c_{2} α_{2}}{4} Ω (\frac{α_{2} \sqrt{D t}}{2})]} . \end{matrix}

(26)

We have checked the accuracy of the rate expressions given by Eqs. (16), (22), (25), and (26) by comparison with the exact results obtained by solving the Smoluchowski equation numerically. Figures 1 and 2 display the effects of the Coulomb interaction, U(r) = r_C/r, and the hydrodynamic interaction, modeled as h(r) = 1 − 3σ/4r, on the first encounter time distribution, respectively. In all calculations, we are using the dimensionless variables: length in σ, time in σ²/D, and energy in k_BT. r_C is the Onsager distance given by r_C/σ = z₁z₂e²/(4πε₀ε_rσk_BT) where z₁e and z₂e are the charges of the reactants, ε₀ is the permittivity of vacuum, and ε_r is the relative permittivity of the medium.

FIG. 1.

View large Download slide

Effects of the Coulomb interaction on the first encounter time distribution. The parameter values are r_C/σ = −4, r₀/σ = 8, and h(r) = 1.

FIG. 2.

View large Download slide

Effects of the hydrodynamic interaction on the first encounter time distribution. The parameter values are r_C/σ = 0, r₀/σ = 8, and h(r) = 1 − 3σ/4r.

In Figs. 1 and 2, the exact numerical results are represented by filled squares. The results obtained from Eq. (16) by numerical inverse Laplace transformation are represented by the solid curves, while the results obtained from the more approximate time-domain expression in Eq. (22) are represented by dotted curves. We also compare the results of the rate expression given by Pedersen and Sibani [Eq. (39) of Ref. 8]. We see that Eq. (16) provides almost exact results for the whole time region.

Figures 3 and 4 display the effects of the Coulomb and the hydrodynamic interactions, respectively, on the bimolecular rate coefficient in the diffusion-controlled limit. The exact numerical results are represented by filled squares, and the results obtained from Eqs. (25) and (26) are represented by the solid and dotted curves, respectively. As for the effects of the interaction potential on k_f(t), Dudko and Szabo constructed a quite simple and accurate time-domain rate expression by the Pade approximation using two exact terms in both the short- and long-time expansions of k_f(t).⁹ The result of their rate expression is represented by the dashed curve in Fig. 3. The rate expression given by Pedersen and Sibani [Eq. (58) of Ref. 8] is just that given by the first two terms of rate expression in Eq. (25). Again, we see that our rate expressions provide very accurate results for the whole time region.

FIG. 3.

View large Download slide

Effects of the Coulomb interaction on the bimolecular rate coefficient in the diffusion-controlled limit, r_C/σ = −4 and h(r) = 1.

FIG. 4.

View large Download slide

Effects of the hydrodynamic interaction on the bimolecular rate coefficient in the diffusion-controlled limit, U(r) = 0 and h(r) = 1 − 3σ/4r.

In summary, we have obtained a very accurate expression of the propagator for diffusive dynamics of a pair of particles interacting via an arbitrary central potential and hydrodynamic interaction. The derivation was based on a new method for solving the linear integral equations. We believe that the procedure presented by Eqs. (9) through (13) can be applied to diverse problems involving linear integral equations; it would provide a convergent solution even in the case where the direct iterative solution method leads to a divergent perturbation series. By using the propagator expression, we have obtained accurate analytic expressions for the first encounter time distribution and the bimolecular rate coefficient that are applicable for the whole time region.

This work was supported by the grants from National Research Foundation (NRF), funded by the Korean Government (Grant Nos. 2009–0074693, 2010–0001631, and KRF-C00180).

REFERENCES

1.

U.

Balucani

and

M.

Zoppi

,

Dynamics of the Liquid State

(

Oxford University Press

,

Oxford

,

1994

).

Google Scholar

2.

S. A.

Rice

, in

Diffusion-Limited Reactions, Comprehensive Chemical Kinetics

, edited by

C. H.

Bamford

,

C. F. H.

Tipper

, and

R. G.

Compton

(

Elsevier

,

Amsterdam

,

1985

), Vol.

25

.

Google Scholar

3.

N.

Agmon

,

Phys. Rev. Lett.

67

,

1366

(

1991

).

https://doi.org/10.1103/PhysRevLett.67.1366

Google Scholar

Crossref

PubMed

4.

K. M.

Hong

and

J.

Noolandi

,

J. Chem. Phys.

68

,

5163

(

1978

).

https://doi.org/10.1063/1.435636

Google Scholar

Crossref

5.

M. R.

Flannery

,

Phys. Rev. Lett.

47

,

163

(

1981

).

https://doi.org/10.1103/PhysRevLett.47.163

Google Scholar

Crossref

6.

A. D.

Polyanin

and

A. V.

Manzhirov

,

Handbook of Integral Equations

, 2nd ed. (

Chapman and Hall/CRC

,

Boca Raton

,

2008

).

Google Scholar

Crossref

7.

G.

Wilemski

and

M.

Fixman

,

J. Chem. Phys.

58

,

4009

(

1973

).

https://doi.org/10.1063/1.1679757

Google Scholar

Crossref

8.

J. B.

Pedersen

and

P.

Sibani

,

J. Chem. Phys.

75

,

5368

(

1981

).

https://doi.org/10.1063/1.441981

Google Scholar

Crossref

9.

O. K.

Dudko

and

A.

Szabo

,

J. Phys. Chem. B

109

,

5891

(

2005

).

https://doi.org/10.1021/jp044433q

Google Scholar

Crossref

PubMed

2011

American Institute of Physics

Communication: Propagator for diffusive dynamics of an interacting molecular pair

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Communication: Propagator for diffusive dynamics of an interacting molecular pair

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only