After inserting the above series in Eq. 5, taking the necessary scalar products to solve in different perturbation orders [1]}, [2]}, [3]}, we obtain the analytic expressions for the mixing coefficients and the propagation constant corrections:
\(a_{(l,m)(n,k)}^{(1)}=\frac{\epsilon \delta H_{(n,k)(l,m)}}{(\beta _{l,m}^2-\beta _{n,k}^2)},\) 
\(a_{(l,m)(n,k)}^{(2)}=\epsilon ^2\bigg ( \sum _{r,s}\frac{\delta H_{(n,k)(r,s)}\delta H_{(r,s)(l,m)}}{(\beta _{l,m}^2-\beta _{n,k}^2)(\beta _{l,m}^2-\beta _{r,s}^2)}-\frac{\delta H_{(n,k)(l,m)}\delta H_{(l,m)(l,m)}}{(\beta _{l,m}^2-\beta _{n,k}^2)^2}\bigg ),\) 
\(\beta ^{2(1)}_{l,m}=\epsilon \delta H_{(l,m)(l,m)},\) 
\(\beta ^{2(2)}_{l,m}=\epsilon ^{2}\sum _{n,k}\frac{\delta H_{(l,m)(n,k)}\delta H_{(n,k)(l,m)}}{(\beta _{l,m}^2-\beta _{n,k}^2)}.\) 
