Derivation in LSZ Reduction Formula In deriving the LSZ formula, a crucial step is to show
$$\langle|a_{p}^{\dagger}|\rangle=-i\int dx^0  \int \mathrm{d}^{3} x \partial_{0}\langle | e^{-i p\cdot x} \overleftrightarrow{\partial_{0}} \phi(x)|\rangle=-i\int d^4x e^{-ip\cdot x}(\square+m^2)\langle|\phi|\rangle$$
where for conciseness, I omitted the arguments in $\langle|,|\rangle$.
But, from our experience with QFT, the creation operator can equally well be represented by the field operator $\phi(x)$ and its canonical conjugate $\Pi(x)$.
$$a_{p}^{\dagger}=-i \int d^{3} x e^{-i p \cdot x}\left[\Pi(x)+i E_{p} \phi(x)\right].$$
Then, using Hamilton's equations for fields
$$\begin{aligned}
&\partial_{0} \phi(x)=\frac{\delta \mathcal{H}}{\delta \Pi(x)}=\Pi(x), \\
&\partial_{0} \Pi(x)=-\frac{\delta \mathcal{H}}{\delta \phi(x)}=\left(\nabla^{2}-m^{2}\right) \phi(x)
\end{aligned}$$
we would have
$$\begin{align}
\partial_{0}(e^{-ip\cdot x}\left[\Pi+i E_{p} \phi\right]) =&- ie^{-ip\cdot x}E_p \left[\Pi+i E_{p} \phi\right]+e^{-i p\cdot x}\left[(\nabla^2-m^2)\phi+iE_p \Pi\right]\\
=&e^{-ip\cdot x}(E_p^2+\nabla^2-m^2)\phi=e^{-ip\cdot x}(-\partial^2/\partial t^2+\nabla^2-m^2)\phi
\\=&-e^{-ip\cdot x}(\square+m^2)\phi
\end{align}.$$
Notice the extra minus sign!
On the other hand
$$\partial_{0} \left[e^{-i p\cdot x} \overleftrightarrow{\partial_{0}} \phi(x)\right]=e^{-ip\cdot x}(-E_p^2-\nabla^2+m^2)\phi=e^{-ip\cdot x}(\square+m^2)\phi$$
as desired. This is the correct result to be used.

As is apparent, I'm rather careless about these calculations, so hopefully nothing goes seriously messed. Assuming that is the case, these two expressions are clearly different by a minus sign! What is the origin for this discrepancy and how should we interpret this result?
This question emerges as I work my way through the QFT book by Gelis (section 1.5 specifically), the book might contain more relevant details, in case one finds my re-expression of materials makes no sense at all.
 A: You have just shown that
$$\partial_0\left[e^{-ipx}\overleftrightarrow{\partial_0}\phi(x)\right]=e^{-ipx}\left(\Box + m^2\right)\phi(x) \tag{1}$$
Now you just need to show that
$$a^{\dagger}_p = -i \int dx_0 \int d^3x\, \partial_0 \left[ e^{-ip\cdot x} \,\overleftrightarrow{\partial_0} \phi(x)\right] \tag{2}$$
To do this, just use the expression for $a^{\dagger}_p$ in terms of $\phi(x)$ and $\Pi (x)$ that you gave. Then use the Hamilton equations that you quoted, plugging in for $\Pi (x)$.

Edit:
I believe what you are actually trying to calculate is
$$a^{\dagger}_p(\infty)-a^{\dagger}_p(-\infty) = -i \int dx_0 \int d^3x\, \partial_0 \left[ e^{-ip\cdot x} \,\overleftrightarrow{\partial_0} \phi(x)\right] \tag{3}$$
Using that $p\cdot x = E x^0 - \vec x \cdot \vec p\,\,$, so $\,\,\partial_0 e^{-ip\cdot x} = -iE e^{-ip\cdot x}$,
\begin{align}
a_{p}^{\dagger}(t)&=-i \int d^{3} x \,e^{-i p \cdot x}\left[\Pi(x)+i E_{p} \phi(x)\right] \\
&= -i \int d^{3} x \, e^{-i p \cdot x}\left[\partial_0 \phi(x)+i E_{p} \phi(x)\right] \\
&= -i \int d^{3} x \,\left[e^{-ip\cdot x}\partial_0 \phi(x)-\left(\partial_0 e^{-ip\cdot x}\right)\phi(x)\right] \\
&= -i \int d^{3} x \,\left[e^{-ip\cdot x}\overleftrightarrow{\partial_0} \phi(x)\right] \tag{4} 
\end{align}
This therefore implies the eq (3), by integrating a total derivative.
$$a^{\dagger}_p(\infty)-a^{\dagger}_p(-\infty) = -i \int_{-\infty}^{\infty} dx^0\,\partial_0 \int d^{3} x \,\left[e^{-ip\cdot x}\overleftrightarrow{\partial_0} \phi(x)\right]$$
Similarly if we take the adjoint of the above equation, and use that $\phi(x)$ is Hermitian, we get
$$a_p(\infty)-a_p(-\infty) = i \int_{-\infty}^{\infty} dx^0\,\partial_0 \int d^{3} x \,\left[e^{ip\cdot x}\overleftrightarrow{\partial_0} \phi(x)\right]$$
A: I might have sorted this question myself. The Hamilton's equations are derived under the assumption that $\phi$ is a free scalar field, which is solved by the Klein-Gordon equation, so the minus sign does not matter at all. For fields with interaction, we should probably avoid using the Hamilton's equations directly, instead the exponential term does satisfy the Klein-Gordon equation, so there will be a term
$$\phi\frac{\partial^2}{\partial {x_0}^2} e^{ip\cdot x}=\phi(\nabla^2-m^2)e^{ip\cdot x}$$
Integrate by parts, and notice that combined with another term
$$e^{ip\cdot x}\frac{\partial^2}{\partial {x_0}^2}\phi$$
we would get
$$e^{ip\cdot x}(\square+m^2)\phi$$
without any ambiguity.
