$\rho_{\mu}$ are independent field variables that follow [@PhysRev.131.2326]

$$\label{fieldeqns}
  \partial_{\mu} \rho_{\mu\nu}-m^2 \rho_{\nu}+ 2i \kappa m \rho_{\mu} F_{\mu\nu}=0,\quad\quad
  \rho_{\mu\nu} = \partial_{\mu} \rho_{\nu} - \partial_{\nu} \rho_{\mu}.$$

A variation of the Lagrangian with respect to the field $\rho_{\mu}$ give us the equations of motion that in the momentum space read

$$\left((p_{\mu}^2  + m^2)\delta_{\mu\nu} -p_{\mu} p_{\nu}  - 2  i \kappa m F_{\mu \nu}\right)\rho_{\mu} = 0.$$

Thus, the boson propagator is

$$\label{propagator}
D_{\mu\nu}^{-1}=((p_{\mu}^2  + m^2)\delta_{\mu\nu}-p_{\mu} p_{\nu}  - 2  i \kappa m F_{\mu \nu}).$$

Considering the magnetic field uniform, constant and in $p_3$ direction $\textbf{B}=B\textbf{e}_3$ one can start from Eq. (\[fieldeqns\]) and obtain the generalized Sakata-Taketani hamiltonian for the six component wave equation of the system [@PhysRev.131.2326; @PhysRevD.89.121701] following the same procedure of Ref. [@PhysRev.131.2326]. The hamiltonian reads

$$\label{hamiltonian}
H = \sigma_3 m + (\sigma_3 + i \sigma_2) \frac{\textbf{p}^2}{2 m} -
    i \sigma_2 \frac{(\textbf{p}\cdot\textbf{S})^2}{m}
    -(\sigma_3 - i \sigma_2) \kappa \textbf{S} \cdot \textbf{B},$$

with $\textbf{p}=(p_{\perp},p_3)$ and $p_{\perp}=p_1^2 + p_{2}^2$. $\sigma_{i}$ are the $2\times2$ Pauli matrices, $S_{i}$ are the $3\times3$ spin-1 matrices in a representation in which $S_3$ is diagonal and $\textbf{S} = \{S_1,S_2,S_3\}$[^5].

The eigenvalues of (\[hamiltonian\]) are $$\varepsilon(p_3,p_{\perp}, B,s)=\sqrt{m^2+p_3^2+p_{\perp}^2-2\kappa s B\sqrt{p_{\perp}^2+m^2}},\label{spectrum}$$

where $s=0, \pm 1$ are the spin eigenvalues. Although we are dealing with neutral bosons, 