Let me first answer your questions:
Q1: No, the partial transposition acts on one part of the subsystem. To make this more clear, let me give a real definition of the criterion:
Let $\mathcal{H}_1,\mathcal{H}_2$ be two Hilbert spaces, then the partial transpose on the second system is defined via
$$ ^{PT}: \mathcal{B}(\mathcal{H}_1\otimes \mathcal{H}_2)\to \mathcal{B}(\mathcal{H}_1\otimes \mathcal{H}_2);\quad A\otimes B\mapsto A\otimes B^T $$
Obviously, this only defines the partial trace on product states, but the rest is done by linearity. Now a qubit (!) state is entangled if and only if the partial trace on any subsystem is not positive. However, this is equivalent to saying that a state is entangled if and only if the partial trace on the second subsystem is nonnegative, because the partial trace on the first subsystem would map $A\otimes B\mapsto A^T\otimes B=(A\otimes B^T)^T$, i.e. if you compose the partial trace on the second subsystem with the transposition on the whole system, you obtain the partial trace on the first subsystem. Since the transposition leaves the eigenvalues invariant, if a state is not positive under the partial transposition of any one subsystem, it will remain nonpositive under the partial transposition of the other system.
Let me emphasize however (as you stated) that this criterion is not sufficient for proving separability on any larger system than a system comprised of qubits and qutrits!
Q2: Purity and entanglement are completely different concepts in the sense that purity is about a single system, whereas entanglement is about a bipartition of the system, so I wouldn't really expect a nice rank characterization. However, if you have a pure state $\rho=|\psi\rangle\langle \psi|$, then you can actually say something, because such a pure state is entangled if and only if its reduced density matrix is pure - and since this has a characterization with matrix ranks, you have such a connection. This connection is however false for mixed entangled/separable states and I don't know whether such a connection exists.
Let me give a short proof of the claim:
Let $|\psi\rangle =\sum_{i=1}^k \lambda_i |\psi^1_i\rangle|\psi^2_i\rangle$ be an arbitrary mixed state in its Schmidt-decomposition (note that the state is separable iff $k=1$. Now the reduced density matrix of the first subsystem is easily seen to be:
$$ \rho_{red}= \sum_{i=1}^k \lambda_i^2 |\psi^1_i\rangle\langle \psi^1_i|$$
which has as rank the Schmidt-rank of $\rho$. Since a state is pure iff its density matrix is rank one, the reduced density matrix is pure iff the state was separable.
Actually, this characterization is exactly Exercise 2.78 in Nielsen&Chuang!
Q3: This was not actually stated, but you say that you have come across this condition that $\rho$ is pure iff it is rank one. Let me clarify the issue a bit: By definition, $\rho$ is pure iff it is of the form $|\psi\rangle\langle \psi|$ for some $\psi$. Obviously, this matrix is rank one, since it has one eigenvector with eigenvalue 1 (namely $|\psi\rangle$ - granted the state is normalized) and all other eigenvectors are 0 (all other states of an orthonormal basis comprising $|\psi\rangle$). On the other hand, if a state is rank 1, this means that it has only one eigenvector $|\psi\rangle$ with eigenvalue 1 and all other eigenvalues are zero. Since the density matrix is positive, via the spectral decomposition this implies that $\rho=|\psi\rangle\langle \psi|$.
Your Problem: Let me first point out a few misconceptions and misnotations and problems with the way you approach your problem. Later, I will give a solution.
You write that $\tau:\mathcal{H}_1\to\mathcal{H}_1$ is a linear map. This means that given the states $|x_i\rangle$, the map $\tau$ is actually represented by a matrix, which I am going to denote by $T$. Of course, $TT^*$ is then a positive operator and you could say that it is a "state", but I don't think one can ascribe the same phyiscal meaning. Then you write:
$$\sum_{i,j} |\tau x_i \rangle \langle \tau x_j | |y_i\rangle\langle y_j| = \sum_{i,j} |\rho x_i \rangle \langle x_j| |y_i\rangle \langle y_j|$$
This line is actually wrong. You can't just pull the $\tau$ to the other side like you want to. The problem is that the matrix on the right hand side is most likely not even Hermitian, while the matrix on the right hand side most definitely is! Therefore, everything that follows is wrong.
Still, there is one other problem that I'd like to address: You want to prove separability by showing that the partial transpose is positive. Let me repeat again: This is not possible. If the partial transpose is not positive, the state is entangled, but the other direction (If the state is entangled, the partial transpose is not positive, or else: If the partial transpose is positive, the state is separable) does not hold in dimensions other than $2\times 3$ and $2\times 2$. You sometimes do talk about "qubits", which would imply that $\mathcal{H}_1=\mathcal{H}_2=\mathbb{C}^2$, but you did not really specify this.
The final problem, and this is what I was getting at before: You did not specify the range of the indices of your sum in $|\Psi\rangle$. If it is the dimension of $\mathcal{H}_1$, the statement you want to prove is essentially correct, if it is not, however, the statement you want to prove is totally wrong.
My solution:
So we are given a state $|\Psi_0\rangle=\sum_{i=1}^k |x_i\rangle|y_i\rangle$ (note that the state is in Schmidt decomposition!). We want to say something about the separability of the state $|\Psi\rangle=\sum_{i=1}^k |\tau x_i\rangle|y_i\rangle$.
By definition $\tau:\mathcal{H}_1\to \mathcal{H}_1$, i.e. if we suppose that $\{x_i\}_{i=1}^n$ is a basis, we can define $\tau$ via $\tau(x_i)=\sum_j \lambda_{ji} x_j$. So far, nothing new.
We now ask: When is $\Psi$ separable? First note that whatever $\tau$ is, $|\Psi\rangle$ is still pure, because the map is defined on the level of pure states, e.g. we will not always have to go to density matrices.
Now let's have a look at the problem for different numbers of elements in the sum, i.e. for different $k$ in $|\Psi_0\rangle$. First cover the easy cases:
k=1: We have $|\Psi_0\rangle=|x_1\rangle|y_1\rangle$, which means that we start out with a separable state. In this case,
$$|\Psi\rangle=|\tau x_1\rangle|y_1\rangle=\left(\sum_j \lambda_{1j}|x_j\rangle\right)\otimes |y_1\rangle $$
is still separable. This is completely independent of $\tau$.
k=n: This means we have essentially a maximally mixed state. To properly normalize it, we should write $|\Psi_0\rangle=\frac{1}{\sqrt{n}} \sum_{i=1}^n |x_i\rangle|y_i\rangle$. We now apply $\tau$ and we want the outcome to be separable, i.e. we want:
$$ |\Psi\rangle=\frac{1}{\sqrt{n}} \sum_{i=1}^n |\tau x_i\rangle|y_i\rangle=\frac{1}{\sqrt{n}} \sum_{i=1}^n \left(\sum_{j=1}^n \lambda_{ji}|x_j\rangle\right)|y_i\rangle \stackrel{!}{=} |\phi\rangle|\varphi\rangle $$
for some states $|\phi\rangle$ and $|\varphi\rangle$. This last bit is the criterion for separability. Now note that we have the sum and by multilinearity of the tensor product, the only way to get $|\phi\rangle$ and $|\varphi\rangle$ as above would be to have $|\varphi\rangle=\sum_{i=1}^n |y_i\rangle$, i.e. the first part must be independent of $i$. This means $\sum_{j=1}^n \lambda_{ji}|x_j\rangle$ must be equal to some state $|\phi\rangle$ independent of $i$, which means (since the $x_i$ form a basis) that $\tau$ must map $|x_i\rangle$ to $\lambda_i |\phi\rangle$ (there may be a factor depending on $i$ because $\lambda_i|\phi\rangle \otimes |\varphi\rangle = |\phi\rangle \otimes \lambda_i|\varphi\rangle$ always).
This means: $|\Psi\rangle$ is separable if and only if $\tau$ maps any state of the Hilbert space $\mathcal{H}_1$ to a state proportional to $|\phi\rangle$. The proportionality factor is dictated by linearity - in fact, most states must be mapped to zero. More to the point, the matrix representing $\tau$ must have rank one, because the image is one-dimensional. This means that the matrix $\tau\tau^*$ is also rank one if and only if $|\Psi\rangle$ is separable.
$1<k<n:$ This case lies in between the two extreme cases above. It is more cumbersome, because (if I am not mistaken) the result will be that $|\Psi\rangle$ is separable iff $\tau$ maps the vector space spanned by $\{|x_i\rangle\}_{i=1}^k$ to some arbitrary or fixed state $|\phi\rangle$, while the complement spanned by $\{|x_i\rangle\}_{k+1}^n$ (if we assume that all $x_i$ are orthogonal) is arbitrary.
Neat question! Here is an almost complete answer; I encourage someone to edit my answer to tighten up the math at the end, but all of the steps are clear.
Inspired by another inequality for fidelity, let's resolve the identity using states that are Fourier transforms of eigenstates $\hat{\mathbf{B}}$; i.e., using states that are displaced by $e^{-i x t \hat{\mathbf{B}}}$:
$$\mathbb{I}=\int dx |x\rangle\langle x|;\quad e^{-i y t \hat{\mathbf{B}}}|x\rangle=|x+ y t\rangle; \quad \langle x|y\rangle=\delta(x-y).$$ These can always be found by using $|x\rangle=\int dp e^{-ixp}|p\rangle$ with $|p\rangle$ the delta-orthogonal eigenstates of $\hat{\mathbf{B}}$.
We use that the trace norm can always be rewritten using $\mathrm{Tr}(|A|)=\mathrm{Tr}(A U)$ for some unitary $U$ and the usual $|A|=\sqrt{A A^\dagger}$. The fidelity $F(\rho,\sigma;t)=||\sqrt{\rho(t)}\sqrt{\sigma(t)}||_1=\mathrm{Tr}(|\sqrt{\rho(t)}\sqrt{\sigma(t)}|)$ can then be written as
$$F(\rho,\sigma;t)=\mathrm{Tr}(\sqrt{\rho(t)}\sqrt{\sigma(t)}U)$$ for some $U$. Inserting the resolution of identity and then using a Cauchy-Schwarz inequality $|\mathrm{Tr}( A^\dagger B)|\leq \sqrt{\mathrm{Tr}(A^\dagger A)\mathrm{Tr}(B^\dagger B)}$ yields
\begin{aligned}
F(\rho,\sigma;t)=&\int dx \mathrm{Tr}(\sqrt{\rho(t)}|x\rangle\langle x|\sqrt{\sigma(t)}U)\\
=&\int dx \mathrm{Tr}(\sqrt{\rho(t)}\sqrt{|x\rangle\langle x|}\sqrt{|x\rangle\langle x|}\sqrt{\sigma(t)}U)\\
&\leq \int dx \mathrm{Tr}(\sqrt{\rho(t)}\sqrt{|x\rangle\langle x|}\sqrt{|x\rangle\langle x|}\sqrt{\rho(t)})\mathrm{Tr}(\sqrt{|x\rangle\langle x|}\sqrt{\sigma(t)}UU^\dagger \sqrt{\sigma(t)}\sqrt{|x\rangle\langle x|})\\
&= \int dx \langle x|\rho(t)|x\rangle\langle x|\sigma(t)|x\rangle.
\end{aligned} For this we simply needed that $\rho$, $\sigma$, and $|x\rangle\langle x|$ are Hermitian and that $\mathrm{Tr}(\sqrt{\rho}|x\rangle\langle x|\sqrt{\sigma}U)\leq |\mathrm{Tr}(\sqrt{\rho}|x\rangle\langle x|\sqrt{\sigma}U)|$.
Next, we expand $|\psi\rangle=\int dx \psi(x)|x\rangle$ in our chosen basis and use $e^{it \hat{\mathbf{B}}}|x\rangle=|x-t\rangle$ to compute
$$\langle x|\rho(t)|x\rangle=\int dy p_i (y) \langle x-yt|\psi\rangle\langle \psi|x-yt\rangle=\int dy p_i(y) |\psi(x-yt)|^2$$ and $$\langle x|\sigma(t)|x\rangle=\int dz p_j(z) |\psi(x-zt)|^2.$$
Putting these together, our inequality becomes
$$F(\rho,\sigma;t)\leq \int dx \int dy\int dz p_i(y) p_j(z) |\psi(x-yt)|^2 |\psi(x-zt)|^2$$ and the goal is to prove that the integral on the right-hand side vanishes in the limit of $t\to \infty$ when $p_i(x)$ and $p_j(x)$ have nonoverlapping compact support.
Now is where the hand-waving begins, but I'm sure someone can make this rigorous. Define the support of $p_i(x)$, $p_j(x)$, and $|\psi(x)|^2$ respectively as $S_i$, $S_j$, and $S_\psi$. The integrand is only nonzero for regions where $y\in S_i$, $z\in S_j$, $x-yt\in S_\psi$, and $x-zt\in S_\psi$. Since $S_i$ and $S_j$ are disjoint, we must have $|y-z|>0$ for the integrand to be nonzero (ignoring regions of zero measure). This means that the the integrand can only be nonzero if $S_\psi$ has two "spaced apart" regions with $x^\prime\equiv x-yt\in S_\psi$ and $x^\prime+(y-z)t\in S_\psi$. In the limit of $t\to\infty$, the integrand is only nonzero if there are two "infinitely spaced apart" regions $S_\psi$.
So, for any state $|\psi\rangle$ whose support in the Fourier transform basis of $\hat{\mathbf{B}}$ does not have two regions that are infinitely far apart, the relationship
$$\lim_{t\to\infty} F(\rho,\sigma;t)=0$$ holds.
Best Answer
Starting from a coherent state, which possesses a positive Wigner function and thus has $N_{|\alpha\rangle} = 0$, and evolving $|\alpha\rangle$ to a final state $|\psi\rangle$, for which $N_{|\psi\rangle} \geq 0$, indicates that some non-classicality, as evidenced by the emergence of negativity, was generated during the unitary process. This suggests that there is no classical theory capable of fully explaining the statistics of your data. However, caution is warranted as the negativity is not a necessary and sufficient condition for non-classicality.
For example, in the paper "Quantum catalysis in cavity QED", the authors discuss the possibility of generating non-classical states of light using a specific mechanism. They employ two figures of merit to characterise these states: the negativity of the Wigner function and second-order coherence. While neither figures of merit is individually necessary and sufficient to define non-classicality, their combined use is particularly effective for the task at hand.
Typically, the Wigner logarithmic negativity (WLN) is used instead of the standard negativity. This is defined as
$$\mathsf{W} \left ( \rho \right) := \log \left( \int \! \mathrm{d} x \,\mathrm{d} p \, \left| W_\rho \left(x,p \right) \right| \right),$$
where $W_\rho$ is the Wigner function of the state $\rho$. This quantity is related to the negativity via
$$\mathsf{W} = \log(N+1).$$
The WLN is an additive monotone (take a look at Appendix C of arXiv.1804.06763) since the Wigner function of separable states can be factorised. The aim of this function is to quantify a measure of non-Gaussianity, which has full monotonicity under Gaussian protocols.
Finally and more towards your last question, I would look at the average of WLN after $U[\theta(t)]$. This indeed would tell us the average of non-classicality generated by this protocol.