Here we briefly explore the impact of the presence of \(SO\)  interaction. The \(SO\)  induced correction to the scalar propagation constant, denoted \(\delta \beta ^{2(SO)}_{(l,m)}\) , is of the order of \(l\Delta /a^2\)  [1]}, [2]}, [3]}, being of opposite sign for the \(-l\)  state as compared to the \(l\)  state. This then implies that the diagonal elements of the \(2~x~2\)  symmetric matrix in the \((l,-l)\)  subspace (see Section 3.1 and Eq. A. 20), hitherto equal, will tend to become unequal, the change for one diagonal element being of the opposite sign compared to the other. Consequently, this can affect the accuracy of the derived results. In order for our results to be as accurate as possible, we therefore impose the stringent condition that the splitting of the degenerate spatial modes, \(2\gamma _{l,m}\)  due to ellipticity \(\epsilon \)  (Eqs. 21-23) far exceeds the splitting due to the \(SO\)  interaction, i.e., the \((l,-l)\)  degeneracy is still considered broken primarily through the elliptic effects.
A simple analysis based on the requirement \(\delta \beta ^{2(SO)}_{(l,m)}<<\gamma _{l,m}\)  as well as the sufficiency of the first order perturbation theory considered in our derivations (see Appendix B, [4]}), yields, in an approximate manner, bounds on \(\epsilon \)
\(\Big (\frac{lV}{2k^2a^2 n_1^2}\Big )^{1/l}\Big (\frac{2}{V}\Big )<< \epsilon <<\Big (\frac{2}{V}\Big ).\) 
