  For a pair of coupled rectangular random matrices we consider the squared
singular values of their product, which form a determinantal point process. We
show that the limiting mean distribution of these squared singular values is
described by the second component of the solution to a vector equilibrium
problem. This vector equilibrium problem is defined for three measures with an
upper constraint on the first measure and an external field on the second
measure. We carry out the steepest descent analysis for a 4 $\times$ 4
matrix-valued Riemann-Hilbert problem, which characterizes the correlation
kernel and is related to mixed type multiple orthogonal polynomials associated
with the modified Bessel functions. A careful study of the vector equilibrium
problem, combined with this asymptotic analysis, ultimately leads to the
aforementioned convergence result for the limiting mean distribution, an
explicit form of the associated spectral curve, as well as local Sine, Meijer-G
and Airy universality results for the squared singular values considered.
