Conceptually, without any other input about the physics of the problem, the only way I can think to pick the QNs of the expanded subspace is to pick them randomly (perhaps some gaussian distribution of QNs). An alternative way could be the following, since you are doing a comparison calculation anyway:
1. Do a hybrid 2-site and 1-site TDVP calculation.
2. Then do a purely 1-site TDVP, but get the QNs for the subspace expansion from the hybrid calculation done in 1. that you will compare to.
Then, the calculations are being performed in the same MPS manifold, so it seems like the most fair comparison.
I have indeed heard that about DMRG, but had not tried it myself. It would be good for ITensor to support both 1-site and 2-site DMRG calculations, and we are also interested in supporting TDVP as soon as possible.