## 2 Statement of the Modal Modification Problem
### Topology Optimization  
The objective of the topology optimization process is to reduce the thermoelastic dissipation of the microbeam resonator. The thermoelastic dissipative ratio \(\varphi_{i}\) can be calculated by dividing the imaginary part of the eigenvalue of the system into its real part [6; 43]:  
\[\varphi_{i}=2\bigg{|}\frac{\Im(\lambda_{i})}{\Re(\lambda_{i})}\bigg{|} \tag{8}\]  
where \(\Im(\bullet)\) and \(\Re(\bullet)\) denote the imaginary and real part of \((\bullet)\), respectively. An iterative bidirectional evolutionary structural optimization (BESO) method for thermoelastic dissipation of a first-order mode shape can be constructed as follows:  
\[ \begin{split}\text{\emph{minimize}}&\{\varphi_{1} \}\\ \text{\emph{subject to}}&(\lambda_{i}\mathbf{B}- \mathbf{A})\mathbf{\Phi}_{i}=\mathbf{0}\\ &\sum_{e=1}^{N_{e}}\Omega_{e}p_{e}-l_{v}\Omega_{0}\leq 0\\ & p_{e}=p_{\text{min}}\text{ or }1,\quad(e=1,2,...N)\end{split} \] (9a) where \[p_{e}\] is the design variable (elemental pseudo-density) and \[p_{\text{min}}\] is a minimum positive number near 0 to avoid singularity in the finite element analysis process. \[\Omega_{e}\] and \[\Omega_{0}\] are the elemental and global volumes, respectively, \[l_{v}\] is the volume constraint factor, and \[N_{e}\] is the total number of the finite element mesh. The combined governing equation (Equation ( 9a )) is the main constraint of the optimization problem, along with the standard BESO topology optimization volume constraint in Equation ( 9b ) and design variable constraint in Equation ( 9c ). According to the definition of thermoelastic dissipation (8) and the complex eigenvalue problem (7), the sensitivity of the objective function \(\varphi_{i}\) with respect to design variable \(p_{e}\) can be obtained through the chain rule as follows [44]:  
\[\frac{\partial\varphi_{i}}{\partial p_{e}}=2\frac{\Im\big{(}\frac{\partial \lambda_{i}}{\partial p_{e}}\big{)}\Re(\lambda_{i})-\Im(\lambda_{i})\Re\Big{(} \frac{\partial\lambda_{i}}{\partial p_{e}}\Big{)}}{\left[\Re(\lambda_{i}) \right]^{2}} \tag{10}\]  
where the sensitivity of the \(i\)th eigenvalue of the system can be obtained through  
\[\frac{\partial\lambda_{i}}{\partial p_{e}}=\mathbf{\
abla}_{i}^{\text{T}}\Big{(} \frac{\partial\mathbf{A}}{\partial p_{e}}-\lambda_{i}\frac{\partial\mathbf{B }}{\partial p_{e}}\Big{)}\mathbf{\Phi}_{i} \tag{11}\]  
where \(\mathbf{\
abla}_{i}\) is the left eigenvector of the system. The process of the gradient-based iterative optimization method can be discretized into several steps, as shown in Figure 2. The main task of the iterative optimization process is to obtain the sensitivity numbers of elements through the modal modification analysis process. Marked in blue in Figure 2, the solution of the eigenvalue problem and calculation of global sensitivity numbers are the most time- and computational effort-consuming steps of the modal modification analysis process. As mentioned before, the complex eigenvalue problem for large-scale models is difficult and time-consuming to solve. Due to the change in the material distribution of the structure during the iterative process of topology optimization, the system matrices and eigenvalue problems are altered in each iteration. According to Equations (10) and (11), the sensitivity of the objective function with respect to the design variable is dependent on the mode shape of the system. This leads to the fact that the optimization process requires repeated complex modal analysis of the system in each iteration, which will further prolong the existing time-consuming task. On the other hand, the observation of the sensitivity equation (Equation (10)) shows that for large-scale models with massive degrees of freedom, solving the sensitivity number of the whole system at one time could be difficult (or even impossible). For this, the sensitivity number of each element has to be solved at the element scale and aggregated to a global matrix.