Multiple orthogonal polynomials with respect to two weights on the step-line are considered. A connection between different dual spectral matrices, one banded (recursion matrix) and one Hessenberg, respectively, and the Gauss–Borel factorization of the moment matrix is given. It is shown a hidden freedom exhibited by the spectral system related to the multiple orthogonal polynomials. Pearson equations are discussed, a Laguerre–Freud matrix is considered, and differential equations for type I and II multiple orthogonal polynomials, as well as for the corresponding linear forms are given. The Jacobi–Piñeiro multiple orthogonal polynomials of type I and type II are used as an illustrating case and the corresponding differential relations are presented. A permuting Christoffel transformation is discussed, finding the connection between the different families of multiple orthogonal polynomials. The Jacobi–Piñeiro case provides a convenient illustration of these symmetries, giving linear relations between different polynomials with shifted and permuted parameters. We also present the general theory for the perturbation of each weight by a different polynomial or rational function aka called Christoffel and Geronimus transformations. The connections formulas between the type II multiple orthogonal polynomials, the type I linear forms, as well as the vector Stieltjes–Markov vector functions is also presented. We illustrate these findings by analyzing the special case of modification by an even polynomial.

A spectral Favard theorem is derived for bounded banded lower Hessenberg matrices that possess a positive bidiagonal factorization. The rich knowledge concerning the spectral and factorization properties of oscillatory matrices forms the basis of this theorem, which is formulated in terms of sequences of multiple orthogonal polynomials of types I and II, associated with a set of positive Lebesgue-Stieltjes measures. Additionally, a multiple Gauss quadrature is established, and the corresponding degrees of precision are determined. The spectral Favard theorem finds application in the context of Markov chains with transition matrices having (p + 2) diagonals, extending beyond the birth and death scenario, while still maintaining a positive stochastic bidiagonal factorization. In the finite case, the Karlin–McGregor spectral representation is provided, along with the demonstration of recurrent behavior and explicit expressions for the stationary distributions in terms of the orthogonal polynomials. Analogous results are obtained for countably infinite Markov chains. In this case, the Markov chain may not be recurrent, and its characterization is expressed in relation to the first measure. The ergodicity of the Markov chain is explored, taking into consideration the presence of a mass at 1, which corresponds to eigenvalues associated with the right and left eigenvectors.