Here is an example which answers several of your questions. It uses 2 different $7 \times 7$ matrices $A_G$ and $A_H$ whose entries are $0$ and $1$. These zeros are OK for your question because $A_G^4$ and $A_H^4$ are positive matrices (they are adjacency matrices of graphs each with an odd cycle)
Above are two graphs. Their respective adjacency matrices share some but not all eigenvalues. They do have the same largest eigenvalue and share an eigenvector for that eigenvalue.
The characteristic polynomial of the adjacency matrix of $G$ is $$(x+1)^2(x^2-3)(x^3-2x^2-3x+2)$$ while that of $H$ is $$(x-1)(x+1)(x^2+2x-1)(x^3-2x^2-3x+2)$$ The Perron-Frobenius eigenvalue of each comes from the shared cubic whose largest root is $\lambda \approx 2.8136.$ An eignvector for the adjacency matrix (or any other vector in $\mathbb{R}^7$ in this case) can be viewed as a weighting of the vertices of the graph. This is illustrated above where $a=\lambda-1$ and $b=\lambda^2-2\lambda-1$. This particular vector has entries which sum to $\lambda^2+1$. Dividing through by this and simplifying makes the entries $\frac{1}{\lambda^2+1}=\frac{6+\lambda-\lambda^2}{8}$ four times, $\frac{\lambda-1}{\lambda^2+1}=\frac{\lambda-2}{4}$ twice, and $\frac{\lambda^2-2\lambda-1}{\lambda^2+1}=\frac{\lambda^2-2\lambda-2}{2}$ once. Hence $H(x)=\mathbb{Z}[\frac{1}{2},\frac{\lambda}{8},\frac{\lambda^2}{8}].$ So the $\mathcal{B}_k$ for $A_G$ also contains $A_H^k$ as well as all the matrices $A_G^jA_H^{k-j}$
So I don't see any hope of classifying $\mathcal{B}_k$. I'm sure that there are also examples which don't have the same eigenvector and examples where there is a matrix $M$ with Perron-Frobenius eigenvalue $\lambda^k$ eigenvalue but such that $M$ is not the $k$th power of any matrix with $\lambda$ as an eigenvalue. Here I wanted an example that was small, simple, and did not involve two matrices with exactly the same spectrum.
I claim that any algebraic integer $\lambda$ is the eigenvalue of a nonnegative matrix with integer coefficients. This answer relies on Doug Lind's answer here. Let $\lambda_1$, $\lambda_2$, ..., $\lambda_r$ be the Galois conjugates of $\lambda$. Let $\ell = \max(|\lambda_i|)$. Choose an integer $L$ large enough that $L^n \cdot (L-2)/(L-1) > \ell^n \cdot \ell/(\ell-1)$ for all positive integers $n$. I claim there will be a nonnegative integer matrix with eigenvalues $(\lambda_1, \lambda_2, \ldots, \lambda_r, L,L,\ldots, L)$ where there are $r$ copies of $L$.
Let $\mu_1$, $\mu_2$, ..., $\mu_{2r}$ be $(\lambda_1, \lambda_2, \ldots, \lambda_r, L,L,\ldots, L)$. According to Lind's answer, all we have to check is that, for all $n>0$,
$$\frac{1}{n} \sum_k \sum_{d|n} \mu(n/d) \mu_k^d$$
is a nonnegative integer.
It's an integer: Let $g(x) = \prod_{i=1}^{2r} (x-\mu_i)$. Since $\lambda$ is an algebraic integer, $g$ is a monic polynomial with integer coefficients. So there is a $(2r) \times (2r)$ integer matrix $B$ with characteristic polynomial $g$ -- the companion matrix.
As explained on the second paper Lind links, if $g$ is the characteristic polynomial of an integer matrix, then there is a graph theoretic interpretation of the above quantity which is manifestly integral. (I gave a slightly incorrect description of this interpretation before; I'll stick to just citing for now.)
It's nonnegative
We have
$$\sum_{d|n} \mu(n/d) L^d \geq L^n - \sum_{-\infty < d<n} L^d = L^n - \frac{L^{n-1}}{1-1/L}=L^n \frac{L-2}{L-1}.$$
and
$$\left| \sum_{d|n} \mu(n/d) \lambda^d \right| \leq \sum_{- \infty < d \leq n} |\lambda|^d = \frac{|\lambda|^n}{1-1/|\lambda|} = |\lambda|^n \frac{|\lambda|}{|\lambda|-1}$$
So, for every $i$, $\sum_{d|n} \mu(n/d) L^d$ is a positive integer which is greater than the norm of $\sum_{d|n} \mu(n/d) \lambda_i^d$.
So the real part of
$$\sum_{i=1}^r \left( \sum_{d|n} \mu(n/d) L^d + \sum_{d|n} \mu(n/d) \lambda_i^d \right)$$
is nonnegative. By Galois symmetry, the sum is real, so it is a nonnegative real, as desired.
Disclaimer: I haven't read the paper Lind cites.
In response to Turbo's question: Any eigenvalue of an nonnegative integer matrix is an eigenvalue of a $(0,1)$ matrix. Choose $N$ larger than any matrix entry and replace each matrix entry $a_{ij}$ with an $N \times N$ block of the form
$$\begin{pmatrix}
1 & 1 & \cdots & 1 & 0 & 0 \cdots & 0 \\
1 & 1 & \cdots & 1 & 0 & 0 \cdots & 0 \\
1 & 1 & \cdots & 1 & 0 & 0 \cdots & 0 \\
\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \\
1 & 1 & \cdots & 1 & 0 & 0 \cdots & 0 \\
\end{pmatrix}$$
where there are $a_{ij}$ columns of $1$'s.
If $(v_1, v_2, \dots, v_n)^T$ is an eigenvector of the old matrix, then $(v_1, v_1, \dots, v_1, v_2, v_2, \dots, v_2, \dots, v_n, v_n, \dots, v_n)^T$ is an eigenvector of the new matrix with the same eigenvalue.
Best Answer
The answer to the first question is "yes", see p 21 of
http://arxiv.org/pdf/math/0703532v1
On that same page, if you look at conditions (6-8), that should answer your second question with a little work...