## 3 Linearized System
### Slip Boundary Condition  
Either we have \(\alpha=0\), in which case no wave solution exists, or there is a non-trivial solution with a dispersion relation for the phase velocity as a function of wavenumber \(c\left(k\right)\) implied by the singularity of the matrix, \(\Delta=det\left(\mathbf{Z}\right)=0\). Expanding the determinant \(\Delta\) algebraically using the permutation rule gives rise to 720 complex exponential terms upon eliminating \(H\) through application of the kinematic condition on the upper surface. Thus, the determination of the dispersion relation through solving \(\Delta=0\) symbolically is not a trivial undertaking. In contrast, the computation is readily tractable numerically, apart from a particular difficulty in finding the zeros of \(\Delta\). Numerical computations with fixed precision arithmetic are unlikely to return exactly zero. Therefore, a small number below a given threshold is typically taken to be zero. When the determinant of a matrix approaches zero, its condition number becomes large and thus the matrix becomes ill-conditioned, meaning that computations could involve high levels of numerical error. Consequently, the threshold level required to identify zeros of \(\Delta\) can be difficult to determine a priori. This motivates our desire to seek closed form symbolic approximations to the dispersion relation. Such approximations usually start with the long wave limiting cases. Presuming that \(c\) is an analytic function of \(k\) allows us to write the Taylor's series expansion of \(\Delta\left(c,k\right)\):  
\[\Delta\left(c,k\right)=\sum\Delta_{n}\left(c,0\right)k^{n}. \tag{24}\]  
Truncating the series at \(n=N\) yields an approximate dispersion relation implicitly by requiring the approximate determinant to vanish. We note that, with the matrix \(\mathbf{Z_{slip}}\) as written, the first non-trivial contribution comes at \(n=3/2\), with  
\[\Delta_{3/2}=-8 iPr^{7/2}c^{4}\left(c^{2}-\left(G+M\right)Pr\right) \tag{25}\]  
having a quadruple root at \(c=0\) and real roots  
\[c=\pm\sqrt{\left(G+M\right)Pr}. \tag{26}\]  
These are the long wave asymptotic limits that demonstrate that the Marangoni effect is additive to the gravity wave effect in contributing to the phase velocity. The capillary effect, however, is classically weaker by \(O(k^{2})\)[12]. To isolate these modes, it is convenient to translate the phase velocity to the frame of reference of the leftward (or rightward) moving wave and consider \(\Delta\left(\delta c+\sqrt{\left(G+M\right)\Pr},k\right)=0\). Since, by construction, \(\Delta\) is analytic in both \(k\) and \(c\), it is useful to develop the \(k\)-power series of the \(c\)-truncation of \(\Delta\), i.e.,  
\[\Delta\left(\delta c+\sqrt{\left(G+M\right)\Pr},k\right)\approx\Delta\left( \sqrt{\left(G+M\right)\Pr},k\right)+\delta c\Delta^{\left(1\right)}\left( \sqrt{\left(G+M\right)\Pr},k\right)=0. \tag{27}\]  
Developing Equation (27) to \(O(k^{4})\) yields  
\[\delta c=q_{1}k+q_{2}k^{2}+q_{3}k^{3}, \tag{28}\]with  
\[q_{1} = -\frac{i}{6}\left(M-12\right)\,\mathrm{Pr},\] \[q_{2} = -\frac{\sqrt{\mathrm{Pr}}\left(-45+3 GK\left(5+2M\,\mathrm{Pr} \right)+K\left(30M+M^{2}+1440\,\mathrm{Pr}-240M\,\mathrm{Pr}+16M^{2}\,\mathrm{ Pr}\right)\right)}{90K\sqrt{G+M}},\] \[q_{3} = \frac{i}{7560K\left(G+M\right)}(4410\left(12-M\right)\,\mathrm{Pr }+204\mathrm{G}^{2}K\mathrm{MPr}^{2}\] \[+ GK\left(-2520\,\mathrm{Pr}\left(11+6\,\mathrm{Pr}\right)-42M \left(-4-85\,\mathrm{Pr}+156\mathrm{Pr}^{2}\right)+M^{2}\left(8+29\,\mathrm{ Pr}+1059\mathrm{Pr}^{2}\right)\right)\] \[+ K(-15120\,\mathrm{Pr}\left(1+65\,\mathrm{Pr}\right)+2520M\, \mathrm{Pr}\left(-17+91\,\mathrm{Pr}\right)\] \[\Delta_{3/2}=-8 iPr^{7/2}c^{4}\left(c^{2}-\left(G+M\right)Pr\right) \tag{25}\]0  
which is readily interpreted as a long wave instability occurring for a Marangoni number greater than the critical value \(M_{crit}=12\) for long waves \(k<<1\).