MUSIC algorithm derivation
The piece I was missing was the distribution of the initial phase values $\varphi_1$ and $\varphi_2$. It is standard to assume that these are uniformly distributed [^]. This leads to:
$$
\mathbf{R_x} = \mathbf{E[xx^H]} = \mathbf{S \Lambda S^H}
$$
where
$\mathbf{\Lambda} = \begin{bmatrix} A_1^2 &amp; 0 \\ 0 &amp; A_2^2
\end{bmatrix}$
and $\mathbf{S} = \begin{bmatrix}\mathbf s_1 &amp; \mathbf s_2 \end{bmatrix}$ so that
$$
\mathbf{R_y} = \mathbf{S \Lambda S^H} + \sigma^2 \mathbf{I}.
$$

Now note that $\mathbf{R_y}$ is Hermitian symmetric and must possess an eigendecomposition with eigenvalues $\lambda_1 \geq \lambda_2 \geq \lambda_3 \geq \cdots \geq \lambda_L$. However, also note that $\mathbf{R_x}$ has two nonzero eigenvalues, say, $\tilde\lambda_1 \geq \tilde\lambda_2$. Since $\mathbf{R_y} = \mathbf{R_x} + \sigma^2 I$ it is easy to verify that $\lambda_1 = \tilde \lambda_1 + \sigma^2$, $\lambda_2 = \tilde \lambda_2 + \sigma^2$ and $\lambda_i = \sigma^2$ for $3\leq i \leq L$.

In the eigendecomposition of $\mathbf{R_y}$, let $\mathbf e_1$ and $\mathbf e_2$ be the eigenvectors corresponding to $\lambda_1$ and $\lambda_2$ and $\mathbf g_1, \ldots, \mathbf g_{L-2}$ be the eigenvectors corresponding to the eigenvalues $\lambda_3, \lambda_4,\ldots, \lambda_L$. One can now show that $\mbox{span}\{\mathbf e_1, \mathbf e_2\} = \mbox{span}\{\mathbf s_1, \mathbf s_2\}$. This is often called the "signal subspace" and the orthogonal subspace $\mbox{span}\{ \mathbf g_1,\ldots \mathbf g_{L-2}\}$ is called the "noise subspace."

Here's the punchline: To estimate $\omega_1$ and $\omega_2$ we need to solve for $\omega$ in $\mathbf s(\omega)^H \mathbf G \mathbf G^H \mathbf s(\omega) = 0$, where $\mathbf G = \begin{bmatrix} \mathbf g_1 &amp; \ldots &amp; \mathbf g_{L-2}\end{bmatrix}$, and $\mathbf s(\omega) = \begin{bmatrix} 1 &amp; e^{j\omega} &amp; e^{j2\omega} &amp; \ldots &amp; e^{j(L-1)\omega} \end{bmatrix}$.

Reference

[^] P. Stoica, R. Moses, Spectral Analysis of Signals 1st Ed. Ch. 4.