5 a GLYPH<c165> 0 2 4 6 8 0 2 4 6 8 10 12 14 16 x10 -5 1/Q eff Z$_{load}$ =9 GLYPH<c87> 91 GLYPH<c87> 240 GLYPH<c87> 554 GLYPH<c87> GLYPH<c87> I 2 $_{dc}$(mA$^{2}$) b 0 2 4 6 8 0 10 20 30 40 Z$_{load}$ =40 GLYPH<c109> H 22 GLYPH<c109> H 10 GLYPH<c109> H 5.6 nF 51pF 1.0 nF Z =L load Z =R load Z =C load Cooling Heating -Im ( ) GLYPH<c103> GLYPH<c98> Z 0 0.0 0.4 0.8 -100 -80 -60 -40 -20 0 20 40 -0.4 0.0 0.4 Im GLYPH<c103>$_{Z}$ -100 -80 -60 -40 -20 0 20 40 -Im (A ) GLYPH<c98> -2 c Refrigerator Heat pump x10 -5 Re GLYPH<c103>$_{Z}$ I 2 $_{dc}$(mA$^{2}$) of these linear fits, Im β is determined and is plotted against γ$_{Z}$ as symbols in Fig. 4c. As shown by solid and dashed lines in Fig. 4c, an excellent multiple linear regression fit of the data is obtained using the function Im β = Im ( γ$_{Z}$β$_{0}$ ), which yields the fit parameters β$_{0}$ = − 123 . 7 + 64 . 6 i A − 2 and R$_{beam}$ = 439 . 3 Ω. In order to gain insight in the operation of the heat en- gine, a schematic showing the rate equations for the flow of heat and work between electrical source, heat sink, en- gine beam and resonance mode is shown in Fig. 5a. All rates shown in this figure are time averaged, such that all variations which cancel each other during a single cycle are averaged out. The three modes of operation which were experimentally observed in Fig. 2d, 3a, and 3b are shown in Fig. 5b. For Im β < 0, ˙ W$_{he}$ is negative, thus cooling the effective temperature of the resonance mode and operating as a refrigerator. For 0 < Q$_{int}$I 2 $_{dc}$Im β < 1 heat is pumped into the resonance mode, increasing T$_{eff}$ until a new steady state is reached at which ˙ W$_{e}$ + ˙ U$_{i}$ = 0. In this mode the engine beam operates as a heat pump. In Supplementary Discussion A an analytical model is derived for the complex Young’s modulus and spring constant k ∗ eff and it is used to derive an estimate for β$_{0}$ , which yields β$_{0}$ = − 132 + 157 i A − $^{2}$. A finite element simulation of the full geometry including the anisotropy of the silicon crystal results in β$_{0}$ = − 111+58 i A − $^{2}$. Both the excellent fits in Fig. 3 and 4, and the quantitative agreement between the measured and simulated values of β$_{0}$ support the proposed feedback mechanism in Fig. 1d. FIG. 4: Measurements and fits to confirm the operation mechanism proposed in Fig. 1d. a,b Inverse effective Q -factor 1 /Q eff versus DC current I 2 dc for different values of the load impedances Z load indicated in the graph. Dashed lines are fits using equation (2). Symbols on the x -axis correspond to the minimum current level at which self-sustained oscillation was observed. a Resistive load impedances Z load . b Capacitive and inductive load impedances Z load . c Measured slopes − Im β of the curves in Fig. a , b plotted along the real (left) and the imaginary (right) γ Z -axis for resistive (red circle), capacitive (blue square) and inductive (green diamond) values of Z load . An excellent fit of the data is obtained by plotting − Im ( γ Z β 0 ) with β 0 = − 123 . 7 + 64 . 6 iA $^{$_{−}$2}$, for resistive (dashed lines), capacitive and inductive (solid lines) values of Z load . If Im β is positive the engine beam acts as a heat pump that increases the effective temperature and if Im β is negative the engine beam acts as a refrigerator that cools the effective temperature of the resonance mode. Flow rates of heat and work