We study the semigroups of random Schrödinger operators of the form $\widehat{H}f=-\frac12f''+(V+ξ)f$, where $f:I\to\mathbb F^r$ ($\mathbb F=\mathbb R,\mathbb C,\mathbb H$) are vector-valued functions on a possibly infinite interval $I\subset\mathbb R$ that satisfy a mix of Robin and Dirichlet boundary conditions, $V$ is a deterministic diagonal potential with power-law growth at infinity, and $ξ$ is a matrix white noise. Our main result consists of Feynman-Kac formulas for trace moments of the form $\mathbf E[\prod_{k=1}^n\mathrm{Tr}[\mathrm e^{-t_k\widehat{H}}]]$ ($n\in\mathbb N$, $t_k>0$). One notable example covered by our main result consists of the multivariate stochastic Airy operator (SAO) of Bloemendal and Virág (Ann. Probab., 44(4):2726-2769, 2016), which characterizes the soft-edge eigenvalue fluctuations of critical rank-$r$ spiked Wishart and GO/U/SE random matrices. As a corollary of our main result, we prove that if $V$'s growth is at least linear (this includes the multivariate SAO), then $\widehat{H}$'s spectrum is number rigid in the sense of Ghosh and Peres (Duke Math. J., 166(10):1789-1858, 2017). Together with the rigidity of the scalar SAO, this completes the characterization of number rigidity in the soft-edge limits of Gaussian $β$-ensembles and their finite-rank spiked versions.