# Convection-diffusion-reaction problem

I seek to solve to the system

$$\frac{\partial \phi_{a}}{\partial t} = D_{a} \frac{\partial^{2} \phi_{a}}{\partial x^{2}} - v_{a} \frac{\partial \phi_{a}}{\partial x} + \mathfrak{K}_{b}\phi_{b}$$

$$\frac{\partial \phi_{b}}{\partial t} = D_{b} \frac{\partial^{2} \phi_{b}}{\partial x^{2}} - v_{b} \frac{\partial \phi_{b}}{\partial x} - \mathfrak{K}_{b}\phi_{b}$$

in unbounded space, given $\phi_{n} (x,0) = \phi_{n,0}\delta(x)$ where $\phi_{a,0}$, $\phi_{b,0}$, $D_{a}$, $D_{b}$, $v_{a}$,$v_{b}$, $\mathfrak{K}_{b}$ $\in \Re_{+}$ are constants

The solution for $\phi_{b}$ is well known:

$$\phi_{b} = \phi_{b,0}\frac{e^{-\frac{(x-v_{b}t)^2}{4D_{b} t} - \mathfrak{K}_{b} t }}{\sqrt{4 \pi D_{b} t}}$$

To obtain an expression $\phi_a$ one may take a Fourier-Laplace transform:

$$s \hat{\phi}_{a}(k,s) = \phi_{a,0} - (D_{a}k^{2} + iv_{a}k)\hat{\phi}_{a} + \mathfrak{K}_{b} \hat{\phi}_{b}$$

$$s \hat{\phi}_{b}(k,s) = \phi_{b,0} - (D_{b}k^{2} + iv_{b}k + \mathfrak{K}_{b} ) \hat{\phi}_{b}$$

defining $\omega _{a} = D_{a}k^{2} + iv_{a}k$ and $\omega _{b} = D_{b}k^{2} + iv_{b}k + \mathfrak{K}_{b}$, solving the linear system and taking the inverse L.T. we arrive at

$$\tilde{\phi}_{a}(k,t) = \phi_{a,0}e^{-\omega _{a} t} + \mathfrak{K}_{b} \phi_{b,0} \left( \frac{e^{-\omega _{a} t} -e^{-\omega _{b} t}}{\omega_{b}-\omega_{a}} \right)$$

Taking the inverse F.T. leads to the solution:

$$\phi_{a} (x,t) = \phi_{a,0}\frac{e^{-\frac{(x-v_{a}t)^2}{4D_{a} t}}}{\sqrt{4 \pi D_{a} t}} + \frac{\phi_{b,0} \mathfrak{K} _{b}}{2\pi} \int{ \left( \frac{e^{-\omega _{a} t +ikx} -e^{-\omega _{b} t +ikx}}{\omega_{b}-\omega_{a}} \right) dk }$$

To adress $\int{ \mathrm{exp}{(-\omega _{a} t +ikx)} \over {\omega_{b}-\omega_{a}} }dk$ in the above expression, first note:

$$\frac{1}{\omega_{b}-\omega_{a}} = \frac{-i}{(D_b - D_a)(R_1 - R_2)} \left( \frac{1}{k - iR_1} - \frac{1}{k-iR_2} \right)$$

$$R_{1,2}=\frac{1}{2}\left( \frac{v_b - v_a}{D_b - D_a} \pm \sqrt{\left( \frac{v_b - v_a}{D_b - D_a} \right)^{2}+4\frac{\mathfrak{K}_{b}}{D_b - D_a}}\right)$$

distinguishing two cases, $D_a<D_b$ and $D_a>D_b$. Furthermore,

$$-\omega_a t + ikx = \frac{-(x-v_a t)^2}{4D_a t} + \left( k-i\frac{x-v_a t}{2 D_a t} \right) ^2$$

Setting $\kappa_a = \frac{x-v_a t}{2 D_a t}$, $\xi_{1a} = \frac{R_1 - \kappa_a}{\sqrt{D_a t}}$ and $u=\sqrt{D_a t}(k-i\kappa_a)$ we find

$$\int{ \frac{e^{-\omega_a t + ikx}}{k-R_1} dk}= e^{\frac{-(x-v_a t)^2}{4D_a t}} \int{ \frac{e^{-u^2}}{u-i\xi_{1a}}du}= \sqrt{\pi} e^{\frac{-(x-v_a t)^2}{4D_a t}} Z(i\xi_{1a})$$

where $Z(z)$ is the plasma dispersion function. It is known that $Z(z)$ has an alternative representation (for all $\Im \{z\} \neq 0$):

$$Z(z) = i\sqrt{\pi}e^{-z^2}\mathrm{erfc}(-iz)$$

When $D_a<D_b$ $\xi_{1a}$ is real, which leads to the following expression for $\phi_a$:

$$\phi_{a} (x,t) = \phi_{a,0}\frac{e^{-\frac{(x-v_{a}t)^2}{4D_{a} t}}}{\sqrt{4 \pi D_{a} t}} + \frac{\phi_{b,0} \mathfrak{K} _{b}}{(D_b - D_a)(R_1 - R_2)\sqrt{4\pi}} \left( e^{-\frac{(x-v_{a}t)^2}{4D_{a} t}} \left( e^{\xi_{1a} ^{2}} \mathrm{erfc}(\xi_{1a})-e^{\xi_{2a} ^{2}} \mathrm{erfc}(\xi_{2a})\right) -e^{-\frac{(x-v_{b}t)^2}{4D_{b} t} -\mathfrak{K}_{b} t } \left( e^{\xi_{1b} ^{2}} \mathrm{erfc}(\xi_{1b})-e^{\xi_{2b} ^{2}} \mathrm{erfc}(\xi_{2b})\right) \right)$$

I cannot find a major fault with this derivation; however, when I plot the above against time, I find that the density is very low density between the two main peaks (the principal solutions for $\phi_a$ and $\phi_b$) especially for low values of $\mathfrak{K} _{b}$. Can anyone see any obvious faults with this method?

When $R_{1,2}$ are not real (i.e. $D_a>D_b$ for sufficiently large $\mathfrak{K} _{b}$) then one can show that $e^{z^2}\mathrm{erfc}(z)-e^{{z^*}^2}\mathrm{erfc}(z^*)$ is purely imaginary, leading to a real-valued expression for $\phi_a$, as we would expect; however that expression also becomes oscillatory, which bearing in mind the physical process this equation should model, isn't exactly sensible.

I would really appreciate some help on this probelm. Have I made any obvious mistakes in this derivation?

• I suspect the the case where $R_{1,2}$ are not real may be treated by decomposing the Faddeeva function to Voigt functions an taking advantage of the fact that $R_1$ and $R_2$ are complex conjugates; and hopefully the imaginary part will vanish. But I haven't explored this option in detail yet. – JMK Jun 2 '14 at 5:10