Is there a natural way to define relative wavefunctions? Suppose I have two free-particles $A,B$ whose positions evolve according to wave functions $\psi_A(x,t)$ and $\psi_B(x,t)$. I am interested in talking about the position vector from $A, B$, and ideally I would like to describe it according to a wave function lets call it $\psi_{B,A}(x,t)$. 
In attempting to find this wave function I did the following: 


*

*I note that $|\psi_{B,A}(x,t)|^2$ strictly describes the probability of the position vector from $A$ to $B$. 

*I note that $$|\psi_{B,A}(x,t)|^2 = \int_{-\infty}^{\infty} |\psi_B(x-u, t)|^2 |\psi_A(u, t) |^2 \partial u $$ from a heuristic argument that for a particular position vector $x$ the probability density of this position vector should be equal to a weighted sum of all the densities of all pairs of vectors $u_1, u_2$ distributed according to B and A respectively, whose difference is $x$. 
But now I want to get an expression for $\psi_{B,A}(x,t)$ not the square of its absolute value, and this is where it gets hard. I tried looking at the equation of (2) and breaking the absolute values into pairs of value and conjugate. Then I hoped to differentiate the expressions according time, and use schrodinger equation. But it seems impossible to me that I will ever be able to generate a linear PDE from this, let alone that I will generate a linear PDE that is purely a function of $x,t, \psi_{A,B}$. 
One Idea:
I'm still toying with this one, but one could make an equation stating that the relative momentum from $A$ to $B$ is distributed according to $\phi_{B,A}$ which we define to be the fourier transform of $\psi_{B,A}$ and moreoever
$$ |\phi_{B,A} (p,t)|^2 = \int_{-\infty}^{\infty} |\phi_{B} (p -s, t)|^2 |\phi_A(s, t)|^2  \partial s$$ 
 and see what that yields.
 A: Why not just redefine your variables? If $\vec{x}_{A}$ is the position of particle $A$ and $\vec{x}_{B}$ that of particle $B$, then
$$\vec{x}_{A}=\vec{X}+\frac{m_{B}}{M}\vec{x} \qquad \text{and}\qquad\vec{x}_{B}=\vec{X}-\frac{m_{A}}{M}\vec{x}$$
where $m_{A}$ and $m_{B}$ are the masses of the particles, $M=m_{A}+m_{B}$ is the total mass of the system, $\vec{x}=\vec{x}_{A}-\vec{x}_{B}$ is the separation vector and
$$
\vec{X}=\frac{m_{A}\vec{x}_{A}+m_{B}\vec{x}_{B}}{M}
$$
is the position of the center of mass.
If $\psi_{A}(\vec{x}_{A})$ describes particle $A$ and $\psi_{B}(\vec{x}_{B})$ describes particle $B$, the composite system is described by:
$$
\psi(\vec{x}_{A},\vec{x}_{B})=\psi_{A}(\vec{x}_{A})\,\psi_{B}(\vec{x}_{B})\quad\text{(distinguishable particles)}\\
\psi(\vec{x}_{A},\vec{x}_{B})=\frac{1}{\sqrt{2}}\ (\psi_{A}(\vec{x}_{A})\,\psi_{B}(\vec{x}_{B})+\psi_{A}(\vec{x}_{B})\,\psi_{B}(\vec{x}_{A}))\quad\text{(indistinguishable particle, bosons)}\\
\psi(\vec{x}_{A},\vec{x}_{B})=\frac{1}{\sqrt{2}}\ (\psi_{A}(\vec{x}_{A})\,\psi_{B}(\vec{x}_{B})-\psi_{A}(\vec{x}_{B})\,\psi_{B}(\vec{x}_{A}))\quad\text{(indistinguishable particle, fermions)}
$$
(one should also interchange the spin labels for the indistinguishable cases). Let us assume for simplicity that the particles are distinguishable. Then you can define a wavefunction depending on the center of mass position and separation vector as
$$
\psi(\vec{X},\vec{x})=\psi_{A}\big(\vec{X}+\frac{m_{B}}{M}\vec{x}\big)\,\psi_{B}\big(\vec{X}-\frac{m_{A}}{M}\vec{x}\big)
$$
and this answers your question. Notice that, in general, the wavefunction cannot be factorized into a $\vec{X}$-dependent and a $\vec{x}$-dependent part: this depends on the specific form of $\psi_{A}$ and $\psi_{B}$. However, you still have a well-defined probability for the variable $\vec{x}$:
$$
\rho(\vec{x})=\int d^{3}X\ |\psi(\vec{X},\vec{x})|^{2}=\int d^{3}X\ \bigg|\psi_{A}\big(\vec{X}+\frac{m_{B}}{M}\vec{x}\big)\,\psi_{B}\big(\vec{X}-\frac{m_{A}}{M}\vec{x}\big)\bigg|^{2}
$$
Notice that this definition - at least for the case at hand - leads to the same result as your own's: define $\vec{U}=\vec{X}+\frac{m_{B}}{M}\vec{x}$ and change variables of integration from $\vec{X}$ to $\vec{U}$: you get
$$
\rho(\vec{x})=\int d^{3}U\ \bigg|\psi_{A}\big(\vec{U}\big)\,\psi_{B}\big(\vec{U}-\frac{m_{B}}{M}\vec{x}-\frac{m_{A}}{M}\vec{x}\big)\bigg|^{2}=\int d^{3}U\ \bigg|\psi_{A}\big(\vec{U}\big)\,\psi_{B}\big(\vec{U}-\vec{x}\big)\bigg|^{2}
$$
where in the second step I only simplified $(m_{A}+m_{B})/M=1$. However the result only holds if the particles are distinguishable: if they are not, then $\psi(\vec{X},\vec{x})$ takes on one of the expressions that I wrote above for the bosonic or fermionic case, and my definition will lead to a different result. The one I wrote is the correct one because it is based on sound principles (composite wavefunctions for in/distinguishable particles) rather than on a heuristic.
On the other hand, it can be shown that if the potential acting on the system $A+B$ is the sum of a $\vec{X}$-dependent potential and a $\vec{x}$-dependent potential, then the eigenfunctions of the Hamiltonian indeed always factorize into a $\vec{X}$-dependent and a $\vec{x}$-dependent part: in this case you can always write the eigenfunctions as
$$
\psi(\vec{X},\vec{x})=\psi_{c.o.m.}(\vec{X})\ \psi_{sep.}(\vec{x})
$$
where each of the $\psi$'s on the R.H.S. separately solves the time-independent Schrödinger equation with the appropriate potential. In this case, $\psi_{sep.}(\vec{x})$ is the wavefunction you are looking for (and it may have specific simmetry properties in the indistinguishable case). For example, if the particles are acted upon by a force that only depends on the interparticle separation $\vec{x}$, you will find that
$$
\psi_{c.o.m.}(\vec{X})=e^{i\vec{P}\cdot \vec{X}/\hbar}
$$
where $\vec{P}$ is the center of mass momentum of the system while $\psi_{sep.}(\vec{x})$ satisfies the Schrödinger equation
$$
\bigg[-\frac{\hbar^{2}}{2\mu}\ \nabla^{2}+V(\vec{x})\bigg]\psi_{sep.}(\vec{x})=E_{int}\ \psi_{sep.}(\vec{x})
$$
with $\mu=m_{A}m_{B}/M$ the reduced mass of the system.
