## 3 Results  
### Power Optimization of System under Inertial Load  
The local optima for dimensionless power, in the domain of \(\gamma\) and \(\xi_{C}\), can be determined from the zeros of the corresponding partial derivatives of the power expression. One approach to finding the analytical solution to the local power optima is to first find the relation between \(\gamma\) and optimal load coefficient from the zeros of \(\partial\overline{P}/\partial\xi_{C}\) (assuming \(\omega_{N}\) and \(B\) are constants):  
\[\xi_{C}=\sqrt{\gamma^{2}+\xi_{E}^{2}+\frac{Q\eta(2\xi_{E}+Q\eta-2\epsilon \gamma)}{(\epsilon^{2}+1)}} \tag{17}\]  
By applying Equation (17) in Equation (11), we find the expression for dimensionless power under the condition of optimal load, given by Equation (18). \[\overline{P}_{IL\_OptLoad}=\frac{Q\eta/4}{\sqrt{\gamma^{2}\big{(}\epsilon^{2} +1\big{)}(\epsilon^{2}+1)+\xi_{E}(\epsilon^{2}+1)(2Q\eta+\xi_{E}(\epsilon^{2 }+1))+\xi_{E}(\epsilon^{2}+1)+Q\eta}} \tag{18}\]  
For a more concise expression, an additional parameter, \(\xi\), is used in Equation (18). This additional parameter is defined by Equation (19). \[\xi=Q\frac{\gamma^{2}-(\eta+1)}{\gamma} \tag{19}\]  
The zeros of the derivative with respect to \(\gamma\), \(\partial\overline{P}_{IL\_OptLoad}/\partial\gamma\), then give the local optimum \(\gamma\). However, without assuming negligible \(\xi_{E}\), we find the analytical root derivation to be prohibitively complex, even for an analytical software tool. We therefore perform a numerical investigation to determine the model characteristics under optimal excitation frequency and load resistance. We find the locally optimal values of \(\gamma\) and \(\xi_{C}\) by finding the intersections of solutions to \(\partial\overline{P}_{IL}/\partial\xi_{C}=0\) and \(\partial\overline{P}_{IL}/\partial\gamma=0\). To compare our results with the previous literature, we can set \(\xi_{E}=0\) in Equation (18) and solve for \(\partial\overline{P}_{IL,OptLoad}/\partial\gamma=0\), in which case we arrive at the same solutions as [PERSON] et al. [3], under the restrictions of positive and real solutions. Under the same condition of \(\xi_{E}=0\), we also find a critical value of the effective electromechanical coupling, \(\eta_{crit}=4\zeta(1+\zeta)\), above which the number of real and positive solutions to \(\partial\overline{P}_{IL,OptLoad}/\partial\gamma=0\) increases from one to three, where one solution is a local minimum. For \(\eta\leq\eta_{crit}\) there is only a single maximum. In our continued investigation, we will use an equivalent critical value, expressed as a quality factor, given by Equation (20). \[Q_{crit}=\left(1+\sqrt{1+\eta}\right)/\eta \tag{20}\]  
\(Q_{crit}\) is used in the same way as \(\eta_{crit}\), i.e., for \(Q\leq Q_{crit}\), there is only one real maximum, and for \(Q>Q_{crit}\), there are two real maxima. For the case of \(Q\gg Q_{crit}\) and \(\xi_{E}=0\), the two possible local power maxima, with respect to \(\gamma\), tend toward the open-circuit (\(\beta\to\infty\), \(\gamma\to\sqrt{1+\eta}\)) and closed-circuit (\(\beta\to\infty\), \(\gamma\to 1\)) resonance frequencies, also called resonance and anti-resonance. In our investigation, we dub the solution tending toward resonance as \(\gamma_{opt\_R}\), and for the solution tending toward anti-resonance, we use \(\gamma_{opt\_AR}\). Correspondingly we use \(\xi_{C,opt\_R}\) and \(\xi_{C\_opt\_AR}\). The results from our numerical investigation show that for \(\xi_{E}>0\), there still exists a critical value of \(Q\) (or equivalently, \(\eta\)), above which we find three viable local optima, although it is significantly modified by the inclusion of \(\xi_{E}\). As we do not have the exact analytical expression for the new boundary, we perform our analysis by varying \(Q\) in factors of \(Q_{crit}\). For this reason, we introduce the parameter \(k\), defined by Equation (21).