This report performs weak inhomogeneity expansion for the continuous Gaussian chain model of a homopolymer in an external field. The aim is to understand where does the Debye function come from. This technique constitutes the main parts of the random phase approximation (RPA). By extending this technique from single-chain formulation to many-chain formulation and regarding the potential field arising from interaction with other chains in the same system as an external field, RPA is equivalent to the weak inhomogeneity expansion.

## Continuous Gaussian Chain

The normalized single partition function for a continuous Gaussian chain under an external field is given by1

$$$Q[w] = \frac{1}{V}\int d\vx\; q(\vx, 1; [w])$$\label{eq:Q}$

where the physical length is scaled by the radius of gyration of a non-perturb Gaussian chain $R_g=\sqrt{Nb^2/6}$. The volume $V$ is also dimensionless scaled by $R_g^3$.

The propagator $q$ in the above equation satisfies the Fokker-Planck equation, also known as the modified diffusion equation (MDE)

$$$\der{q(\vx, s; [w])}{s} = \nabla^2 q(\vx, s; [w]) - w(\vx)q(\vx, s; [w])$$\label{eq:MDE}$

subject to the initial condition

$$$q(\vx, 0; [w]) = 1$$$

Note that to write the MDE in the form as in eq. \eqref{eq:MDE}, $w(\vx)$ is actually $N$ times the external potential field. Assuming the external field is $W(\vx)$, then $w(\vx)=NW(\vx)$

For a given external field $w(\vx)$, we can obtain the propagator by solving eq. \eqref{eq:MDE}.

## Weak Inhomogeneity Expansion

In general, it is impossible to find an analytic solution to the MDE with a general $w(\vx)$. However, a particularly perturbation expansion can be derived when the applied potential field $w(\vx)$ has inhomogeneities that are weak in amplitude. To define such a situation, we introduce the volume average of the potential

$$$w_0 \equiv \frac{1}{V}\int d\vx\; w(\vx)$$\label{eq:w0-w}$

and re-express $w(\vx)$ according to

$$$w(\vx) = w_0 + \epsilon\omega(\vx)$$\label{eq:w-omega}$

which serves to define the inhomogeneous part of the field, $\epsilon\omega(\vx)$. For weak inhomogeneities, a small parameter $\epsilon$ ($\abs{\epsilon}\ll 1$) describes their characteristic amplitude. For the continuous Gaussian chain model of a homopolymer, the MDE and initial condition become

$$$\der{q(\vx, s)}{s} = \nabla^2 q(\vx, s) - w_0q(\vx, s) - \epsilon\omega(\vx)q(\vx, s)$$\label{eq:MDE-w0}$ $$$q(\vx, 0) = 1$$$

where the functional dependence of q on $w(\vx)$ has been suppressed in our notation. The term proportional to $w_0$ on the right-hand side of eq.\eqref{eq:MDE-w0} can be removed by the substitution

$$$q(\vx, s) = e^{-w_0s}p(\vx, s)$$\label{eq:q-p}$

$$$\der{p(\vx, s)}{s} = \nabla^2 p(\vx, s) - \epsilon\omega(\vx)p(\vx, s)$$\label{eq:MDE-omega}$ $$$p(\vx, 0) = 1$$$

A weak inhomogeneity expansion can be developed by assuming that $p(\vx, s)$ can be expressed as

$$$p(\vx, s) \sim \sum_{j=0}^{\infty} \epsilon^j p^{(j)}(\vx, s)$$\label{eq:p-expansion}$

where the $p^{(j)}(\vx, s)$ are independent of $\epsilon$. In eq.\eqref{eq:p-expansion} we adopt the conventional notation $\sim$ to indicate an asymptotic expansion. As such, the infinite series on the right-hand side may be either convergent or divergent. Even when it does not converge, eq.\eqref{eq:p-expansion} can still be useful in truncated form for approximating $p(\vx, s)$ at sufficiently small $\epsilon$.

The $p^{(j)}$ are calculated by inserting eq.\eqref{eq:p-expansion} into eq.\eqref{eq:MDE-omega} and equating terms order by order in $\epsilon$.

### Zero-th Order Solution

At leading order, $O(\epsilon^0)$, we have

$$$\der{p^{(0)}(\vx, s)}{s} = \nabla^2 p^{(0)}(\vx, s)$$\label{eq:MDE-p0}$ $$$p^{(0)}(\vx, 0) = 1$$$

which has the trivial solution

$$$p^{(0)}(\vx, s) = 1$$$

### First Order Solution

At $O(\epsilon^1)$, the corresponding equations (see Appendix A for derivation) are

$$$\der{p^{(1)}(\vx, s)}{s} = \nabla^2 p^{(1)}(\vx, s) - \omega(\vx)p^{(0)}(\vx, s)$$\label{eq:MDE-p1}$ $$$p^{(1)}(\vx, 0) = 0$$$

Provided the system under consideration is unbounded or subject to periodic boundary conditions, this initial value problem is most easily solved by means of spatial Fourier transforms. Defining Fourier transforms in accordance with

$$$\hat{f}(\vk) \equiv \int d\vx\; e^{-i\vk\cdot\vx}f(\vx)$$$

and assuming that the Fourier transform of $\omega(\vx)$ exists, denoted by $\hat{\omega}(\vk)$, one finds that (See Appendix B for derivation)

$$$\hat{p}^{(1)}(\vk, s) = -\hat{h}_2(\vk, s)\hat{\omega}(\vk)$$\label{eq:p1}$

where the carets denote Fourier-transformed quantities and

$$$\hat{h}_2(\vk, s) \equiv \frac{1}{k^2}\left( 1 - e^{-k^2s} \right)$$$

### Second Order Solution

At $O(\epsilon^2)$, the corresponding equations (see Appendix C for derivation) are

$$$\der{p^{(2)}(\vx, s)}{s} = \lap p^{(2)}(\vx, s) - \omega(\vx)p^{(1)}(\vx, s)$$\label{eq:MDE-p2}$ $$$p^{(2)}(\vx, 0) = 0$$$

A similar procedure leads to (see Appendix D for derivation)

$$$\hat{p}^{(2)}(\vk, s) = \frac{1}{V}\sum_{\vk'} \hat{h}_3(\vk, \vk', s)\hat{\omega}(\vk-\vk')\hat{\omega}(\vk')$$\label{eq:p2}$

with

$$$\hat{h}_3(\vk, \vk', s) = \frac{1}{k^2\abs{\vk-\vk'}^2}\left[ 1 - e^{-k^2s} -\frac{k^2}{k^2-\abs{\vk-\vk'}^2} \left( e^{-\abs{\vk-\vk'}^2s} - e^{-k^2s} \right) \right]$$$

### Expansion of $Q$

We can expand the normalized single partition function $Q$ based on above perturbation expansion of the propagator $q$. The propagator being expanded to $O(\epsilon^2)$ is

$$$\begin{split} q(\vx, s) &= e^{-w_0s}p(\vx, s) \\ &\sim \sum_{j=0}^{2} \epsilon^j p^{(j)}(\vx, s) \\ &= e^{-w_0s}\left[ p^{(0)}(\vx, s) + \epsilon p^{(1)}(\vx, s) + \epsilon^2 p^{(2)}(\vx, s) \right] \\ &= e^{-w_0s}\left[ 1 + \epsilon p^{(1)}(\vx, s) + \epsilon^2 p^{(2)}(\vx, s) \right] \end{split}$$$

Therefore we can expand $Q$ as

$$$\begin{split} Q[w] &= \frac{1}{V}\int d\vx\; q(\vx, 1) \\ &\sim \frac{1}{V}\int d\vx\; e^{-w_0}\left[ 1 + \epsilon p^{(1)}(\vx, 1) + \epsilon^2 p^{(2)}(\vx, 1) \right] \\ &= e^{-w_0}\left[ \frac{1}{V}\int d\vx\; + \frac{1}{V}\int d\vx\;\epsilon p^{(1)}(\vx, 1) + \frac{1}{V}\int d\vx\;\epsilon^2 p^{(2)}(\vx, 1) \right] \\ &= e^{-w_0}\left[ 1 + \frac{\epsilon}{V}\int d\vx\; p^{(1)}(\vx, 1) + \frac{\epsilon^2}{V}\int d\vx\; p^{(2)}(\vx, 1) \right] \end{split}$$\label{eq:Q-expan-v1}$

It is more convenient to write the integrals of $p^{(j)}$ in Fourier space, which can be expressed as

$$$\begin{split} \int d\vx\; p^{(1)}(\vx, 1) &= \int d\vx\; \frac{1}{(2\pi)^3}\int d\vk\; \hat{p}^{(1)}(\vk, 1) e^{i\vk\cdot\vx} \\ &= \int d\vk\; \hat{p}^{(1)}(\vk, 1) \left[\frac{1}{(2\pi)^3}\int d\vx\; e^{i\vk\cdot\vx} \right] \\ &= \int d\vk\; \hat{p}^{(1)}(\vk, 1) \delta(\vk) \\ &= \hat{p}^{(1)}(\vzero, 1) \end{split}$$$

and

$$$\begin{split} \int d\vx\; p^{(2)}(\vx, 1) &= \int d\vx\; \frac{1}{(2\pi)^3}\int d\vk\; \hat{p}^{(2)}(\vk, 1) e^{i\vk\cdot\vx} \\ &= \int d\vk\; \hat{p}^{(2)}(\vk, 1) \left[\frac{1}{(2\pi)^3}\int d\vx\; e^{i\vk\cdot\vx} \right] \\ &= \int d\vk\; \hat{p}^{(2)}(\vk, 1) \delta(\vk) \\ &= \hat{p}^{(2)}(\vzero, 1) \end{split}$$$

Now we can express $Q$ in the Fourier space as

$$$Q[w] \sim e^{-w_0}\left[ 1 + \frac{\epsilon}{V}\hat{p}^{(1)}(\vzero, 1) + \frac{\epsilon^2}{V}\hat{p}^{(2)}(\vzero, 1) \right]$$\label{eq:Q-expan-v2}$

The expansion of $Q$ can be further simplified. Firstly, we know that the volume average of the fluctuation of the potential field is 0 because

$$$\begin{split} \frac{1}{V}\int d\vx\; \omega(\vx) &= \frac{1}{V}\int d\vx\; \left[ w(\vx) - w_0 \right] \\ &= \frac{1}{V}\int d\vx\; w(\vx) - \frac{1}{V}\int d\vx\; w_0 \\ &= w_0 - \frac{1}{V} w_0 V \\ &= 0 \end{split}$$$

The value of the Fourier transform of the potential field at zero wave number is equal to this average because

$$$\begin{split} \int d\vx\; \omega(\vx) &= \int d\vx\; \frac{1}{(2\pi)^3}\int d\vk\; \hat{\omega}(\vk)e^{i\vk\cdot\vx} \\ &= \int d\vk\; \frac{1}{(2\pi)^3}\int d\vx\; \hat{\omega}(\vk)e^{i\vk\cdot\vx} \\ &= \int d\vk\; \hat{\omega}(\vk) \left[ \frac{1}{(2\pi)^3}\int d\vx\; e^{i\vk\cdot\vx}\right] \\ &= \int d\vk\; \hat{\omega}(\vk) \delta(\vk) \\ &= \hat{\omega}(\vzero) \end{split}$$$

Therefore, one obtains

$$$\hat{\omega}(\vzero) = 0$$$

Substituting this into eq.\eqref{eq:p1} gives

$$$\hat{p}^{(1)}(\vzero, 1) = 0$$\label{eq:p101-Fourier}$

Secondly, from eq.\eqref{eq:p2} we know

$$$\begin{split} \hat{p}^{(2)}(\vzero, 1) &= \frac{1}{V}\sum_{\vk'} \hat{h}_3(\vzero, \vk', 1)\hat{\omega}(-\vk')\hat{\omega}(\vk') \\ &= \frac{1}{V}\sum_{\vk} \hat{h}_3(\vzero, \vk, 1)\hat{\omega}(-\vk)\hat{\omega}(\vk) \end{split}$$\label{eq:p201-h3}$

In the second line, we only change the summation variable from $\vk’$ to $\vk$. Noting that

$$$\begin{split} \hat{h}_3(\vzero, \vk', 1) &= \lim_{\vk\to 0} \hat{h}_3(\vk, \vk', 1) \\ &= \lim_{\vk\to 0} \frac{1}{k^2\abs{\vk-\vk'}^2}\left[ 1 - e^{-k^2} -\frac{k^2}{k^2-\abs{\vk-\vk'}^2} \left( e^{-\abs{\vk-\vk'}^2} - e^{-k^2} \right) \right] \\ &= \lim_{\vk\to 0} \frac{1}{k^2\abs{\vk-\vk'}^2}\left( 1 - e^{-k^2} \right) - \lim_{\vk\to 0} \frac{1}{k^2\abs{\vk-\vk'}^2} \frac{k^2}{k^2-\abs{\vk-\vk'}^2} \left( e^{-\abs{\vk-\vk'}^2} - e^{-k^2} \right) \\ &= \lim_{\vk\to 0} \frac{1}{k^2\abs{\vk-\vk'}^2}\left[ 1 - (1 - k^2) \right] + \frac{1}{k'^4}\left( e^{-k'^2} - 1 \right) \\ &= \frac{1}{k'^2} + \frac{1}{k'^4}\left( e^{-k'^2} - 1 \right) \end{split}$$$

Let

$$$x = k'^2$$$

Then it becomes

$$$\begin{split} \hat{h}_3(\vzero, \vk', 1) &= \frac{1}{x} + \frac{1}{x^2}\left( e^{-x} - 1 \right) \\ &= \frac{1}{x^2}\left( e^{-x} + x - 1 \right) \end{split}$$$

This is where the well known Debye function comes from, which is defined as

$$$\hat{g}_D(x) = \frac{2}{x^2}\left( e^{-x} + x - 1 \right)$$$

Therefore,

$$$\hat{h}_3(\vzero, \vk', 1) = \frac{1}{2}\hat{g}_D(x)$$\label{eq:h301-final}$

Insert eq.\eqref{eq:p101-Fourier}, \eqref{eq:p201-h3}, and \eqref{eq:h301-final} into eq.\eqref{eq:Q-expan-v2}, we arrive the final expansion of $Q$ in the Fourier space

$$$Q[w] \sim e^{-w_0}\left[ 1 + \frac{\epsilon^2}{2V^2}\sum_{\vk}\hat{g}_D(k^2)\hat{\omega}(-\vk)\hat{\omega}(\vk) \right]$$\label{eq:Q-expan-v3}$

or, by inverting the Fourier transforms,

$$$\begin{split} Q[w] &\sim e^{-w_0}\left[ 1 + \frac{\epsilon^2}{2V^2}\sum_{\vk}\hat{g}_D(k^2)\hat{\omega}(-\vk)\hat{\omega}(\vk) \right] \\ &\sim e^{-w_0}\left\lbrace 1 + \frac{\epsilon^2}{2V^2} \sum_{\vk}\hat{g}_D(k^2) \left[\int d\vx\;\omega(\vx)e^{i\vk\cdot\vx} \right] \left[\int d\vx'\;\omega(\vx')e^{-i\vk\cdot\vx'} \right] \right\rbrace \\ &\sim e^{-w_0}\left\lbrace 1 + \frac{\epsilon^2}{2V} \int d\vx\int d\vx'\;\left[\frac{1}{V}\sum_{\vk}\hat{g}_D(k^2) e^{i\vk\cdot(\vx-\vx')} \right]\omega(\vx)\omega(\vx') \right\rbrace \end{split}$$$

If we define the inverse transform of the Debye function as

$$$\begin{split} g_D(\abs{\vx - \vx'}) &= \frac{1}{V}\sum_{\vk}\hat{g}_D(k^2)e^{i\vk\cdot(\vx-\vx')} \\ &= \frac{1}{(2\pi)^3}\int d\vk\; \hat{g}_D(k^2)e^{i\vk\cdot(\vx-\vx')} \end{split}$$$

Then the expansion of $Q$ in real space can be written as

$$$Q[w] \sim e^{-w_0}\left[ 1 + \frac{\epsilon^2}{2V} \int d\vx \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx)\hat{\omega}(\vx') \right]$$\label{eq:Q-expan-final}$

### Expansion of the Density Operator

A weak inhomogeneity expansion for the segment density operator $\rho(\vx; [w])$ can be obtained in one of two equivalent ways. One way is expressing $\rho(\vx; [w])$ as an integral of the propagator, such as

$$$\rho(\vx; [w]) = \frac{1}{VQ[w]}\int_0^1 q(\vx, 1-s; [w])q(\vx, s; [w])$$$

and substituting the weak inhomogeneity expansion of the propagator eq.\eqref{eq:q-p} and \eqref{eq:p-expansion} into above equation and keeping to a certain order. The other way is performing a direct functional differentiation of eq.\eqref{eq:Q-expan-final} according to

$$$\begin{split} \rho(\vx; [w]) &= \rho(\vx; [W(\vx)]) \\ &= -\Der{\ln Q[W(\vx)]}{W(\vx)} \\ &= -N\Der{\ln Q[W(\vx)]}{NW(\vx)} \\ &= -N\Der{\ln Q[w]}{w} \end{split}$$$

Here we Follow the latter approach. From eq.\eqref{eq:Q-expan-final} we have

$$$\ln Q[w] \sim -w_0 + \ln\left[ 1 + \frac{\epsilon^2}{2V} \int d\vx \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx)\hat{\omega}(\vx') \right]$$$

From eq.\eqref{eq:w0-w} we have

$$$\Der{w_0}{w} = \frac{1}{V}$$$

From eq.\eqref{eq:w-omega} we have

$$$\Der{\omega}{w} = \frac{1}{\epsilon}$$$

Now we can perform the functional differentiation as

$$$\begin{split} \rho(\vx; [w]) &= -N\Der{\ln Q[w]}{w} \\ &= N\Der{w_0}{w} - \Der{}{w}\ln\left[ 1 + \frac{\epsilon^2}{2V} \int d\vx \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx)\hat{\omega}(\vx') \right] \\ &= \frac{N}{V} - N\Der{}{\omega}\ln\left[ 1 + \frac{\epsilon^2}{2V} \int d\vx \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx)\hat{\omega}(\vx') \right]\Der{\omega}{w} \\ &= \rho_0 - N\frac{1}{\epsilon}\Der{}{\omega}\ln\left[ 1 + \frac{\epsilon^2}{2V} \int d\vx \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx)\hat{\omega}(\vx') \right] \\ &= \rho_0 - N\frac{1}{\epsilon}\frac{1}{1 + \frac{\epsilon^2}{2V} \int d\vx \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx)\hat{\omega}(\vx')} \frac{\epsilon^2}{2V} 2\int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx') \\ &= \rho_0 - \rho_0 \epsilon \frac{1}{1 + \frac{\epsilon^2}{2V} \int d\vx \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx)\hat{\omega}(\vx')} \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx') \end{split}$$$

where $\rho_0 \equiv N/V$ is the volume-average segment density of a single chain. If we only retain the first order contribution of $\epsilon$, above equation can be simplified to

$$$\rho(\vx; [w]) = \rho_0 - \rho_0 \epsilon \int d\vx'\; g_D(\abs{\vx - \vx'})\omega(\vx')$$\label{eq:rho-omega}$

Or we can use eq.\eqref{eq:w-omega} to rewrite $\rho$ as a functional of $w$ instead of $\omega$

$$$\label{eq:rho-w} \rho(\vx; [w]) = \rho_0 \left[ 1 - \int d\vx'\; g_D(\abs{\vx - \vx'})(w-w_0) \right]$$$

We can also inverse the above equation to predict what external field should be applied if we want to obtain a certain segment density. Such inversion can be done in the Fourier space

$$$\hat{\phi}(\vk) = - \hat{g}_D(k^2)\hat{\Delta w}(\vk)$$\label{eq:phik-wk}$

where $\phi=\rho/\rho_0 - 1$ is the dimensionless fluctuation of segment density around its volume-average value and $\Delta w=w-w_0$ is the fluctuation of the external field. It is straightforward to inverse it

$$$\label{eq:wk-phik} -\hat{\Delta w}(\vk) = \hat{g}_D^{-1}(k^2) \hat{\phi}(\vk)$$$

## Appendix

### A. Derivation of eq.\eqref{eq:MDE-p1}

At $O(\epsilon^1)$, the expansion for $p(\vx, s)$ is

$$$p \sim p^{(0)} + \epsilon p^{(1)}$$$

Inserting it into eq.\eqref{eq:MDE-omega}, we have

$$$\der{}{s}\left[ p^{(0)} + \epsilon p^{(1)} \right] = \nabla^2 \left[ p^{(0)} + \epsilon p^{(1)} \right] - \epsilon\omega\left[ p^{(0)} + \epsilon p^{(1)} \right]$$\label{eq:MDE-full-A}$

From the zeroth order expansion eq.\eqref{eq:MDE-p0}, we know

$$$\der{p^{(0)}}{s} - \nabla^2 p^{(0)} = 0$$$

It simplifies eq.\eqref{eq:MDE-full-A} to

$$$\der{p^{(1)}}{s} = \nabla^2 p^{(1)} - \omega p^{(0)} - \epsilon\omega p^{(1)}$$$

As $\epsilon \to 0$, the last term in the right-hand side of the above equation vanishes, leading to eq.\eqref{eq:MDE-p1}.

### B. Derivation of eq.\eqref{eq:p1}

We know

$$$p^{(0)}(\vx, s) = 1$$$

which is the trivial solution of zeroth order expansion. Therefore eq.\eqref{eq:MDE-p1} becomes

$$$\der{}{s}p^{(1)}(\vx, s) = \nabla^2 p^{(1)}(\vx, s) - \omega(\vx)$$$

Perform Fourier transforms on both sides of above equation, we reaches

$$$\der{}{s}\hat{p}^{(1)}(\vk, s) = -k^2 \hat{p}^{(1)}(\vk, s) - \hat{\omega}(\vk)$$$

Re-organize above equation into the form

$$$\frac{d\left[ k^2\hat{p}^{(1)} + \hat{\omega} \right]}{k^2\hat{p}^{(1)} + \hat{\omega}} = -k^2 ds$$$

It can be easily solved and the solution is

$$$\ln\left[ k^2\hat{p}^{(1)} + \hat{\omega} \right] = -k^2 s + C$$\label{eq:p1-intermediate-A}$

Let $s=0$, and with $\hat{p}^{(1)}(\vk, 0)=0$ (because $p^{(1)}(\vx, 0)=0$), we then have

$$$\ln \hat{\omega} = C$$$

or

$$$\hat{\omega} = e^C$$$

Write eq.\eqref{eq:p1-intermediate-A} with $\hat{p}^{(1)}$ in the left-hand side and substitute above equation into it, we arrives at

$$$\hat{p}^{(1)} = -\frac{1}{k^2}\left( 1 - e^{-k^2s} \right)\hat{\omega}$$$

which is equivalent to eq.\eqref{eq:p1}.

### C. Derivation of eq.\eqref{eq:MDE-p2}

At $O(\epsilon^2)$, the expansion for $p(\vx, s)$ is

$$$p \sim p^{(0)} + \epsilon p^{(1)} + \epsilon^2 p^{(2)}$$$

Inserting it into eq.\eqref{eq:MDE-omega}, we have

$$$\der{}{s}\left[ p^{(0)} + \epsilon p^{(1)} + \epsilon^2 p^{(2)} \right] = \nabla^2 \left[ p^{(0)} + \epsilon p^{(1)} + \epsilon^2 p^{(2)} \right] - \epsilon\omega\left[ p^{(0)} + \epsilon p^{(1)} + \epsilon^2 p^{(2)} \right]$$\label{eq:MDE-full-p2-A}$

Substitute eq.\eqref{eq:MDE-p0} and \eqref{eq:MDE-p1} into above equation and ignore terms with $\epsilon$, it simplifies to eq.\eqref{eq:MDE-p2}.

### D. Derivation of eq.\eqref{eq:p2}

We copy eq.\eqref{eq:MDE-p2} here

$$$\der{p^{(2)}(\vx, s)}{s} = \lap p^{(2)}(\vx, s) - \omega(\vx)p^{(1)}(\vx, s)$$$

Perform Fourier transforms on both sides of above equation, we arrives

$$$\der{\hat{p}^{(2)}(\vk, s)}{s} = -k^2\hat{p}^{(2)}(\vk, s) - \widehat{\omega p^{(1)}}(\vk, s)$$$

This is a general first order linear partial differential equation with respect to $s$

$$$\der{\hat{p}^{(2)}}{s} + k^2\hat{p}^{(2)} = -\widehat{\omega p^{(1)}}$$$

The solution is

$$$\label{eq:p2-raw-A} \hat{p}^{(2)} = \frac{\int ds\; u(s)[-\widehat{\omega p^{(1)}}] + C}{u(s)}$$$

where $C$ is a constant which can be obtained by applying the initial condition, and

$$$u(s) = \exp\left( \int ds\;k^2 \right) = e^{k^2s}$$\label{eq:us-A}$

From the first order solution, we can find

$$$\begin{split} \widehat{\omega p^{(1)}} &= \hat{w} * \hat{p}^{(1)} \\ &= \frac{1}{V}\sum_{\vk'} \hat{w}(\vk') \hat{p}^{(1)}(\vk-\vk') \end{split}$$$

Substituting eq.\eqref{eq:p1} into above equation, we have

$$$\widehat{\omega p^{(1)}} = -\frac{1}{V}\sum_{\vk'} \hat{h}_2(\vk-\vk',s)\hat{w}(\vk-\vk')\hat{w}(\vk')$$\label{eq:wp1-A}$

To simplify the notation, we define

$$$\alpha \equiv k^2$$$ $$$\beta \equiv \abs{\vk-\vk'}^2$$$

Substitute eq.\eqref{eq:us-A} and \eqref{eq:wp1-A} into eq.\eqref{eq:p2-raw-A}

$$$\begin{split} \hat{p}^{(2)} &= e^{-\alpha s} \left\lbrace \int ds\; e^{\alpha s} \left[ \frac{1}{V}\sum_{\vk'}\frac{1}{\beta}(1-e^{-\beta s})\hat{w}(\vk-\vk')\hat{w}(\vk') \right] + C \right\rbrace \\ &= e^{-\alpha s} \left\lbrace \frac{1}{V\beta}\sum_{\vk'} \left[\int ds\; e^{\alpha s}(1-e^{-\beta s}) \right]\hat{w}(\vk-\vk')\hat{w}(\vk') + C \right\rbrace \\ &= e^{-\alpha s} \left\lbrace \frac{1}{V\beta} \sum_{\vk'} \left[ \frac{1}{\alpha}e^{\alpha s} - \frac{1}{\alpha-\beta}e^{(\alpha-\beta)s} \right]\hat{w}(\vk-\vk')\hat{w}(\vk') + C \right\rbrace \end{split}$$$

Because $p^{(2)}(\vx, 0)=0$, thus $\hat{p}^{(2)}(\vk, 0)=0$, therefore

$$$0 = e^{-\alpha 0} \left\lbrace \frac{1}{V\beta} \sum_{\vk'} \left[ \frac{1}{\alpha}e^{\alpha 0} - \frac{1}{\alpha-\beta}e^{(\alpha-\beta)0} \right]\hat{w}(\vk-\vk')\hat{w}(\vk') + C \right\rbrace$$$

Thus

$$$C = -\frac{1}{V\beta} \sum_{\vk'} \left( \frac{1}{\alpha} - \frac{1}{\alpha-\beta} \right)\hat{w}(\vk-\vk')\hat{w}(\vk')$$$

Substitute it back into previous equation,

$$$\begin{split} \hat{p}^{(2)} &= e^{-\alpha s} \left\lbrace \frac{1}{V\beta} \sum_{\vk'} \left[ \frac{1}{\alpha}e^{\alpha s} - \frac{1}{\alpha-\beta}e^{(\alpha-\beta)s} \right]\hat{w}(\vk-\vk')\hat{w}(\vk') -\frac{1}{V\beta} \sum_{\vk'} \left( \frac{1}{\alpha} - \frac{1}{\alpha-\beta} \right)\hat{w}(\vk-\vk')\hat{w}(\vk') \right\rbrace \\ &= \frac{1}{V} \sum_{\vk'} \left[ e^{-\alpha s}\left( \frac{1}{\alpha}e^{\alpha s} - \frac{1}{\alpha} - \frac{1}{\alpha-\beta}e^{(\alpha-\beta)s} + \frac{1}{\alpha-\beta} \right)\frac{1}{\beta} \right] \hat{w}(\vk-\vk')\hat{w}(\vk') \end{split}$$$

We can define

$$$\begin{split} \hat{h}_3(\vk, \vk', s) &= e^{-\alpha s}\left( \frac{1}{\alpha}e^{\alpha s} - \frac{1}{\alpha} - \frac{1}{\alpha-\beta}e^{(\alpha-\beta)s} + \frac{1}{\alpha-\beta} \right)\frac{1}{\beta} \\ &= \frac{1}{\beta}\left( \frac{1}{\alpha} - \frac{1}{\alpha}e^{-\alpha s} - \frac{1}{\alpha-\beta}e^{-\beta s} + \frac{1}{\alpha-\beta}e^{-\alpha s} \right) \\ &= \frac{1}{\alpha\beta}\left[ 1 - e^{-\alpha s} - \frac{\alpha}{\alpha-\beta}\left( e^{-\beta s} - e^{-\alpha s} \right) \right] \end{split}$$$

which completes the derivation.