  We consider the type I multiple orthogonal polynomials (MOPs) $(A_{n,m},
B_{n,m})$ and type II MOPs $P_{n,m}$, satisfying non-hermitian orthogonality
with respect to the weight $e^{-z^3}$ on two unbounded contours on $\mathbb C$.
Under the assumption that $$ n,m \to \infty, \quad \frac{n}{n+m}\to \alpha \in
(0, 1) $$ we find the detailed asymptotics of the MOPs, and describe the phase
transitions of this limit behavior as a function of $\alpha$. This description
is given in terms of vector critical measures, which are saddle points of an
energy functional comprising both attracting and repelling forces. These
critical measures are characterized by a cubic equation (spectral curve), and
their components $\mu_j$ live on trajectories of a canonical quadratic
differential $\varpi$ on the Riemann surface of this equation, which was object
of study in our previous paper [Adv. Math. 302 (2016), 1137--1232].
  The asymptotic zero distribution of the polynomials $A_{n,m}$ and $P_{n,m}$
are given by appropriate combinations of the components of the vector critical
measure. However, in the case of the zeros of $B_{n,m}$ the behavior is totally
different, and can be described in terms of the balayage of $\mu_2 - \mu_3$
onto certain curves on the plane. These curves are constructed with the aid of
$\varpi$, and their topology has three very distinct characters, depending on
the value of $\alpha$, and are obtained from the critical graph of $\varpi$.
  Once the trajectories and vector critical measures are studied, the main
asymptotic technical tool is the analysis of a $3\times 3$ Riemann-Hilbert
problem characterizing the MOPs.
  We illustrate our findings with results of several numerical experiments, and
formulate some conjectures and empirical observations based on these
experiments.
