Accurate modeling and characterization of photothermal forces in optomechanics

Photothermal effects have been pointed out as prominent sources of forces in optomechanical systems, competing with the standard radiation pressure interactions. In this Article, we derive a novel and accurate model for the prediction of photothermal forces and establish how some previous proposals can be complemented to yield precise results. As a proof-of-concept, we perform numerical and experimental tests on GaAs microdisks cavities and obtain striking agreement with our framework, revealing the importance of considering surface photothermal forces and the effects of multiple thermal modes in microphotonic devices.


INTRODUCTION
Microscale photonic devices revolutionized the field of cavity optomechanics over the past two decades. The ability to selectively control the photon-phonon interaction through the detuning between optical resonances and external laser sources led to novel applications, ranging from nonlinear dynamics [1,2] to the quantum manipulation of mechanical degrees of freedom [3][4][5]. The extreme confinement of the optical fields and small effective masses result in devices with enhanced optomechanical effects, which may arise from distinct and competing mechanisms, such as photothermal (bolometric) forces [6], radiation pressure [7], and piezoelectricity [8]. While the last two are built upon a robust theoretical framework, based on optical and mechanical modal analysis, photothermal forces in optomechanical systems are often treated with phenomenological models that require a complete experimental characterization of the structures as input [6,9,10]. The absence of an accurate and predictive photothermal force model has hindered its understanding and control within photonic devices.
Despite the high optical quality factors (Q > 10 4 ) of typical optomechanical resonators, absorptive losses may significantly impact the dynamics of mechanical modes [11,12]. As depicted in Fig. 1 a), Brownian noise-driven mechanical motion modulates the number of circulating photons in the cavity, which, in association with absorption, drives oscillations in the temperature of the system. Finally, thermallyinduced stresses couple back to the mechanical domain and close a feedback loop that defines the so-called photothermal backaction [13,14]. The finite thermal and optical response times yield forces that are time-delayed relative to mechanical oscillations, allowing for the cooling [15,16] and amplification [9,10] of mechanical normal modes in a range of dielectric and plasmonic [17] resonators. This extra degree of freedom opens up new possibilities, such as thermallymediated optomechanical ground state cooling in the badcavity regime [14].
In this work, we propose and demonstrate a model for photothermal (PTh) forces. We introduce a novel mathematical treatment for the description of the thermal fields that ultimately allows the prediction of the PTh response in devices with arbitrary geometries. It is built upon thermal modal analysis [18] and perturbation theory under a linear diffusive heat transfer regime, overcoming the limitations of several previous models based on a phenomenological treatment of PTh effects. As an example, we perform experiments in a GaAs microdisk cavity that remarkably agree with our predictions.

MODEL
The mechanical system is described by the equation of motion: ρ¨ U = ∇ · σ, where U denotes the displacement field, and σ the stress tensor. In thermo-elasticity, the selfconsistency of this problem requires a constitutive relation linking the stress tensor to the displacement and thermal fields. Since the stress arises solely from elastic deformations [19], it is necessary to split the strain of the system, S = 1 2 (∇ U + ∇ T U) =∇ U, into elastic (S x ) and thermal (S θ ) components as S = S θ + S x . The constitutive relation then reads: σ = c:S x = c:S − c:S θ [20], where c is the stiffness tensor, and ":" denotes the tensor contraction operation. Due to the thermal strain contribution, the free-boundary condition (σ ·n = 0), commonly used in micromechanical devices, leads to a temperature-dependent surface-traction on S that acts as a drive for the mechanical fields, as detailed in section S1 of the Supplemental Material. In some previous formulations of PTh forces in optomechanical systems, this subtlety has been neglected [21] and can lead to inaccurate predictions of the PTh response of optomechanical resonators.
A simple way to account for the boundary conditions is by calculating the PTh force directly from the work done by thermally-induced stresses on a given mechanical mode [18,22,23]. This is done in a linear approximation, where the elastic strain can be decomposed in the mechanical normal modes of the system ( u n ) as S x = ∑ n x n (t)∇ u n ( r) = ∑ n x n (t)S x n ( r), uncoupling the photothermal forces acting on each of the modes. The x n (t) are the normal mode amplitudes and are analogous to the generalized coordinates in analytical me-chanics [24]. The present calculation allows, to first order, direct access to an expression for the lumped PTh force on a mechanical mode [25][26][27][28]: F θ n (t) = ∂ ∂ x n S x :(c:S θ )dV = S x n :(c:S θ )dV.
A similar reasoning and expression was successfully used to estimate thermoelastic damping in MEMS [18], and displayed exceptional agreement with experimental values. Note that a partial differentiation in the amplitude x n (t) is performed, where the index n denotes the mechanical mode in which we are evaluating the PTh force. The first integral in the above expressions resembles the known (elastic) strain energy [29]. This association allows us to interpret it as the energy transferred between thermal and mechanical domains. The lumped photothermal force in Eq. 1 can be rewritten as the sum of a volume and a surface contributions, F θ n = F θ -Vol. n + F θ -Sur. n , where: F θ -Sur. n = u n · (c:S θ ) · d S, from which we verify that the PTh force field is composed of surface ( t θ ) and body ( f θ ) loads given by: t θ ( r,t) = (c:S θ ) ·n.
The volume load f θ was used in past work on optomechanical PTh forces to provide an estimate of their magnitude [21]. We demonstrate here that both surface and volume contributions are generally relevant in microscale devices, and must be considered for accurately describing dynamical backaction in optomechanical systems.
In order to grasp the time-dependence of PTh forces, a constitutive relation between the thermal strain and temperature field δ T ( r,t) must be assumed, S θ = αδ T ( r,t), where α is the thermal expansion tensor. The temporal analysis can be simplified by expanding δ T ( r,t) in multiple thermal modes [30], δT k ( r), with different relaxation constants, τ k , as: δ T ( r,t) = ∑ k θ k (t)δT k ( r), where θ k (t) is the k-thermal mode amplitude. This procedure is described in detail in section S2 of the Supplemental Material. In this framework, the PTh force for the n-th-mechanical can be written as a sum of the contribution from multiple thermal modes as F θ n = ∑ k Λ θ k,n θ k , where Λ θ k,n = S x n : (c : α) δ T k dV . Similarly, the surface and volume contributions can be decomposed in terms of Λ θ -Vol. k,n and Λ θ -Sur. k,n . Due to the optical drive (at frequency ω l ), the evolution of the thermal amplitudes θ k (t) is given by: where the thermal response to the optical heat source is modeled through the thermal relaxation time τ k , the thermal resistance R θ k and the optical absorption rate κ abs . The electromagnetic mode amplitude a(t) is normalized such that |a(t)| 2 is  Lagevin-type force F L induces fluctuations, δ x, in the mechanical position. As a consequence, the optical resonance frequency is shifted by δ ω causing a delayed modulation in the number of circulating photons δ n. If optomechanical backaction is considered, optics and mechanics are coupled through a force F RP , caused by radiation pressure and electrostriction. On the other hand, PTh backaction is obtained if absorption occurs, causing fluctuations in temperature δ θ . Thermal stresses modeled by the PTh force F θ act back on the mechanics inducing cooling or amplification (±) of mechanical modes. b) Geometric parameters for microdisk resonators, along with optical mode and temperature profile induced by absorption. c) Thermally induced deformations on disk structure, calculate through a fully-coupled thermoelastic model, volume PTh, and full PTh forces. d) Displacement field on ther andẑ directions as a function of the radial coordinate. The projections were evaluated at half the thickness (t/2) of the GaAs membrane, depicted by the red dashed curve in b).
the number of circulating photons in the cavity. In this analysis, we neglect thermoelastic damping heat relative to the optical absorption one. Although equations similar to Eq. 6 appeared in the literature [6,12], their use was performed considering a single "effective" thermal relaxation time obtained from a fit of experimental data regarding the magnitude of the PTh forces or from a frequency-resolved measurement of the resonator's thermo-optical response. Those two methods often yield contrasting results [31]; our model indicates that this comes in place due to different thermal modes being relevant to the thermo-mechanical and thermo-optical processes. For instance, in our devices, the thermo-optical response is well represented by the fundamental thermal mode while this approximation for PTh effects yields approximately 200% error in the overall force. We stress that in the case of "effective" thermal relaxation times no formal derivation of Eq. 6 can be given, and its use is only justified from a phenomenological point of view.

NUMERICAL MODELING
Static case -In order to numerically validate the present model, we first consider the case of static thermal deformations on a GaAs on Al 0.7 Ga 0.3 As (250 /2000 nm) microdisk with radius R = 6 µm and pedestal radius R ped = 0.75 µm. The first radial order optical TE mode of the disk is used as a heat source that drives a stationary temperature field in the structure, as shown in Fig. 1 b). Due to the static nature of this problem, thermal modal analysis is not necessary, such that the full thermal field is used in all calculations following.
We use finite element method (FEM) calculations to compare the thermal displacement u θ predicted by the derived PTh force field (PTh) to a fully-coupled thermoelastic model in COMSOL Multiphysics © . The fully-coupled model calculations are carried out in the linear elastic approximation, in consistency with the hypothesis used in our derivation. For completeness, we further calculate the temperature-induced displacement resulting from the volume force (PTh -Vol.) alone, as shown in Fig. 1 c); the thermal displacement field components alongr andẑ directions are shown in Fig. 1d). Our PTh force formulation, which includes both surface and volume contributions, accurately reproduces the fully-coupled model thermal displacement field, with major deformations present near the edge of the disk. This is in stark contrast with the volume-only PTh force calculations, where deformations are mostly confined to the pedestal region. This discrepancy indicates that thermo-mechanical coupling calculations can be critically affected by the existence of the surface PTh force in microphotonic structures.
Dynamic case -To study how surface and volume loads may affect an optomechanical system, we solve the coupledmode equations for optical, thermal and mechanical amplitudes in the frequency domain (S3 of the Supplemental Material), under a continuous-wave optical excitation which allows the linearization of these equations around an average amplitude. This procedure allows one to evaluate and compare the average photothermal force per photon obtained from surface and volume contributions. Neglecting optical resonance frequency shifts due to temperature variations, the lumped bolometric force (PTh) per photon is given by: where χ θ k (Ω n ) is the thermal susceptibility of the k-th thermal mode evaluated at the angular frequency of the n-th mechanical mode.
The surface (PTh -Surface) and volume (PTh -Vol.) contributions are shown in Fig. 2 a) for three different mechanical modes. Here, we include the effects of fluctuations in the optical frequency due to the thermo-optical effect. Calculations involving the volume force alone not only would drastically overestimate the total PTh forces but would also carry a flipped sign with respect to the correct results. In our case, localized losses such as surface absorption do not modify the PTh response appreciably -as shown in S4 of the Supplemental Material -, hence we assume that the absorptivity of GaAs is homogeneous and isotropic throughout the device. This result holds if the PTh response is dominated by low-order thermal modes that are essentially homogeneous in the region of confinement of the optical mode. If this is not the case, a thorough characterization of the nature of absorptive losses is necessary to obtain accurate predictions from the model.
Surface and volume contributions are obtained by replacing Λ θ k,n with Λ θ −Sur.
k,n in Eq. 7. Since the χ θ k (Ω n ) are complex numbers, forces are composed of real and imaginary parts; the latter is largely dominant in the total PTh forces. Physically, this phenomenon is related to the relatively large thermal relaxation times 1/τ k Ω n of the relevant thermal modes, which cause their response to lag behind the mechanical oscillations. The optical absorption rate was chosen to be κ abs /(2π) = 1 GHz following state-of-the-art experiments on GaAs microdisks [32]. The total loss rate (κ = κ e + κ abs + κ non−abs ) is κ/(2π) ≈ 1.93 GHz, with extrinsic coupling rate (i.e. coupling to a waveguide) κ e /(2π) ≈ 0.48 GHz. These numbers yield a loaded quality factor Q opt ≈ 10 5 ; all other parameters are obtained through FEM simulations, where first order TE optical mode was considered.
The summation in Eq. 7 raises a question on the number of thermal modes that must be accounted for to correctly evaluate F θ n . In Fig. 2 b) we consider the 230 MHz mechanical breathing mode, and calculate the contributions of the surface and volume components to the total PTh force per photon as a function of the number of thermal modes considered, ordered by decreasing τ k . The total PTh response -largely dominated by the imaginary component -takes only 6 thermal modes to converge reasonably, whereas the volume and surface forces take ≈ 40 thermal modes. This feature arises from the fact that in high order (k > 6) thermal modes volume and surface terms yield opposite contributions that approximately cancel each other. Physically, this comes in place because the temperature profiles of high order thermal modes are associated with rapid spatial oscillations within the microdisk; as the spatial frequency increases, strong temperature-gradient forces (Eq. 2) are generated. This behavior, however, is not verified for the total force in Eq. 1, in which the gradient operation appears acting on the elastic deformations through the quantity S x . Consequently, as the order of the thermal modes is increased, the integrand in Eq. 2 overcomes its counterpart in Eq. 1, indicating that volume forces should become larger than the total force (in absolute terms). The only way such phenomenon can be observed is if surface and volume contributions counterbalance each other. We illustrate this in Fig. 2c), where the profiles of the total and volume pressures per photon are displayed for three different thermal modes δ T k , k = 1, 3, 50. Remarkably, in the k = 1, 3 modes, the total and volume photothermal pressures are of similar amplitude, whereas for k = 50, the volume component is almost two orders of magnitude larger than its counterpart. Lastly, since volume and surface terms yield opposite contributions even for the dominant low order thermal modes, surface engineering may emerge as a route for the enhancement or even cancellation of PTh forces. Radiation pressure (RP) and photothermal (PTh) backaction-induced mechanical a1) frequency and a2) linewidth shifts for the mechanical breathing mode in Fig. 2. b) Ratio between mechanical damping rate modifications δ Γ θ /δ Γ RP induced by PTh and RP forces evaluated at half optical linewidth (∆ = κ/2). The star marker depicts the device in a1) and a2).
We now turn our attention to the complete optomechanical interaction considered in this work, composed of radiation pressure and photothermal forces. Both contribute independently to an effective optomechanical backaction and must be considered for a correct description of the effects that will be studied in our experiment. We consider the same device as in Fig. 2, with a mechanical breathing mode at 230 MHz and 50 µW incident power. In Figs. 3 a1) and a2) the PTh and RP backaction curves are displayed. For the RP calculations, both photoelastic [33] and moving boundary contributions [34] are considered. While RP dominates the optically-induced frequency shift, cooling and amplification are largely dominated by PTh forces. This is due to slow thermal responses (when compared to the mechanical periods) yielding PTh forces outof-phase with respect to the mechanical oscillations, which favors mechanical linewidth modification processes. This is a key feature that is explored in our experiments. Importantly, for GaAs microdisks, PTh and RP effects add constructively in cooling/amplification processes. In Fig. 3 b) the ratio of PTh and RP cooling at ∆ = 0.5κ is evaluated as a function of κ abs and the disk radius; for these calculations, the pedestal radius R ped = 0.75 µm is kept fixed. Such diagram can be used as a tool for choosing geometries in order to maximize or suppress PTh effects: while larger disks display PTh-dominated dynamical backaction (red region), in smaller disks -where optical and mechanical modes are more tightly confined and with larger overlap -RP interaction prevails (blue region). The marker displays the parameters used in Figs. 3 a1), a2).

EXPERIMENTAL RESULTS
The effectiveness of the thermodynamic description is tested by monitoring the modification on the mechanical linewidth of a cavity optomechanical system consisting of a R ≈ (5.1 ± 0.1)µm, R ped ≈ (0.6 ± 0.1)µm GaAs/Al 0.7 Ga 0.3 As (250 /2000 nm) microdisk. The experimental setup used to characterize both the mechanical and optical spectra is shown in Fig. 4 a). A scanning electron microscope (SEM) image of the fabricated device is also presented. Light emitted by a tunable laser source is coupled in and out the resonator through a tapered fiber loop. The output from the cavity is collected at both fast and slow photodetectors. The fast response is fed into an electrical spectrum analyzer (ESA), while the slow signal is collected by an analog-to-digital converter (DAQ). A Mach-Zehnder interferometer (MZI) and hydrogen cyanide reference gas cell (HCN) are used for the calibration of the cavity's optical response. A detailed description of the experimental setup is found in S5 of the Supplemental Material.
A thorough optical characterization of the device is necessary in order to calibrate both nonlinear losses and thermal frequency shift, both crucial to accurately predict backaction effects at high incident powers. We monitor the optical transmission spectrum of the cavity for various incident powers, as illustrated in Fig. 4 b). The cold-cavity transmission yields intrinsic and extrinsic optical damping rates of κ i /(2π) = 7.0 GHz and κ e /(2π) = 4.2 GHz. Assuming that the coupling to the fiber taper remains constant during the measurements, power-dependent changes in the transmission can be traced back to recover the nonlinear losses and internal optical energy of the resonator (U) [35,36]. Nonlinear frequency shift is directly obtained through joint calibration with the MZI and tracking of the resonance shift (∆ω). Nonlinear losses and frequency shift data are then simultaneously adjusted to polynomial curves to obtain absorptive and nonabsorptive optical dissipation rates as a function of the energy in the resonator. The polynomial approximation is valid for sufficiently low input powers, where only terms up to O(U 3 ) in internal energy are enough to describe our results. A comprehensive guide for this analysis is found in S6 of the Supplemental Material. In Fig. 4 c) we show ∆ω as a function of the internal energy in the resonator. The inset shows the dispersion for low incident powers, critical for determining the portion of the cold-cavity losses with absorptive nature. Fig. 4 d) displays the total optical dissipation rate of the system split in absorptive (κ abs ) and non-absorptive (κ non-abs ) parts.
We measure the back-action effects by monitoring the mechanical mode spectrum through the RF power spectrum, which is recorded for a range of positive (blue) laser-cavity detuning, resulting in spectrograms similar to Fig. 4 e). From a Lorentzian fit (shown in Fig. 4 f)) both mechanical frequency (δ Ω) and linewidth (δ Γ) changes are obtained and the latter is compared with the predictions of the PTh and RP models previously discussed. For the tested device, RP yields negligible contribution, demonstrating the role of distinct backaction mechanisms in explaining the observed phenomena. The optical mode excited in our measurements is identified through its free-spectral range (FSR), consistent with the 6-th order TE optical mode (S7A of the Supplemental Material).
The estimated photothermal response of the system is obtained through a combination of FEM simulations for a mechanically anisotropic GaAs microdisk (S7B of the Supplemental Material) and the experimental nonlinear loss and dispersion described above. Importantly, our FEM resultswhere the parameters R θ k and Λ θ k were evaluated -require only thermal and thermo-elastic material properties as input, all of which are well established in the literature [37]. Fig. 4 g) exhibits the comparison between the measured maximal δ Γ and its theoretical estimate as a function of the incident powers on the cavity. Those values are obtained from measurements of δ Γ as a function of the laser to cold-cavity detuning, ∆ , (i.e. ∆ = 0 refers to the cold-cavity resonance frequency), exemplified in Figs. 4 h1)-h4) for four different input powers. Details on the collection and analysis of data regarding the mechanical response are found in S8 of the Supplemental Material, along with the treatment of the stiffening of the mechanical oscillator (δ Ω), which is dominated by a static temperature softening of GaAs [38], red-shifting the mechanical frequency up to −20 kHz. As expected, the uncertainty in δ Γ becomes more relevant at low input powers, since in that case the optomechanical transduction is less efficient.
Excellent agreement between our prediction (curves) and experiment (markers) is found over the whole range of detuning measured in Figs. 4 h1)-h4). Despite the encouraging nature of our results, we stress that, as discussed in S4 of the supplementary material, in the microdisk geometry it is hard to distinguish bulk from surface absorption. Experiments with more sophisticated resonators, in which PTh effects are sensitive to the spatial distribution of the absorption, provide an interesting route for testing further our predictions.

CONCLUSION
In summary, we have proposed and verified experimentally a model for the photothermal forces acting on cavity optomechanical systems derived through thermodynamic considerations. The theoretical estimates were shown to display remarkable agreement with our measurements, thus providing a solid route to design the thermo-optomechanical response in nanomechanical resonators. Also, the modal treatment for the thermal response is a significant step towards thermal engineering in the broad field of nanophotonics, and paves the way for a new class of experiments where those effects are tailored to interest. Finally, although GaAs based devices were taken as an example, we emphasize photothermal forces can be appreciable in other platforms and geometries and that the content of can provide insight in those cases.

ACKNOWLEDGMENTS
The authors would like to acknowledge CCSNano-UNICAMP for providing the micro-fabrication infrastructure and CMC Microsystems for providing access to MBE epitaxy and the GaAs wafers.

DATA AVAILABILITY
FEM and scripts files for generating each figure are available at Ref. [39]. Additional data that support the findings of this study are available from the corresponding author upon   Tab.I summarizes all experimental and simulated parameters used throughout the main and supplemental text. Each parameter used on our model is also presented with its value, source (the way they were obtained), and the section where it is discussed.

S1. DERIVATION OF THE PHOTOTHERMAL COUPLING
We discuss two different approaches to derive the photothermal coupling; the first is built from thermodynamic arguments and the second directly from the mechanical equation of motion. Our analysis assumes that all transport processes are quasi-static, in the sense that local equilibrium is assumed at each step. This hypothesis allows us to safely use definitions of quantities such as temperature and is valid within the mechanical frequencies considered in this work. Local heat transport is ultimately related to thermal phonons relaxation times which range around 10 −11 s, much faster than the mechanical oscillations with periods around 10 −9 s.

A. Thermodynamic derivation
The total strain tensor (S) in a thermoelastic body is given by: where S θ and S x are the thermal and elastic components of the strain, respectively. The symbol∇ = 1 2 (∇ + ∇ T ) denotes the symmetric gradient tensor operation and U is the displacement field. The thermal strain tensor is defined as: where α is the rank 2 thermal expansion tensor and δT is the temperature field of the body. In an isotropic medium α can be simplified to, αI, where α is the thermal expansion coefficient and I is the 3 × 3 identity matrix. Denoting c as the stiffness tensor, the stress tensor is obtained through the following constitutive relation: We argue that the additional stresses related to photothermal forces are necessarily related to the term − c : S θ . Evaluating the work per unit volume δw done during an arbitrary deformation δ U we get: where X j (j = 1, 2, 3) denotes each of the cartesian coordinates and summation over repeated indexes is used. Considering the work in the entire volume where the deformation takes place we have: Using Gauss's theorem and integration by parts: For a free surface (σ ·n = 0), the second integral on the RHS vanishes and we get: where we have identified ∂Ui ∂Xj with the total strain. We are only interested in forces that effectively transfer energy from the thermal to the mechanical domain over a stationary cycle [5], specifically, in the form of work. Substituting Eq. S3 in Eq. S7 we get: Integrating this expression over a mechanical cycle and noting that the first term on the right-hand side is an exact differential -and therefore vanishes -, we conclude that c ijkl S θ kl dS ij is responsible for feeding energy into the mechanical domain. Aiming at finding a linear theory, underpinned by modal analysis, we keep only terms to first order in temperature. The final expression for the work per unit volume, done by photothermal forces, is given by: Integrating over the entire volume we obtain the total photothermal work: The elastic strain can be decomposed in the mechanical modes u n ( r) of the system as S x = n x n (t)S x n , where x n (t) is the n-th mode amplitude (i.e. the u n ( r) are normalized such that its global maximum is set to be 1) and S x n =∇u n ( r) . Using the modal amplitudes x n as generalized coordinates in the sense of analytical mechanics, the generalized force on the n-th acoustic mode is: Or simplifying: and explicitly in terms of mechanical modes and the temperature field: (S13)

B. Derivation from the mechanical equation of motion
We start with the mechanical equation of motion: with the boundary condition: In view of Eq. S3, these equations can be recast as: with the boundary condition: c:S ·n = c:S θ ·n. (S17) The thermal strain in the above equations can be pictured to induce both an external volume load f θ = −∇ · c:S θ , and a surface traction t θ = c:S θ ·n.
In order to perform the modal decomposition of the photothermal force it is important to consider the mechanical eigenmodes of our device. These are solutions of: with free (untractioned) boundary conditions: where the strain is given by S x n =∇ u n , and Ω n is the angular mechanical frequency of the n-th mechanical mode. Projecting Eq. S16 in a mechanical mode u n and performing a volume integration, we get: ρ u n ·¨ U dV = u n · ∇ · c:S dV − u n · ∇ · c:S θ dV.
Integration by parts of the first term on the right-hand side yields: u n · ∇ · c:S dV = u n · c:S ·n dS − U · c:S x n ·n dS + U · ∇ · c:S x n dV.
Using Eqs. S17 and S19 to rewrite the surface integrals and also using Eq. S18 to recast the volume integral on the right-hand side, Eq. S21 reads: Substituting Eq. S22 back in Eq. S20: Finally, since the first term in the RHS is only defined inside the device volume, we can then perform a modal expansion of the displacement field as U = n x n (t) u n ( r). Using the orthogonality relations for the mechanical modes, ρ u n · u m dV = m eff,n δ n,m , where m eff is effective mass, we get: m eff,n ẍ n + Ω 2 n x n = u n · c:S θ ·n dS − u n · ∇ · c:S θ dV, Usually, a damping term Γ nẋn is heuristically introduced on the left hand side, but will be omitted for the moment. The surface and the volume contributions to the photothermal force can be assembled in a single volume integral: Which agrees with the thermodynamic treatment previously discussed. In previous works the surface load, which is evident in our analysis, was ignored.

S2. THERMAL MODE ANALYSIS
The heat equation for diffusive heat transport reads: where c p , ρ and k th are respectively the specific heat, mass density and thermal conductivity of the medium.Q is the heat generation in the structure, which will be associated with optical absorption.
In the absence of heat sources, this equation is separable and admits solutions of the form δT ( r, t) = δT ( r)e −t/τ leading to a generalized eigenvalue equation for the thermal modes: Since the operators ∇ · k th ∇) and c p ρ are both Hermitian and the latter is positive definite, the equation above is associated with a complete set of eigenmodes, δT ( r, t) = k θ k (t)δT k ( r). Applying this expansion to Eq. S26 yields: Assuming fixed temperature or insulating boundary conditions, one derives orthogonality relations for the thermal eigenmodes: which allows the time evolution of any thermal mode to be written as: Applying the thermal mode expansion over the expression for the lumped photothermal force, Eq. S13, allows us to write F θ n (t) = k F θ k,n (t), where: where the coefficient Λ θ k,n has units of N/K. We are interested in heating due to optical absorption. Assuming that only one optical mode is excited and denoting the local absorptive (optical) loss rate as κ loc ( r, t), the total power dissipated per unit volume is given by: where |a(t)| 2 is the number of photons circulating in the optical cavity, e( r) is the electromagnetic mode profile, ε 0 and ε r are the vacuum permittivity and dielectric function, respectively. We would like to express Eq. S30 in terms of the experimentally accessible total absorptive loss rate κ abs (t) (the time dependency arises if nonlinear absorption is considered) instead of κ loc ( r, t). These quantities are related by: Since κ abs (t) is part of a lumped model, it is clear that with its knowledge alone, one cannot solve for κ loc ( r, t). Instead, if linear absorption is considered, we simplify matters by assuming κ loc ( r, t) is a constant κ loc,0 within an absorbing volume V abs . In that case, we may solve for κ loc,0 in terms of the total linear absorption κ abs = η 0 κ 0 , where κ 0 denotes the cold-cavity (linear) intrinsic losses of the device, and η 0 its absorptive fraction. We have: Although assuming homogeneous absorption within the absorbing volume is certainly an approximation, it is generally sufficient to correctly capture optical absorption even when it takes place near the surface. For example, in case the main mechanism of absorption is mediated by surface imperfections, we should expect an exponential decay in κ loc ( r) as we move away from the surface. In that case, we could simply treat κ loc and V abs as effective quantities.
In order to incorporate nonlinear absorption mechanisms such as two-photon (TPA) and free-carrier (FCA) absorption, the spatio-temporal dependency of κ loc ( r, t) must be related to the fields in the resonator. For phenomena that scale linearly with the electromagnetic energy in the resonator, e.g. TPA [6], we expect κ loc ( r, t) to have the following form: since the fields themselves generate absorption. Analogously, for phenomena that scale quadratically with the electromagnetic energy in the resonator, e.g FCA: where nonlocal effects such as carrier diffusion are neglected.
In both cases above, C 2 and C 1 are approximated as material-dependent constants which are only appreciable inside the volume of the cavity, defined by V cav . They can be expressed in terms of the experimentally accessible quantities such as the first (κ abs,2 (t)) and second (κ abs,2 (t)) order absorption rates: Writing κ loc ( r, t) = κ loc,0 + κ loc,1 ( r, t) + κ loc,2 ( r, t) in Eq. S32 and using the definitions of C 1 and C 2 above, equation S30 can be written in terms of the photon energy ω l and the zeroth, first, and second order thermal resistances, yielding:θ with: Eqs. S40, S41, S42 can be evaluated in a straightforward way by any electromagnetic/thermal solver. Mechanical vibrations may also lead to heating through thermoelastic damping, however, due to its negligible dynamical contributions in the devices discussed here, this effect will be neglected. Importantly, if δT k ( r) is approximately constant in the region of overlap with the optical mode, all three expressions above yield similar values. This approximation is well justified for the first few low-order thermal modes, which turn out to be the most relevant for the case treated in the main text.
In practice, not all of the nonlinear absorptive losses (κ 1 and κ 2 ) are necessarily converted to heat, as implied by Eq. S39. Nonlinear absorption is strongly tied to carrier dynamics inside the device. As we are dealing with a direct band-gap semiconductor, we expect that carriers generated through virtual processes such as TPA to display both radiative and non-radiative recombination routes. That is, only a fraction of the nonlinear loss rate acts as a heat source in our resonator. As a consequence, we are led to define the ratios η i of the ith order nonlinear loss rate that contributes to the heating of the resonator. With that in mind, the generalized evolution of the amplitude θ k (t) is given by:θ

S3. FREQUENCY DOMAIN DESCRIPTION OF THE PHOTOTHERMAL COUPLING
A. Single-mode Photothermal coupling We start from the coupled equations for the temporal evolution of the optical, mechanical and temperature fields. In this subsection, we assume that the dynamics is dominated by one mode in each of the previously cited domains, therefore, we drop all mode-related indices. For simplicity, in this and in the next subsection we treat the case of linear absorption, i.e. κ abs,1 (t) = κ abs,2 (t) = 0. The equations read: where ∆(t) = ω l − ω c + G θ θ + G x x, is the relative detuning between the cavity and the laser, including optomechanical and thermal frequency pullings, ω c is the (optical) resonance frequency and ω l the external laser frequency. F L , F RP and F θ are respectively the thermal Langevin, the radiation pressure and the photothermal (bolometric) forces, G x and G θ are respectively the optomechanical and thermo-optical coupling rates. We also choose a normalization for the electromagnetic modes such that |a| 2 is the number of circulating photons in the cavity, consequently, α in is the rate of incident photons on the resonator, driven by an external source. We may derive expressions for G θ and G x using first-order perturbation theory on Maxwell's equations [4,7]: where Eq. S45a is related to the thermo-optic effect, related to temperature induced changes in the dielectric response of the medium. The integrals in Eq. S45b model the photoelastic and moving boundary optomechanical couplings related to a mechanical mode u, where p is the photoelastic tensor. The moving boundary contribution is a function of ∆ε = ε 1 − ε 2 and ∆ε −1 = (ε 1 ) −1 − (ε 2 ) −1 , which are related to the permittivities of the guiding (ε 1 ) and surrounding (ε 2 ) materials. In the spirit of deriving an entirely linear theory, we look for linearized solutions to Eqs. S44 as a(t) = a 0 + δa(t), x(t) = x 0 + δx(t), and θ(t) = θ 0 + δθ(t). Using F RP = G x |a| 2 , F θ = Λ θ θ and collecting the fluctuating terms up to first order we obtain the following dynamical equations for thermo-optomechanics (initially disregarding the Langevin term): where ∆ 0 absorbs the static frequency shifts given by G θ θ 0 + G x x 0 . In frequency space: since [δa] * (ω) = [δa(−ω)] * and for real variables such as δθ and δx, δx(ω) = [δx(−ω)] * , we may eliminate δa(ω) and [δa] * (ω), obtaining: ( where Ψ(ω, ∆ 0 ) = [ 1 (∆0+ω)+iκ/2 + 1 (∆0−ω)−iκ/2 ], a combination of the cavity's optical susceptibility and its conjugate, and χ θ (ω) = 1 1/τ −iω is the thermal susceptibility. Eliminating δθ(ω) in Eq. S48b, we readily obtain the dressed optomechanical self-energy. The term on the RHS which is being multiplied by G θ yields the correction acting on the radiation pressure feedback, therefore the total RP modification to the linear response of the system is: where we defined Σ RP (ω) = (G x ) 2 Ψ(ω, ∆ 0 )|a 0 | 2 , the optically induced mechanical inverse susceptibility in the absence of thermal feedback. The correction above is negligible for the GaAs microdisks considered in the main text, however it might become relevant in other systems where mechanical frequencies are comparable to or smaller than the inverse of thermal relaxation times. In that case, thermal dispersion affects the optical spectrum within the mechanical timescale, impacting optomechanical transduction. The term proportional to Λ θ in Eq. S48b yields the dressed bolometric force, which is this work's object of study. The photothermal contribution to the inverse mechanical susceptibility is given by: where we defined the bare bolometric contribution to the inverse mechanical susceptibility as: The dressed inverse mechanical susceptibility χ −1 m,eff (ω), including RP and photothermal contributions is finally given by: where χ −1 m (ω) = m eff Ω 2 − ω 2 − iΓω is the inverse of the bare mechanical susceptibility. The expressions above for Σ RP eff and Σ θ eff may be simplified in many practical scenarios. In agreement with the experimental case reported in the main text, we assume the bad-cavity regime (Ω κ/2) and τ 1/Ω. Evaluating all frequency responses at Ω , which is the driving frequency given by the mechanical oscillations, and also taking the detuning ∆ 0 = κ/2, we note that χ θ ≈ i/Ω, and Ψ ≈ 2/κ. For typical microdisks, τ /(2π) ≈ 10 µs, Ω/(2π) ≈ 250 MHz, κ abs /(2π) ≈ 1 GHz, ω l /(2π) ≈ 200 THz, R θ ≈ 4 × 10 4 , G θ /(2π) ≈ 10 GHz, κ/(2π) ≈ 10 GHz and |a 0 | 2 ≈ 10 4 . With these numbers we see that: Confirming that the correction due to the thermo-optical dispersion is negligible in our case. One last approximation can be made in order to compare Σ θ and Σ RP . Considering only cooling/amplification under the same approximations as before, the ratio of photothermal (δΓ θ ) and RP (δΓ RP ) mechanical amplification, evaluated at ∆ 0 = κ/2 is given by: which can be used as a simple test for the relevance of photothermal forces in devices that operate within the regimes above.

B. Multimode Photothermal coupling
We now generalize our discussion to the case where several thermal modes couple relevantly to optical and mechanical modes. The appropriate equations in frequency space are given by: Eliminating δa(ω) from the equations for the thermal modes we get: where we defined: By inspection, one may show that δθ l (ω) ∆0) . With that result we get: and the total photothermal modification to the inverse mechanical susceptibility is: where the bare photothermal contribution is given by: In analogy with the single thermal mode case, the corrected RP induced inverse susceptibility is given by: From the symmetries of the equations above, it is useful to define "effective thermal response" functions, at a given frequency ω: in terms of which, expressions S59 and S61 can be recast as: which share the same form as their single mode counterparts. We note that the summations for h θ 1 (ω) and h θ 2 (ω) are well behaved since R θ k /τ k is independent of τ k and χ θ k is a monotonically decreasing function of τ k , guaranteeing gradually decreasing importance of high-order thermal modes (τ k → 0). From the equations above it is straightforward to compare the photothermal and radiation pressure forces amplitudes. The ratio Σ θ /Σ RP yields: Since the denominator of the RHS essentially gives the RP force per photon, the numerator may be interpreted likewise, but regarding the photothermal force. Notice however, that due to non-zero relaxation times, the two forces display a non-trivial phase relation, here captured by the complex part of h θ 1 (ω). The approximation obtained in the previous section can be generalized to the multimode case. Again, considering the bad-cavity limit and neglecting thermal dispersion, i.e. h θ 2 (Ω) 1:

C. Nonlinear Multimode Photothermal coupling
In this subsection, we summarize the main results concerning the usage of Eq. S55, without neglecting the nonlinear losses in general. That is, we use Eq. S39 and make κ → κ(|a| 2 ) in the evolution of a. That is, the total nonlinear intrinsic loss, κ i (t), is given by κ i (t) = κ 0 + κ 1 (t) + κ 2 (t). The coupled equations of the optical, mechanical and thermal responses read:ȧ where A fundamental difference between the present case and the analysis in past sections is that, when considering the frequency domain version of the above equations, the Fourier transform of the nonlinear losses must also be performed, since those also display a time dependency, substantially modifying the final results.
In our model for the nonlinear losses we neglected all but the dynamics of the energy in the resonator. That is, the only time dependency of the first and second order dissipation is given through the number of photons |a(t)| 2 . Since those losses are intrinsically related to carrier dynamics, this is explicitly an approximation. Its validity requires that any timescale associated to carrier relaxation (i.e. its recombination rate) is must faster than the period of mechanical oscillations, in such a way that we may assume that their population responds instantly to energy variations. As stated in past work on this subject, typical recombination rates for GaAs range around 10 ps [8], much faster than the ≈ 4 ns associated with mechanical oscillations.
Under the aforementioned hypothesis, we write κ 1 (t) = κ 1 |a(t)| 2 and κ 2 (t) = κ 2 |a(t)| 4 , where the κ i are constants. This form allows us to linearize the losses according to a(t) → a 0 + δa(t). After linearizing and writing Eqs. S66 in frequency space, they read: where we defined the average lossesκ j = κ j |a 0 | 2j andκ 0 = κ 0 We follow the same procedure as before: we eliminate θ(ω) and δa(ω) from the equation for the mechanical response, from which we recover a dressed susceptibility for the mechanical mode. The PT-induced and the RPinduced modifications to the mechanical susceptibility are: where h θ 1,dyn (ω) and h θ 2,dyn (ω) are, respectively, the effective responses for the photothermal force and for the thermooptic dispersion: and Ψ dyn is the modified optical response due the dynamic of the nonlinear losses: The equations above were used in the analysis of the experiment reported in the main text.

S4. LOCALIZED ABSORPTION
In this section we use FEM simulations to investigate how surface localized absorption affects the photothermal effects on our microdisk devices. An absorptive layer with thickness t abs is defined within the GaAs membrane composing the microdisk illustrated in Fig. S1 a). Importantly, the absorption is taken to be constant over this layer. The radius of the disk is 6 µm and the pedestal diameter is 1.5 µm, the latter is measured in the region of contact of the AlGaAs/GaAs layers. The geometry of this device agrees with the one used in Fig. 2 in the main text. In the following analysis we consider the first order TE optical mode and the first order mechanical breathing (Ω/(2π) = 230 MHz) mode of the microdisk.
In the main text we study the case of an homogeneously absorbing GaAs membrane (no absorption layer), in which the total PTh force per mode tends to become negligibly small for higher order thermal modes (k > 5), this happens due to the self-cancellation of the surface and volume contributions. This is shown in Fig. S1 b) where each thermal mode's contribution to the imaginary and real parts total PTh force per photon was evaluated. As the thermal frequency -defined as τ −1 /(2π) -increases, the red T-shaped bars in the imaginary part of such force become negligibly small. Nevertheless, we verify appreciable -contrary -surface and volume forces for frequencies as high as 1 GHz, in agreement with the discussion in the main text. The total real component of the force (horizontal dashed line) is evaluated to be approximately two orders of magnitude smaller than its quadrature and its related red bars are barely seen when compared to the dark blue and cyan bars.
In Fig. S1 c1) and c2) we evaluate the consequences of PTh backaction on the mechanical mode's linewidth and frequency, respectively, for several different absorption layers, with t abs = 15 , 62.5 , 110 nm and also for the case of homogeneously distributed absorption in the GaAs membrane. The total surface and volume contributions are appreciably altered by the absorptive layer thickness, remarkably the overall response is essentially unaffected. This is understood under the light of the results of Fig. S1 b). Since only low-order thermal modes are important for the device's overall response and those are essentially homogeneous over the electromagnetic mode's confinement region, their response is unaffected due to the presence of an absorption layer. As the surface and volume forces are still appreciable for high-order modes -and those display rapid spatial oscillations -, their response is significantly modified by localized absorption. We conclude that the differences verified in the navy/cyan curves in Fig. S1 c1) and c2) are due to high-frequency thermal modes that barely contribute to the total response of the system.

S5. EXPERIMENTAL SETUP
Here we discuss the experimental setup used in our two experiments: the characterization of static optical nonlinearities (Section S7), and the measurement of optomechanical effects (Section S8).
Light emitted from a tunable laser source is sent into a (99/1) beam splitter. The 99% output is directed to the main measurement setup, while the remaining power is used in an auxiliary frequency calibration setup, shown in Fig. S2, composed of a Mach-Zehnder interferometer and a referenced HCN gas cell. In the main measurement setup, light passes through a phase modulator (PM), used in the calibration of the optomechanical coupling rate; a polarization controller (PC) and a variable attenuator (ATT1) used respectively to control the polarization and intensity of the light that will be coupled to the cavity.
Light is coupled in and out the microdisk resonator through a tapered fiber loop. The loop is supported by a carefully designed parking lot surrounding the cavity, allowing us to probe the microdisks response without ever touching it, which would degrade both optical and mechanical quality factors. The tapered fiber output is split between two photodetectors: a slow one, linked to the DAQ (Data AcQuisition system) and used to probe the device's static optical response; and a fast one, linked to the ESA (Electrical Spectrum Analyzer), for the measurement of high frequency fluctuations of the optical signal related to the brownian mechanical motion in the optomechanical resonator. The output variable attenuator, ATT2, placed right before the slow photodetector, is used to compensate attenuation changes given by ATT1, allowing us to maintain the electronic gain settings and incident power on the slow detector fixed, avoiding both issues of signal saturation and detection nonlinearities.
In the frequency calibration branch of our setup, light is evenly split between a 4 GHz FSR Mach-Zehnder interferometer (a relative frequency reference), and the HCN cell (an absolute frequency reference), see Fig. S2. The transmission spectra of both references are captured with the DAQ. Its features, such as the distance between peaks in the Mach-Zehnder and the extinctions of the HCN spectrum, are used as frequency references to construct a reliable frequency axis for the cavity's transmission spectrum.
The optical power is measured at the input and output of the tapered fiber loop, allowing us to characterize the propagation losses of the fiber taper, −3.86 dB, which are both due to the taper non-adiabaticity, and due to the scattering at the region coupled to the cavity. Considering the loss is approximately symmetric in the tapered fiber, that is, the losses between the coupling region and the input and output regions are approximately equal, we can estimate the incident power at the cavity region from the input power in the fiber.

S6. FITTING OF THE NONLINEAR OPTICAL RESPONSE
Our model predicts the photothermal effects given the absorption rate is known as a function of the cavity internal energy (number of circulating photons). Multiple linear and nonlinear processes contribute to absorption in typical semiconductor platforms, most of them are hard to determine theoretically, such as surface induced linear/two-photon absorption. Given such difficulty, we experimentally characterize the absorption rate. It is important to stress that we are not interested in carefully discriminating between the contributions of different nonlinear optical effects, but are rather focused on obtaining the heating intensity as a function of the cavity internal energy. In that sense, a phenomenological interpretation of our fitting is appropriate.
We perform measurements of the optical transmission spectrum as a function of the incident power on the cavity. For each incident power, we extract data points containing the optical frequency shift, the optical loss variation, and the cavity internal energy at resonance. Lastly, we concurrently fit such data points using a polynomial model that accounts for different kinds of optical nonlinearities.
All data and scripts discussed in this section can be accessed in Ref. [9].

A. Data collection
Once the fiber loop is coupled to the cavity and propagation losses are calibrated, the measurements are automated. The attenuation in ATT1 is varied from 28 dB to 8 dB in steps of 1 dB. For each attenuation, the laser frequency is swept 5 nm around the optical resonances, which are centered about 1565 nm; the transmission is then recorded and saved. The entire experiment took approximately 7 min, or 19 s per power level.

B. Extracting lumped parameters from the optical transmission spectra
In Fig. S3 a) we display the normalized transmission spectrum of our microdisk for the second lowest incident power, it is centered around the optical doublet used in all the experiments reported in the main text. The frequency axis is the same as in Fig. 4 in the main text, its definition is ∆ /(2π) = ω/(2π) − ω 0 /(2π), where ω is the absolute frequency and ω 0 = 191.5418 THz is approximately the frequency of the higher frequency local minimum of the transmission for the low incident power data, as seen in Fig. S3 a) .
In microdisks, doublets are present due to a scattering-induced coupling between degenerate clockwise and counterclockwise propagating modes. We fit this transmission spectrum to a model given by two equal lorentzian extinction curves: From the fit we can obtain the cold-cavity parameters describing the doublet: the splitting between the resonances (δ/(2π) = 11.3 GHz), the intrinsic losses (κ i /(2π) = 7.3 GHz), and the extrinsic losses (κ e /(2π) = 4.2 GHz, related to the coupling to the tapered fiber). It is also useful to define the total loss, κ = κ i + κ e . As the incident power increases, nonlinear effects change both the frequency and the linewidth of the optical doublet, giving rise to transmission spectra such as the red curve in Fig. S3 b). Results from the fit of such high power transmission spectra are unreliable, even for a simple model where nonlinear effects are polynomial functions of the internal energy. For that reason, we follow a method similar to the one used in [6] where the nonlinear effects are extracted from the transmission local minima, indicated by the circle and square markers in Fig. S3 b). For our analysis we collect the extinction (depth), E(P in ), and the frequency, δω E (P in ), of the local minima for all transmission spectra. In Fig. S3 c) we plot the variation of δω E (P in ) while in Fig. S3 d) the extinction.
In order to extract information regarding the nonlinear effects we make some assumptions for the case of an optical doublet: 1) the frequency splitting between the two resonances and the extrinsic loss, are independent of the cavity internal energy and 2) nonlinear effects are approximately equal for both resonances. Using such assumptions and the cold cavity parameters, the transmission becomes: where κ nl and δω are the, yet unknown, nonlinear loss and nonlinear dispersion.
Considering that the cold-cavity parameters are known from the fit and that δω simply shifts the frequency axis, the extinction of the local minima, E, in Eq. S73 are a function only of the nonlinear loss, E(κ nl ). As the cavity-taper coupling remains in the undercoupled regime throughout all measurements, we can restrict our analysis to a range in which this function is bijective. Inverting it we use the extinction data in Fig. S3 d) to obtain the nonlinear losses.
With the nonlinear losses we can determine the internal energy of the cavity at the extinction minima. In order to simplify our calculations we consider that the energy at the extinction minima is approximately equal the energy at the resonance points (∆ = 0 and ∆ = −δ): As our doublet is nearly resolved this approximation is good, the difference between such points is of approximately ∆ E − δω 1 GHz, much smaller than both the linewidth and the splitting of the resonances, rendering an error smaller than 1% at the final cavity energy. For the determination of the nonlinear dispersion we approximate the variation of the extinction minima, ∆ E (P in ), Fig. S3 d), by the variation of the resonance frequency itself, δω . Our approximation is good, rendering an error smaller than 250 MHz for the higher incident power. This happen because of two reasons, 1) we assumed δ is fixed and 2) the nonlinear loss variation is small when compared to κ and δ.
Since we assume that the linear losses and the nonlinear effects are equal for both modes of the doublet, data used in our analysis is obtained through an averaging process of the information collected for each of the modes in the doublet given an incident power.
It is worth mentioning that in our determination of the nonlinear effects we took the lowest power measurement as a frequency reference, leading to two issues: 1) despite being small, the nonlinear dispersion/losses are different from zero in our reference, and 2) the extinction data for low optical powers is considerably noisy, see Fig.S3c), such that we obtain negative values for the nonlinear loss for some of the experimental points, indicating an uncertainty at the cold-cavity losses. Both these issues can be solved with the introduction of offsets, δω nl,0 and κ nl,0 in our nonlinear model.

C. Fitting of the polynomial model
In the fitting process we consider a polynomial model for both the nonlinear losses and the nonlinear dispersion. First we fit the nonlinear losses, and then we use this information in the fitting of the nonlinear dispersion. This method is efficient since in our model we consider only the nonlinear dispersion due to the thermo-optic effect.
In our model we consider nonlinear losses scaling both linearly and quadratically with the internal energy, U : The free parameter, κ nl,0 , is an offset to the nonlinear losses used to correct the cold cavity intrinsic loss. The first order nonlinear loss coefficient, κ 1 , is related to two-photon absorption (TPA). It is calculated through FEM simulations where the material properties in [2] are applied to the model present in [6]. The second order nonlinear loss coefficient, κ 2 , is added as a phenomenological free parameter. First order nonlinear effects (related to TPA) are reliably negligible with respect to the quadratic effects, as seen in Fig. 4d) where the nonlinear loss displays a parabolic shape. This has also been verified by allowing the coefficient κ 1 to vary up to one order of magnitude during the fitting process, yielding remarkably similar results to those obtained using the FEM-estimated κ 1 .
With the knowledge of the nonlinear losses, we proceed to fit the nonlinear dispersion, given by: The free parameters in this step are δω nl,0 , the nonlinear dispersion offset, κ abs , the absorption loss rate, and η nl , the absorptive fraction of the nonlinear losses. The coefficient G θ s R θ s is the product of the thermo-optic dispersion by the thermal resistance, for static heating sources, it is calculated with FEM simulations. We choose to neglect the change on the thermal resistance due the nonlinear losses, because simulations show that for the static case it is negligible.
As seen in Fig. 4 c) and d), our model fits well the experimental data, with κ nl,0 = −1.0 ± 0.1 GHz, B = 1.0 ± 0.1 ns, δω nl,0 = −5.8 ± 0.4 GHz, κ abs = 5.050 ± 0.006 GHz and η nl = 0.83 ± 0.01. As a final test of the validity of our model we "simulate" our experiment by solving the transcendental equation for the internal energy of the cavity. We use the cold cavity parameters and our nonlinear model, obtaining the transmission as a function of the laser detuning. The agreement between the experimental transmission spectra and its "simulated" counterpart, Fig. 4 b), is nearly perfect.

A. Optical mode identification
In order to identify the optical mode that was excited in our experiment, we analysed the optical spectrum of our microdisk as shown in Fig. S4 a). Due to the symmetry of our device, modes can be classified through an azimuthal modal number m. Using a reference mode with frequency ω 0 and modal number m 0 , optical frequencies may be conveniently described as [10]: where µ = m − m 0 , D 1 /(2π) is the Free Spectral Range (FSR) and the D n /(2π), n > 1, give the FSR dependency on µ n−1 .  For the measured optical modes FSR ≈ 3 THz. We compare this value with simulations that incorporate both geometric and material dispersions for our device. Our measurements are found to be consistent with the 6-th radial order TE optical mode shown in Fig. S4 b). In Fig. S4 c1) (d1) the m dependency with the frequency is plotted for TE (TM) optical modes between 180 and 210 THz. For completeness, in Fig. S4 c2) (d2) we also show the residual dispersion D int = ω µ − ω 0 − µD 1 for the simulated TE (TM) modes. Results are summarized in Table II

B. Mechanical anisotropy
We incorporate the known mechanical anisotropy of GaAs in a 3D numerical model of the microdisk. This is done by correctly prescribing values for the stiffness tensor coefficients used in our simulations: (using Voigt notation) c 11 = 119 GPa, c 12 = 53.4 GPa, and c 44 = 59.6 GPa [4]. The symmetries of the system allow for correct calculations of the acoustic eigenmodes by considering only 1/4 of the full structure. All overlap integrals between optical, thermal, and mechanical modes were carefully evaluated by considering isotropic models for the thermal and optical responses. The former are calculated from 2D axisymmetric simulations for the device, which are then extruded into fully 3D fields. Only by taking in account the GaAs anisotropy we were able to reach an agreement between the simulation and the experimental data for the single photon optomechanical coupling, g 0 /2π = 29 kHz.

S8. OPTOMECHANICAL RESPONSE
In this section we discuss in detail how the experimental data on the optomechanical response was collected, analyzed and compared with the results expected from our model.

A. Data collection
For greater precision, all optomechanical measurements were performed by fine-tuning the laser frequency with a piezoelectric actuator. This choice limits the measurable optical frequency span to 30 GHz around a center frequency, which is chosen to render the peak modification on the mechanical linewidth within the measured range. Furthermore, the thermally-induced optical bistability restricts our experiments to the blue side of the optical resonance for the higher incident powers used in this work.
In this experiment we use a fast photodetector, with a 800 MHz bandwidth, to measure the mechanical fluctuations induced in the optical signal. The electrical signal of the photodetector is monitored with an Electrical Spectrum Analizer (ESA), which returns the Power Spectral Density (PSD) of its input signal. Fig.4 f ) in the main text is an example of a PSD where the fluctuations due to the fundamental breathing mode of our microdisk clearly rise above the detector's noise floor. As the optical fluctuations are very weak, each ESA measurement must be integrated for some seconds.
For each attenuation a regular, static, optical transmission spectrum is measured using the standard laser frequency sweep -as described in section S6 -, and a PSD map varying the laser frequency with the piezoelectric actuator, Fig.4e). A PSD map is composed of 100 PSD taken at equally spaced laser frequencies inside the 30 GHz piezo range. The laser spends approximately 3 s in each frequency, in such a way that building a PSD map for one attenuation takes approximately 5 min and 30 s. The attenuation is varied from 20 dB to 12 dB in steps of 1 dB. The entire experiment takes around 1 h.
As this experiment is considerably longer than the previous one, the tapered fiber position drifts moderately with respect to the microdisk, changing the extrinsic loss of our optical resonances. This is clearly seen in Fig. S5 a) where the extinctions of the optical spectra taken in the optomechanical characterization (slow experiment) and in the optical nonlinear characterization (fast experiment) are compared. For lower incident powers the extinctions of both are similar, but as the P in increases the extinction of the slow experiment becomes smaller than the one in the fast experiment. As the slow experiment is performed just after the fast one, this indicates that the taper position slowly drifts closer to the cavity as a function of time, increasing its optical coupling to the cavity.
We use the extinction information to correct for the extrinsic loss. Since the static optical transmission is recorded for each of the input powers in the fast experiment, we are able to recover its resonance frequency shift and compare it to the data obtained in the slow experiment. This gives us direct knowledge of the energy (at the extinction minima) in the resonator during the fast experiment. We use this information to determine the nonlinear loss, in such a way that the extinction becomes now a function only of the extrinsic loss. In a similar way to what was previously done for the determination of the nonlinear loss, we derive a numerical bijective relation between the extrinsic loss and the extinction, which we can invert to obtain the extrinsic loss, as shown in Fig. S5b). This estimate of the extrinsic loss during the optomechanical characterization is applied to our model for the photothermal effects.

B. Mechanical linewidth variation
Fitting the PSD to a model taking into account the mechanical susceptibility and a constant noise background, as in Fig. S5 c), we can determine the mechanical frequency and linewidth as a function of the laser frequency.
In order to study photothermal effects we analyzed the mechanical linewidth variation. In Fig. S5 f) two distinct behaviors are verified. For ∆ between 10 GHz and −12 GHz (clear region), the mechanical linewidth varies smoothly and presents very small uncertainty bars, while for ∆ < −12 GHz (shaded region), uncertainty increases considerably. As seen in Fig. S5 d), the clear region is found to match the blue detuned lateral of the optical resonance, where the high optomechanical transduction, shown in Fig. S5 e), leads to a large signal to noise ratio (SNR) in the measurement of the mechanical linewidth. Due to its larger precision, we focus our analysis on the clear region. The shaded region is related to the transmission minima, where the optomechanical transduction approaches zero, leading to a low SNR in the PSD.
Feeding our model for the photothermal effects with the optical nonlinearities, and with the extrinsic loss estimate, we can predict effects such as the mechanical linewidth variation. Fig.4g) and h) of the main text show our predictions as a function of both the laser detuning and the incident power at the cavity, matching very well with the experimental data. The value of the parameters fed to our model are listed in Tab.I as well as a brief explanation of how they were obtained. Here we must stress that, differently from other works, we do not fit the optomechanical data, instead we predict it independently using FEM simulations and the characterization of the optical nonlinearities in the resonator. As discussed in the main text, the use of multiple thermal modes is crucial to the level of accuracy obtained here.
The geometrical dimensions (R = 5.1 ± 0.1 µm and R ped = 0.6 ± 0.1 µm) of our device are roughly determined from optical microscopy data. The value for the radius used in all simulations is fine-tuned matching the mechanical frequency from FEM simulations with the experimental mechanical frequency.
As mentioned in the main text, the mechanical frequency modification is dominated by thermalization of the cavity. Qualitatively, the heat generated by the optical absorption softens the GaAs membrane, which in turn downshifts the mechanical frequencies. This competes with the positive mechanical frequency shift predicted due to optomechanical backaction (for blue detunings), as shown in Fig. 3 (a) of the main text. Plots for the measured transmission spectrum and mechanical frequency shift δΩ are found in Figs. S6 (a) and (b), where one readily observes that the trend followed by δΩ strikingly contrasts with the one followed by δΓ (Fig. 4 (g), main text). In fact, the transmission and δΩ display similar patterns. This is a direct consequence of the increase in the circulating power when the laser approaches the optical resonance, increasing the overall temperature in the resonator. Dataset:FEM and scripts files for generating each figure are available at Ref. [9].