  The singular values squared of the random matrix product $Y = G_r G_{r-1}
\cdots G_1 (G_0 + A)$, where each $G_j$ is a rectangular standard complex
Gaussian matrix while $A$ is non-random, are shown to be a determinantal point
process with correlation kernel given by a double contour integral. When all
but finitely many eigenvalues of $A^*A$ are equal to $bN$, the kernel is shown
to admit a well-defined hard edge scaling, in which case a critical value is
established and a phase transition phenomenon is observed. More specifically,
the limiting kernel in the subcritical regime of $0<b<1$ is independent of $b$,
and is in fact the same as that known for the case $b=0$ due to Kuijlaars and
Zhang. The critical regime of $b=1$ allows for a double scaling limit by
choosing $b = (1-\tau/\sqrt{N})^{-1}$, and for this the critical kernel and
outlier phenomenon are established. In the simplest case $r=0$, which is
closely related to non-intersecting squared Bessel paths, a distribution
corresponding to the finite shifted mean LUE is proven to be the scaling limit
in the supercritical regime of $b>1$ with two distinct scaling rates. Similar
results also hold true for the random matrix product $T_r T_{r-1} \cdots T_1
(G_0 + A)$, with each $T_j$ being a truncated unitary matrix.
