Accurate modeling and characterization of photothermal forces in optomechanics

Microscale photonic devices have 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–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,⁶ radiation pressure (RP),⁷ and piezoelectricity.⁸ 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⁴) 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, thermally induced 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¹⁷ resonators. This extra degree of freedom opens up new possibilities, such as thermally mediated optomechanical ground state cooling in the bad-cavity regime.¹⁴ 

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¹⁸ 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.

The mechanical system is described by the equation of motion $ρU⃗̈=∇⋅σ$⁠, where $U⃗$ denotes the displacement field and σ denotes the stress tensor. In thermo-elasticity, the self-consistency 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,¹⁹ it is necessary to split the strain of the system, $S=12(∇U⃗+∇TU⃗)=∇̂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^θ,²⁰ 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 Sec. S1 of the supplementary material. In some previous formulations of PTh forces in optomechanical systems, this subtlety has been neglected²¹ and can lead to inaccurate predictions of the PTh response of optomechanical resonators.

The volume load $f⃗θ$ was used in past work on optomechanical PTh forces to provide an estimate of their magnitude.²¹ 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,³⁰  $δ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 Sec. S2 of the supplementary material. In this framework, the PTh force for the nth-mechanical can be written as a sum of the contribution from multiple thermal modes as $Fnθ=∑kΛk,nθθk$⁠, where $Λk,nθ=∫Snx:(c:α)δTkdV$⁠. Similarly, the surface and volume contributions can be decomposed in terms of $Λk,nθ-Vol.$ and $Λk,nθ-Sur.$⁠.

In order to numerically validate the present model, we first consider the case of static thermal deformations on a GaAs on Al_0.7Ga_0.3As (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 following calculations.

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(a); the thermal displacement field components along $r̂$ and $ẑ$ directions are shown in Fig. 1(d). 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.

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 supplementary 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 $Λk,nθ−Sur.$ and $Λk,nθ−Vol.$ 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.³² The total loss rate (κ = κ_e + κ_abs + κ_non-abs) is κ/(2π) ≈ 1.93 GHz, with an extrinsic coupling rate (i.e., coupling to a waveguide) κ_e/(2π) ≈ 0.48 GHz. These numbers yield a loaded quality factor Q_opt ≈ 10⁵; all other parameters are obtained through FEM simulations, where the 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 $Fnθ$⁠. 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 six 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 a phenomenon can be observed is if surface and volume contributions counterbalance each other. We illustrate this in Fig. 2(c), 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. Finally, since the 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.

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 3(a2), the PTh and RP backaction curves are displayed. For the RP calculations, both photoelastic³³ and moving boundary contributions³⁴ 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 out-of-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 to 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 a 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) and 3(a2).

The effectiveness of the thermodynamic description is tested by monitoring the modification on the mechanical linewidth of a cavity optomechanical system consisting of an R ≈ (5.1 ± 0.1) µm, R_ped ≈ (0.6 ± 0.1) µm GaAs/Al_0.7GA_0.3As (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 supplementary 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 The 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 non-absorptive 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³) in internal energy are enough to describe our results. A comprehensive guide for this analysis is found in S6 of the supplementary 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. Figure 4(d) displays the total optical dissipation rate of the system split in absorptive (κ_abs) and non-absorptive (κ_non-abs) parts.

We measure the backaction 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 sixth order TE optical mode (S7A of the supplementary 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 supplementary material) and the experimental nonlinear loss and dispersion described above. Importantly, our FEM results—where the parameters $Rkθ$ and $Λkθ$ were evaluated—require only thermal and thermo-elastic material properties as input, all of which are well established in the literature.³⁷  Figure 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)–4(h4) for four different input powers. Details on the collection and analysis of data regarding the mechanical response are found in S8 of the supplementary material, along with the treatment of the stiffening of the mechanical oscillator (δΩ), which is dominated by a static temperature softening of GaAs,³⁸ 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)–4(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.

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. In addition, the modal treatment for the thermal response is a significant step toward 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 that photothermal forces can be appreciable in other platforms and geometries and that the content of this work can provide insight into those cases.

See S1 of the supplementary material for details on the derivation of the photothermal coupling using thermodynamics and the mechanical equation of motion. S2 provides the thermal mode analysis of the GaAs microdisk. The frequency domain description of the photothermal coupling is provided on S3 for single-mode, multimode, and nonlinear multimode photothermal coupling. In S4, we discuss the impact of the localized absorption for our proposed model and our current experimental implementation. A detailed description of the experimental setup is provided in S5 followed by all the fitting functions used for the nonlinear optical response on S6. Section S7 discusses the modal identification for both the optical and mechanical modes by comparing experimental and finite element simulations. Finally, in S8, a detailed discussion of the optomechanical response is presented.

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. This work was supported by the São Paulo Research Foundation (FAPESP) through Grant Nos. 2019/09738-9, 2020/06348-2, 2017/14920-5, 2016/18308-0, 2017/19770-1, 2018/15580-6, 2018/15577-5, and 2018/25339-4; Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES) (Finance Code No. 001); Conselho Nacional de Desenvolvimento Científico e Tecnológico through Grant Nos. 425338/2018-5, 310224/2018-7, and 465469/2014-0; Financiadora de Estudos e Projetos (Finep); and the Natural Sciences and Engineering Research Council (NSERC) of Canada.