Let $d_i=I\{y_{i1}>0\;and\;y_{i2}>0\}$ be equal to 1 if the statement inside the brackets is true and 0 otherwise. Then I believe the likelihood would be
$$
\mathscr{L} = \prod_{i=1}^n \boldsymbol{\phi}_{\Sigma}\bigg(
\begin{matrix}
y_{i1}-x'_{i1}\beta\\
y_{i2}-x'_{i2}\beta
\end{matrix}\;
\bigg)^{d_i}\bigg[\Phi_{\Sigma}\bigg(
\begin{matrix}
-x'_{i1}\beta\\
-x'_{i2}\beta
\end{matrix}\;
\bigg)\frac{\phi_{\sigma_1}(y_{i1}-x'_{i1}\beta)}{\Phi{\sigma_1}(-x'_{i1}\beta)}\bigg]^{(1-d_i)}
$$
Where $\Sigma$ is the entire covariance matrix, $\boldsymbol{\phi}_{\Sigma}$ the centered bi-variate normal pdf of covariance $\Sigma$ ,$\Phi_{\sigma_2}$ the centered normal cdf of variance $\sigma_2^2$, and $\phi_{\sigma_1}$ the centered normal pdf of variance $\sigma_1^2$. We can also write it as
$$
\mathscr{L} = \prod_{i=1}^n \boldsymbol{\phi}_{\Sigma}\bigg(
\begin{matrix}
y_{i1}-x'_{i1}\beta\\
y_{i2}-x'_{i2}\beta
\end{matrix}\;
\bigg)^{d_i}\bigg[\Phi_{\Sigma}\bigg(
\begin{matrix}
-x'_{i1}\beta\\
-x'_{i2}\beta
\end{matrix}\;
\bigg)\frac{1}{\sigma_1}\frac{\phi_{1}(\frac{y_{i1}-x'_{i1}\beta}{\sigma_1})}{\Phi_{1}(\frac{-x'_{i1}\beta}{\sigma_1})}\bigg]^{(1-d_i)}
$$
The term $\Phi_{\Sigma}\bigg(
\begin{matrix}
-x'_{i1}\beta\\
-x'_{i2}\beta
\end{matrix}\;
\bigg)$ is the joint probability that both $y_{i1}$ and $y_{i2}$ are less than 0, which is the mass point creating bias in MV least squared regression. Note that if the covariance $\sigma_{12}=0$,
$$\Phi_{\Sigma}\bigg(
\begin{matrix}
-x'_{i1}\beta\\
-x'_{i2}\beta
\end{matrix}\;
\bigg)=\Phi{\sigma_1}(-x'_{i1}\beta)\Phi{\sigma_2}(-x'_{i2}\beta)$$
and
$$\boldsymbol{\phi}_{\Sigma}\bigg(
\begin{matrix}
y_{i1}-x'_{i1}\beta\\
y_{i2}-x'_{i2}\beta
\end{matrix}\;
\bigg)=\phi_{\sigma_1}(y_{i1}-x'_{i1}\beta)\phi_{\sigma_2}(y_{i2}-x'_{i2}\beta)$$
Unfortunately I do not have a source for you. Let me know if there are any problems or if further explanation is needed. Note that if the question goes beyond just how to write the likelihood to actually maximizing it via computational methods, things get much more complicated.
I am not so sure about the last two cases. I am assuming you mean that the error terms in the latent model are iid $f_1,f_2$. First I restate $y_{1i}$ and $y_{2i}$ as piece-wise functions of $y^*_{1i}$ and $y^*_{2i}$.
$$
y_{1i} = \left\{
\begin{array}{lr}
y^*_{1i} & if\;\; y^*_{2i}>0\\
0 & if\;\; y^*_{2i} \leq 0
\end{array}
\right.\;\;\;\;\;\;\;\;
y_{2i} = \left\{
\begin{array}{lr}
1 & if\;\; y^*_{2i}>0\\
0 & if\;\; y^*_{2i} \leq 0
\end{array}
\right.
$$
(or at least that's my interpretation of your notation). Then I create some indicator function. $d_i=I\{y_{2i}>0\}$.
- We observe $(y_{1},y^*_{2})$:
In this case notice that when $d_i=1$, the conditional probability of observing any pair $(y_{1i},y^*_{2i})$ is simply the conditional joint pdf, which I will call $f_{1,2}$, as a function of the $(y_{1i},y^*_{2i})$,$(x_{1i},x_{2i})$ and $(\beta_{1},\beta_{2})$.
When $d_i=0$ the joint probability would be
$$
Pr(Y_{1i}=0\,|\,Y^*_{2i}\leq0)\frac{f_2(Y^*_{2i}-x'_{2i}\beta_2)}{Pr(Y^*_{2i}\leq 0)}=Pr(Y^*_{2i}\leq 0)\frac{f_2(y^*_{2i}-x'_{2i}\beta_2)}{Pr(Y^*_{2i}\leq 0)}=f_2(y^*_{2i}-x'_{2i}\beta_2)
$$
The term $\frac{f_2(y^*_{2i}-x'_{2i}\beta_2)}{Pr(y^*_{2i}\leq 0)}$ is the left-of-zero truncated distribution of $u_{2i}$. So we write the likelihood:
$$
\mathscr{L} = \prod_{i=1}^n f_{1,2}\bigg(
\begin{matrix}
y_{i1}-x'_{i1}\beta\\
y_{i2}-x'_{i2}\beta
\end{matrix}\;
\bigg)^{d_i}f_2(y^*_{2i}-x'_{2i}\beta_2)^{(1-d_i)}
$$
- We observe $(y_{1},y_{2})$:
So when $d_i=1$ the joint probability is
$$
Pr(Y^*_{2i}>0\cap Y_{1i}=y_{1i} )=Pr(Y^*_{2i}>0\,|\, Y_{1i}=y_{1i})Pr(Y_{1i}=y_{1i})=
$$
$$
Pr(U_{2i}>-x'_{2i}\beta_2\,|\, Y_{1i}=y_{1i})Pr(Y_{1i}=y_{1i})=
$$
$$
\bigg(\int_{-x'_{2i}\beta_2}^{\eta}\frac{f_{12}(y_{1i}-x'_{1i},t)}{f_1(y_{1i}-x'_{1i}\beta_1)}dt\bigg)f_1(y_{1i}-x'_{1i}\beta_1)=\int_{-x'_{2i}\beta_2}^{\eta}f_{12}(y_{1i}-x'_{1i}\beta_1,t)dt
$$
for some upper bound $\eta$ which may or may not be $\infty$.
When $d_i=0$ the only information we have to construct a probability is
$$
Pr(y_i \leq 0)=F_2(-x'_{2i}\beta_2)
$$
So the Likelihood will be
$$
\mathscr{L} = \prod_{i=1}^n \bigg[\int_{-x'_{2i}\beta_2}^{\eta}f_{12}\bigg(\begin{matrix}
y_{i1}-x'_{i1}\beta\\
t
\end{matrix}\;
\bigg)dt\bigg]
^{d_i}F_2(-x'_{2i}\beta_2)^{(1-d_i)}
$$