Wightman Function for complex scalar field - timelike separations? For a complex scalar field $\Phi$, the field has the expansion
$$
\Phi(x^0,\mathbf{x}) = \int \frac{d^{3}\mathbf{p}}{\sqrt{ 2 E_{\mathbf{p}} (2\pi)^3 } }\ \bigg[ e^{- i E_{\mathbf{p}}x^0 + i \mathbf{p} \cdot \mathbf{x} } a_{\mathbf{p}} +  e^{+ i E_{\mathbf{p}}x^0 - i \mathbf{p} \cdot \mathbf{x} } \bar{a}_{\mathbf{p}}^{\ast}  \bigg]
$$
where $E_{\mathbf{p}} = \sqrt{ |\mathbf{p}|^2 + m^2 }$ and $ a_{\mathbf{p}}$, $a_{\mathbf{p}}^{\ast}$ are the particle ladder operators, and $\bar{a}_{\mathbf{p}}$, $\bar{a}_{\mathbf{p}}^{\ast}$ are the antiparticle ladder operators.
With the canonical commutation relations (meaning $[a_{\mathbf{p}},a_{\mathbf{k}}^{\ast}] = [\bar{a}_{\mathbf{p}},\bar{a}_{\mathbf{k}}^{\ast}] = \delta^{(3)}(\mathbf{k} - \mathbf{p})$, etc.), the Wightman function becomes the following:
$$
\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \int \frac{d^{3}\mathbf{p}}{(2\pi)^3} \frac{1}{2E_{\mathbf{p}}} e^{- i E_{\mathbf{p}} (x^0-y^0) + i \mathbf{p} \cdot ( \mathbf{x} - \mathbf{y} ) }
$$
In Weinberg's Volume I (Chapter 5.2), he evaluates this function by picking the separation to be space-like with $(x-y)^2 = - (x^0-y^0)^2 + |\mathbf{x} - \mathbf{y}|^2 > 0$. This results in
$$
\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \frac{m}{4\pi^2 \sqrt{ (x-y)^2 } } K_{1}\big( m \sqrt{ (x-y)^2 } \big)
$$
where $K_1$ is the $1^{\mathrm{st}}$ modified Bessel function of the second kind. 
I have two questions:
1. How to compute this function for time-like separations? If $(x-y)^2 <0$ is it simply true that $\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \frac{m}{4\pi^2 \sqrt{ - (x-y)^2 } } K_{1}\big( m \sqrt{ - (x-y)^2 } \big)$?
2. Weinberg says that for spacelike separations $(x-y)^2 > 0$, this function is symmetric under $(x-y)\mapsto (y-x)$. Meaning $\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \langle 0 | \Phi^{\ast}(y) \Phi(x) |  0  \rangle$. I see that this is obvious from the formula involving the Bessel $K_1$, but how does one see this from the integral representation $\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \int \frac{d^{3}\mathbf{p}}{(2\pi)^3} \frac{1}{2E_{\mathbf{p}}} e^{- i E_{\mathbf{p}} (x^0-y^0) + i \mathbf{p} \cdot ( \mathbf{x} - \mathbf{y} ) }$? To me this is not obvious from the integral representation for non-zero $x^0-y^0$...Is it still true that the Wightman function is symmetric under $(x-y)\mapsto (y-x)$ for timelike separations?
EDIT; a typo. I actually am interested in the  Wightman function $\langle 0 | \Phi^{\ast}(x) \Phi(y) |0 \rangle$ (note that $\langle 0 | \Phi(x) \Phi(y) |0 \rangle = 0$)
 A: The fundamental ingredient to properly understand the situation in your second question is the transformation law of the scalar field under general Lorentz transformations. So, let $\Lambda$ be a generic Lorentz transformation such that it rotates a coordinate system $x$ into $x'= \Lambda x$. Let $U(\Lambda)$ be a (anti)unitary representation of $\Lambda$, then we have the following transformation law for a sclar field:
$$U^{\dagger}(\Lambda) \Phi(x')U(\Lambda)= \Phi(x)$$
In the case of parity tranformations $U(\Lambda) \equiv P$, the field may acquire an intrinsic parity factor $P \Phi(x_0,\mathbf{x}) P^{\dagger} = \eta \Phi(x_0,-\mathbf{x})$ in front of it with $\eta=\pm1$, if the field is scalar is $+1$, if it pseudo-scalar, like a Pion it is $-1$ .
Time reversal symmetry is symilar to parity as $T \Phi(x_0,\mathbf{x}) T^{\dagger} = \xi \Phi(-x_0,\mathbf{x})$, where $\xi=\pm1$, fortunately for the causal structure of the Universe, $T$ is anti-unitary $\langle \psi_1,T \psi_2\rangle = \langle T^{\dagger}\psi_1,\psi_2\rangle^{*}$ which is something that gets lost in Dirac notation.
Now come back to the  definition  of the Wightman function 
$$
W(x-y)=\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle \tag 1
$$
If you insert the combination of unitary operators $U^{\dagger}U=1$ into the right hand side of (1), you will get
$$
\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \langle 0 |U^{\dagger}U \Phi^{\ast}(x)U^{\dagger}U \Phi(y) U^{\dagger}U|  0  \rangle = \langle 0 | \Phi^{\ast}(x') \Phi(y') |  0  \rangle \tag 2
$$
since the vacuum state is invariant under Lorentz transformations, i.e. $U | 0 \rangle = | 0 \rangle$. Therefore
$$W(x-y)=W(x'-y')$$
This is important because it means that we can do a Lorentz transformation and choose some new coordinates $x'$ and $y'$  without changing the Wightman function. All results which are true in a given frame are  true in any other frame according to Eq. 2. Incidentally, this is also true for parity :
$$
\langle 0 | \Phi^{\ast}(x_0,\mathbf{x}) \Phi(y_0,\mathbf{y}) |  0  \rangle = \langle 0 | \Phi^{\ast}(x_0,-\mathbf{x}) \Phi(y_0,-\mathbf{y}) |  0  \rangle\tag 3
$$
now, remember that a Wightman function is a function $W=W(x-y)$ only of the difference of the coordinates.
Therefore a parity transformation is exactly the same as exchanging $\mathbf{x}$ with $\mathbf{y}$ in the integral:
$$\langle 0 | \Phi^{\ast}(x_0,-\mathbf{x}) \Phi(y_0,-\mathbf{y}) |  0  \rangle =\int \frac{d^{3}\mathbf{p}}{(2\pi)^3} \frac{1}{2E_{\mathbf{p}}} e^{ -i E_{\mathbf{p}} (x^0-y^0) -i \mathbf{p} \cdot ( \mathbf{x} - \mathbf{y} ) } =  \langle 0 | \Phi^{\ast}(x_0,\mathbf{y}) \Phi(y_0,\mathbf{x}) |  0  \rangle  $$
 and so
$$\langle 0 | \Phi^{\ast}(x_0,\mathbf{x}) \Phi(y_0,\mathbf{y}) |  0  \rangle =  \langle 0 | \Phi^{\ast}(x_0,\mathbf{y}) \Phi(y_0,\mathbf{x}) |  0  \rangle  $$
which means that regardless of the spacetime separation of $x$ and $y$ you have
$$W(x_0-y_0,\mathbf{x} - \mathbf{y})=W(x_0-y_0,\mathbf{y} - \mathbf{x})$$
if parity is a good symmetry for the theory.
Time reversal, being anti-unitary, would not simply give the exchange between $x_0$ and $y_0$ but it also would conjugate the whole Wightman function (which is easy to check flips the coordinates but also their sign). Therefore, you see that in general there is no four dimensional symmetry for the exchange between $x$ and $y$.
The special case of spacelike separations, in which the coordinate exchange symmetry is four dimensional, can be more easily understood (I hope) in the context of parity because it is always possible to find an inertial frame of reference such as $x_0'=y_0'$ and therefore spatial symmetry coincides with the four dimensional one. More rigorously you can say that only for space like separations there exists a continuous Lorentz transformation that brings $(x-y)\rightarrow -(x-y)$. This statement can be proven by extending the Lorentz group into the complex plane (complex angles, etc). This should also hint at the fact that the solution to your first question can't be the Bessel function of the second kind. 
If you want to prove this spatial symmetry directly using the integral representation, the first step should be to make the following coordinate transformation $\mathbf{p}\rightarrow -\mathbf{p}$ to obtain $\langle 0 | \Phi^{\ast}(x_0,\mathbf{x}) \Phi(y_0,\mathbf{y}) |  0  \rangle $ as it should be. But $\mathbf{p}\rightarrow -\mathbf{p}$ is just a parity transformation on momenutm, so we need to prove that the measure of the integral is Lorentz invariant so that it won't change under a parity transformation.
To do this we use the following identity
$$\int \frac{d^{3}\mathbf{p}}{(2\pi)^3} \frac{1}{2E_{\mathbf{p}}} f(p)=\int d^{4}p \delta(p^2-m^2) \theta(p^0) f(p)$$
with $p^2=-p_\mu p^\mu$. To prove the identity it suffices to do the following
$$\delta(p^2-m^2)=\delta(p_0^2-\mathbf{p}^2-m^2) \rightarrow p^0=\pm E_p$$
therefore
$$\delta(p^2-m^2)=\frac{\delta(p_0-E_p)}{2 E_p}+\frac{\delta(p_0+E_p)}{2 E_p}$$
the $\theta(p_0)$ removes the negative $p_0$ and there you have it.
The Proof of Lorentz invariance is rather easy now, $d^4p = |det(\Lambda)|d^4p'=d^4p'$, since $|det(\Lambda)|= 1$. The Dirac delta being a function of Lorentz scalars is invariant. The Theta is invariant since $p_\mu p^\mu = -m^2$ is a timelike vector and so the sign of $p^0$ is conserved under Lorentz transformations.
So, the measure is invariant, we can change $\mathbf{p}\rightarrow -\mathbf{p}$ and prove the result from the integral point of view. 
Concerning your first question about the explicit shape of the soultion. Consider that the Wightman function is a Lorentz invariant, so it must be a function only of (x-y)^2. Consider the space-like case $x_0'=y_0'$, the integral on the right hand side of Eq. 1 becomes
$$
\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \int \frac{d^{3}{p}}{(2\pi)^3} \frac{1}{2E_{{p}}} e^{ i {p} \cdot ( \mathbf{x} - \mathbf{y} ) } = \frac{4 \pi}{(2\pi)^3}\int \frac{{p}^2d{p}}{(2\pi)^3} \frac{1}{2E_{{p}}} \frac{\sin(p |\mathbf{x} - \mathbf{y} |)}{p  \mathbf{x} - \mathbf{y} }
$$
whit $|\mathbf{x} - \mathbf{y}| = \sqrt{({x} - {y})^2 }$ which can be proven to be 
$$
\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \frac{m}{4\pi^2 \sqrt{ (x-y)^2 } } K_{1}\big( m \sqrt{ (x-y)^2 } \big) 
$$
as you did write.
Now consider the timelike case and let's choose a reference frame where $\mathbf{x}=\mathbf{y}$:
$$
\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \int \frac{d^{3}{p}}{(2\pi)^3} \frac{1}{2E_{{p}}} e^{ -i E_p  ( x_0 - y_0 ) } 
$$
the measure $d^3p= d\Omega p^2 dp$ where $d\Omega$ are the angular variables. Now let's do the following trick
$$E = \sqrt{p^2+m^2}\rightarrow d E = \frac{p}{\sqrt{p^2+m^2}} dp\rightarrow E dE= pdp$$
$$
\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \frac{2 \pi}{(2\pi)^3}\int dE\ \sqrt{E^2-m^2}\ e^{ -i E ( x_0 - y_0 ) } 
$$
with $x_0-y_0=\sqrt{-(x-y)^2}$ (assuming $x_0>y_0$). This is can be found in a table to be the Hankel function 
$$\langle 0 | \Phi^{\ast}(x) \Phi(y) |  0  \rangle = \frac{i m}{8 \pi \sqrt{-(x-y)^2}} H^{(2)}_{1}(m\sqrt{-(x-y)^2})$$
For timelike separations it is a completely different function from the one found in the spacelike case.
A tedious albeit complete demonstration of these integral formulas can be found in "W. Greiner and J. Reinhardt, Quantum Electrodynamics, Springer (1992)".
