## 2 Method Description
### Iterative Computation of the Taylor Polynomial  
The constant term of the Taylor approximation of \(\lambda_{i}(\mathbf{q},\mathbf{y})\) and \(\boldsymbol{\phi}_{i}(\mathbf{q},\mathbf{y})\) is simply \(\lambda_{i}^{0}=\lambda_{i}(\mathbf{q}_{0},\mathbf{y}_{0})\) and \(\boldsymbol{\phi}_{i}^{0}=\boldsymbol{\phi}_{i}(\mathbf{q}_{0},\mathbf{y}_{0})\), respectively, which can be computed by solving the direct eigenvalue problem (2) for \((\mathbf{q},\mathbf{y})=(\mathbf{q}_{0},\mathbf{y}_{0})\). The coefficients of the first-order Taylor polynomial of the \(i\)-th eigenvalue \(\lambda_{i}\) and eigenvector \(\boldsymbol{\phi}_{i}\) are obtained by substituting the formulas (3) and (4) for \(p=1\) in (2) and retaining only the constant and linear terms. Therefore, the following equation is obtained  
\[ \begin{split}\mathbf{K}^{0}\boldsymbol{\phi}_{i}^{0}+\sum_{| \alpha|+|\beta|=1}\Bigl{[}\partial_{\mathbf{q}}^{\alpha}\partial_{\mathbf{y}}^{ \beta}\mathbf{K}(\mathbf{q}_{0},\mathbf{y}_{0})\boldsymbol{\phi}_{i}^{0}+ \mathbf{K}^{0}\partial_{\mathbf{q}}^{\alpha}\partial_{\mathbf{y}}^{\beta} \boldsymbol{\phi}_{i}(\mathbf{q}_{0},\mathbf{y}_{0})\Bigr{]}\delta\mathbf{q}^{ \alpha}\delta\mathbf{y}^{\beta}=\lambda_{i}^{0}\mathbf{M}^{0}\boldsymbol{\phi}_ {i}^{0}+\\ \sum_{|\alpha|+|\beta|=1}\Bigl{[}\partial_{\mathbf{q}}^{\alpha} \partial_{\mathbf{y}}^{\beta}\lambda_{i}(\mathbf{q}_{0}.\mathbf{y}_{0}) \mathbf{M}^{0}\boldsymbol{\phi}_{i}^{0}+\lambda_{i}^{0}\partial_{\mathbf{q}}^{ \alpha}\partial_{\mathbf{y}}^{\beta}\boldsymbol{\lambda}(\mathbf{q}_{0},\mathbf{y} _{0})\boldsymbol{\phi}_{i}^{0}+\lambda_{i}^{0}\mathbf{M}^{0}\partial_{\mathbf{q}}^{ \alpha}\partial_{\mathbf{y}}^{\beta}\boldsymbol{\phi}_{i}(\mathbf{q}_{0},\mathbf{y} _{0})\Bigr{]}\delta\mathbf{q}^{\alpha}\delta\mathbf{y}^{\beta}.\end{split}  \tag{5}\]Since the constant terms \(\mathbf{K}^{0}\mathbf{\phi}_{i}^{0}=\lambda_{i}^{0}\mathbf{M}^{0}\mathbf{\phi}_{i}^{0}\) cancel due to Equation (2), the coefficients of the right-hand and left-hand sides can be equated and solved. There are in total \(n_{q}+n_{y}\) multi-indices \(\alpha,\beta\) such that \(|\alpha|+|\beta|=1\). For any \(\alpha,\beta\) the following system of equations is obtained  
\[\left[ \begin{array}{cc}\lambda_{i}^{0}\mathbf{M}^{0}-\mathbf{K}^{0}&\mathbf{ M}^{0}\mathbf{\phi}_{i}^{0}\end{array} \right]\left\{ \begin{array}{cc} \tilde{\mathbf{q}}_{\alpha}^{a}\tilde{\mathbf{\phi}}_{i}^{b}\mathbf{\phi}_{i}(\mathbf{\mathrm{ q}}_{0},\mathbf{\mathrm{y}}_{0})\\ \tilde{\mathbf{q}}_{\alpha}^{a}\tilde{\mathbf{\phi}}_{\beta}^{b}\lambda_{i}(\mathbf{ \mathrm{q}}_{0},\mathbf{\mathrm{y}}_{0})\end{array} \right\}=\left[\tilde{\mathbf{q}}_{ \alpha}^{a}\tilde{\mathbf{q}}_{\beta}^{b}\mathbf{K}(\mathbf{\mathrm{q}}_{0},\mathbf{ \mathrm{y}}_{0})-\lambda_{i}^{0}\tilde{\mathbf{q}}_{\alpha}^{a}\tilde{\mathbf{\phi}}_{ \beta}^{b}\mathbf{M}(\mathbf{\mathrm{q}}_{0},\mathbf{\mathrm{y}}_{0})\right]\mathbf{\phi}_ {i}^{0}. \tag{6}\]  
Since there are \(n_{\mathrm{dof}}\) equations and \(n_{\mathrm{dof}}+1\) unknowns, the system is under-determined and can be solved. An efficient way consists of computing the pseudo-inverse of \(\left[ \begin{array}{cc}\lambda_{i}^{0}\mathbf{M}^{0}-\mathbf{K}^{0}&\mathbf{ M}^{0}\mathbf{\phi}_{i}^{0}\end{array} \right]\), which does not depend on \(\alpha,\beta\) and has to be computed just once for all the systems. If the coefficients of \(\lambda_{i}^{1}(\mathbf{\mathrm{q}},\mathbf{\mathrm{y}})\) and \(\mathbf{\phi}_{i}^{1}(\mathbf{\mathrm{q}},\mathbf{\mathrm{y}})\) are known, it is possible to compute the coefficients of \(\lambda_{i}^{2}(\mathbf{\mathrm{q}},\mathbf{\mathrm{y}})\) and \(\mathbf{\phi}_{i}^{2}(\mathbf{\mathrm{q}},\mathbf{\mathrm{y}})\) by simply substituting (3) and (4) for \(p=2\) into (2) and neglecting the terms of order higher than 2. The obtained equations can be used to compute the coefficients of the second-order terms, which are the only unknowns.