Model and phase-diagram analysis of photothermal instabilities in an optomechanical resonator

A study of the phototermal instabilities in a Fabry–Perot cavity is reported, where one mirror consists of a silicon-nitride membrane coated by the molecular organic semiconductor tris(8-hydroxyquinoline) aluminum and silver layers. We propose a theoretical model to describe the back-action associated with the delayed response of the cavity field to the radiation pressure force and the photothermal force. For the case under investigation, the photothermal force response occurs on a timescale that is comparable to that of mirror oscillations and dominates over the radiation pressure force. A phase diagram analysis has been performed to map the stability of the static solution as a function of the control parameters. The model equations are integrated numerically and the time history is compared to experimental measurements of the transmitted field and displacement of the membrane. In both experimental and theoretical data we observe large amplitude oscillations when the cavity length is scanned at a low speed compared to the growth rate of the instability. The perturbation is found to evolve through three regimes: sinusoidal oscillations, double peaks and single peaks followed by a lethargic regime. When the cavity length is scanned in opposite directions, dynamical hysteresis is observed, whose extension has a power law dependence on the scanning rate.


Introduction
The study of the interaction between a light wave and the mechanical vibration of an object is known as optomechanics. The object dimensions can span from the nanoscale to macroscopic scales. The corresponding range of applications is very wide, including quantum computation with optomechanically coupled hybrid qubits, as well as macroscopic motion control such as in gravitational wave detection. Depending on the characteristics of the optomechanical system, it is possible to observe a variety of phenomena promising for both fundamental investigations [1] and applications [2]. The basic mechanism of optomechanics is the radiation pressure associated to the exchange of momentum between photons and the controlled material [3,4]. Experiments with dielectric membranes or metal/dielectric bilayers are based on photothermal processes, i.e. on membrane deformation by laser heating [5]. Moreover, the interplay of photothermal and radiation pressure forces is responsible for possible effects down to the quantum regime [6][7][8]. For a GaAs semiconductor membrane [9], the optomechanical coupling is provided by the non-radiative decay of photo-excited excitons.
The coupled light-object system is characterized by nonlinear interactions. As a consequence, in presence of a strong light driving the system, the optomechanical operation may reach regimes of bistability, instabilities and chaos. These regimes received a large attention [10][11][12][13], in some cases for the possibility to exploit the very large mechanical vibrations produced by the instabilities, but mostly because of the detrimental effects of instabilities on optomechanical sensors and the related need to avoid them. The first theoretical investigation of dynamical multistability was reported in [14] where the backaction associated to radiation pressure was analyzed in a cavity formed by a mirror and a cantilever. The study of self-sustained oscillations in optomechanical systems is addressed in [15], where the appearance of chaos is investigated in the classic and quantum regime. Following the early observations of the self-induced oscillations in [16], the presence of coexisting attractors and chaotic canard explosion was studied in [17]. In [18,19] the attractor diagram is investigated both experimentally and analitically; Hopf bifurcation and period doubling have been observed in [17,20]. The occurrence of transverse instabilities in an optomechanical set-up was considered in [21].
The present work investigates the instabilities and the phase portraits for an optomechanical cavity containing, as mirror, a three-layer membrane composed by the molecular organic semiconductor tris (8-hydroxyquinoline) aluminum (Alq 3 ), silver, and silicon nitride [22,23]. This material exhibits quite a lot of features similar to inorganic semiconductors and has been widely investigated lately for a large variety of optical applications [24,25], since they are easy to prepare in thin films and to integrate in liquid or polymeric environments. Its use in optomechanics would increase the application range of the material and is relevant for the possible production of optomechanical coupling by direct electronic processes [9] in organic materials. In the present work we study the optomechanical response of the cavity for a wavelength that is outside the Alq 3 absorption band and the main contribution to light absorption is due the silver layer.
As original feature of our investigation we interpret our experimental observations through the analytical and numerical solution of the full model for photo-induced forces. Our theoretical model describes the backactions associated to the delayed response of the cavity field to both photothermal and radiation pressure forces, where the membrane motion is described by an integro-differential equation containing an exponential delayed force whose delay is given by the time response of the photothermal process. For the case under investigation, the photothermal force response occurs on a timescale that is comparable to that of mirror oscillations. This differs from [17], which analyzed the slow-fast dynamics in a cavity where the vibrating mirror oscillation frequency was 10 4 higher than the photo-thermal response. The present work is motivated by the need to develop a more general model to describe our system, which is not based on the timescale separation between the photothermal and radiation pressure contributions. The approach can have application to a wider range of optomechanical systems, regardless of the timescale of the different contributions with respect to the optical field timescale. Our numerical analysis is based on the cavity parameters measured previously in [22,23], with a good agreement between simulations and measurements. The instability domains are determined by the system response to a slow modulation of the laser frequency at fixed input power. At variance with most previous experimental investigations of instabilities [16,[18][19][20], our investigation takes place in a 'dynamic' regime of the optomechanical system where the scanning rate is fast for the cavity parameters and slow for the membrane damping time constant. The experimental set-up enables measuring directly both the temporal evolution of the light field and of the membrane motion. The attractor phase portrait is reconstructed from the recorded temporal evolution using the standard technique based on the introduction of a higher dimension embedding space, where each point is obtained from a scalar set of measurements via a time delay. The instability investigation developed in our work is used to evaluate the amplitude of the photothermal optomechanical coupling, which accounts for the total photothermal force arising in our multilayer membrane, therefore providing an independent and more precise determination of the photothermal process strength. The latter is complementary to our previous determination based on the observation of the optomechanical cooling of the three-layer membrane [22,23].
Section 2 introduces the model of optomechanical system with the cavity, the membrane mirror motion, and the forces acting on it. Section 3 derives the differential system, based on five first order equations and written in dimensionless units, to be solved for the interpretation of the experimental results. Section 4 derives the fixed point solution and the bistability parameter range. Section 5 investigates theoretically the instabilities threshold conditions and their phase diagram for the measured experimental parameters. The following section 6 reports the instability pulse shapes observed experimentally on the cavity transmittance and on the membrane deflection. Phase portraits are reconstructed from those observables. The observation of dynamical hysteresis represented by the dependence of the instability range on the direction of the scanned cavity frequency is reported. Section 7 concludes our work.

Cavity
We consider a cavity with length L composed by the coated silicon-nitride membrane and the mirror as in figure 1. The cavity properties are determined by the reflectivity, absorption and transmission for mirror and membrane. The non-absorbing mirror has intensity reflection coefficient R 1 and transmission T R 1 The membrane displacement L D modifying the cavity length is produced by the total applied force F with an equation of motion given by where M g is the mechanical damping coefficient, M w the oscillation frequency in absence of damping, and m eff the membrane effective mass.
Following [1,14,27,28], the dynamical equation for the intracavity field amplitude α, normalized to the square root of the intracavity photon number at frequency L w , is where κ is the cavity intensity decay rate, is introduced, with c the speed of light.

Forces
The radiation pressure force F rp , determined by the momentum carried by each photon reaching the membrane, is given by the power divided by the light speed for each impinging laser beam [29]. For the configuration of figure 1 the radiation pressure force is F c P P P P ct For a large finesse the P in and P R contributions may be neglected, leading to the reported force expression on the basis of the photon number. The photothermal/bolometric force F ph is due to the absorbed intracavity power. Following the treatment of [16] we introduce a Λ parameter describing the ratio between the bolometric and radiative forces acting on the mirror. The finite time response of the forces plays an essential role. Reference [5] has introduced the ph t temporal response directly into the bolometric force through an exponentially damped delay function and the time derivative of the force, i.e. the so-called rigidity. Thus we write our time dependent photothermal force as where we define t=0 as the initial operation time, and this approach is valid for t ph t  . Notice that a similar expression containing a delay time is not required for the radiation pressure force, because the 1 k cavity delay time is included into the equation for the α temporal evolution. Cavity configurations with P in laser input power entering from right. P R is the reflected power, the transmitted power P T is monitored on the left.

Parameters
The parameters of our systems used in the numerical solutions are given in table 1. The large majority of those parameters are derived from our previous investigation of the membrane static displacement and laser cooling [22,23]. That does not apply to the Λ parameter characterizing the photothermal process that was derived in those references, 1000 L »with a factor two error bar. Table 1 reports the Λ value that provides the best agreement between the present experimental results and the instability simulations.

Optomechanical model
Using the dimensionless quantities We have here introduced the  dimensionless coupling constant [14] m c c P 4 . 9 Here the c g k term is replaced by a c rad radiative coefficient, introduced within [22]. This coefficient takes into account the electric field phase shifts on reflection and transmission from each membrane layer and provides a more precise determination of the radiation pressure force of equation (5). Notice that time-independent contribution to the force depends on 1 + L ( ), therefore positive in absence of photothermal processes, but changing sign with 1 L < -.
This system of partial differential equations represents our model for cavity optomechanics in the presence of radiative and photothermal forces. The Λ and ph t parameters determine the physics of the photothermal processes.
Equations (8) are written more conveniently introducing two additional variables, the v time derivative Therefore we obtain the following system of first order differential equations, including the temporal evolution of *  , conjugate of ò: This system written in dimensionless units is convenient for the following analyses based on the phase-diagram and on numerical solution. Introducing the variables , r i   , which are respectively the real and imaginary part of ò, the above equations become

Fixed points
The time-independent stationary solutions, i.e., fixed points, for equations (12), denoted as v x , , , , , are given by the following combined equations: where we have introduced the Z driving power The solution of the above equations is determined by cubic equations satisfied by both the x s and s 2  | | variables [30]. The cubic equation for the x s solution is written as E x 0 This third-order condition leads to a regime of optical bistability between two fixed points. Reference [30] has shown that equation (20) for E x s 3 ( )has three real roots only if the modulus of driving power exceeds a threshold value and is contained within an interval of the driving power depending on the detuning c d . In addition, the detuning should be larger than the threshold value. The c d bistability values are negative (red detunings) for positive values of Z and positive (blue detunings) when Z changes its sign. Owing to the large negative value of Λ, Z is negative for our case. The associated bistability diagram is shown in figure 2 The following characteristic equation can be written for equations (21): a a a a a 0, 22 , the coefficients are given by The stability boundary, derived applying the Routh-Horwitz criterion, is given by (see [31]) a a a a a a a a a a .

Growth/decay rate
In order to obtain additional information on the time-dependent solution for different parts of the map, we have integrated numerically the set of differential equations (8) at given values of c d and P in and different initial conditions for x 0 t = ( )and the corresponding 0 0 Within the instability regions the electric field amplitude increases at a fast rate until a saturation is reached as shown in figure 4(a). On the contrary the membrane oscillation amplitude increases with lower rate, of the order of 10 s −1 , as shown in figure 4(b), the growth rate being determined by the m g damping rate. Notice that the oscillation amplitude is larger than the static x s value, derived from the fixed point solution, by a factor 10-10 2 , corresponding to peak values x = 1-10 within the instability regime. As an additional test of the stable solutions, we have imposed, for c d outside the instability region, a large initial value x 0 5 t = = ( ) for the initial membrane displacement and verified the evolution towards a steady state solution with a decay rate having the same order of magnitude as the growth rate of the unstable solutions.    figure 5. This asymmetric scan of the instability regime was reported also by Metzger et al [16] and by Marino and Marin [17]. A fair agreement between experiment and model exists for the positive slope, but for the negative slope the experimental results show a reduced extension of the unstable region. As we approach the detuning value for which the periodic attractor would disappear, the experimental noise and the thermal fluctuations are more likely to push the system out of the periodic attractor into the the fixed point solution basin of attraction, causing the observed collapse to the fixed point solution.
Notice that the instability regimes predicted by the simulations occur for detuning values different from those predicted by applying the Routh-Horwitz criterion of section 5.2. This difference appears from a comparison between the horizontal line in the phase diagram of figure 3, extending from 8.45 inst,min , and the actual instability limits, defined by the behavior of t 2  | ( )| . We define the instability detuning intervals respectively as ( inst,min The extremal detuning values of the instability regions are plotted in figure 6(a) as a function of the absolute value of the cavity detuning sweep rate. It appears that inst,max d + , the value for which the numerical solution becomes unstable during a positive slope sweep, depends strongly on b | |, while inst,min d -, corresponding to the return to the stable solution at the end of a negative slope sweep, loosely depends on b | |. The apparent discrepancy between the instability region predicted by the linear stability analysis and what observed both in the time dependent numerical integration, when we apply a negative sweep, and in the experiments, has a simple explanation. The linear stability analysis, as carried out in the manuscript, refers to the stability of a (static) fixed point, and tells us when this static solution is no longer stable: it cannot predict stable dynamical attractors (the large periodic oscillations), which we observe as we sweep down the frequency (b 0 < ). These oscillations are due to the nonlinearity present in the system. In other words, the observed (much wider) hysteresis region is a joint result of nonlinearity and sweeping, and it cannot be described by a static analysis. Besides the analytical interpretation, from a physical point of view these nonlinearity-born periodic oscillations pump energy in the system, supporting a more energetic optomechanical state. The above response of our optomechanical system has a strong analogy with the dynamical hysteresis investigated within [32][33][34] where it was demonstrated that the bistability range depends on the sweep rate of the explored parameter. Notice that such hysteresis was at first presented as a modification of the bistability range for a bistable system. However [34] pointed out that dynamical hysteresis appears also in systems not bistable, but with an operation regime close to a bistability range. This situation applies to our system. In that reference it was pointed out that the deviation of the hysteresis range from its static value depends on the sweep rate following a power law, with the exact power depending on the specific system. For the present analysis, we define a dynamical hysteresis range H inst,min inst,min inst,max inst,max  additional attractors forcing the system into an unstable regime, as pointed out in [17], will be tested in future works.

Pulses
In the experiment and simulation at the given input laser power (P 0.5 in = mW) and scanning the c d cavity detuning, as shown in figure 5, the instability appears as a sinusoidal oscillation at a frequency close to m w , as observed in [17,18] and as predicted by equation (25). This oscillation behavior takes place at the Hopf bifurcation, both when entering the instability region from low and high detunings. This sinusoidal oscillation regime, denoted as regime 1, is shown in figures 7(a) and (b), with the experimental results for the cavity transmitted power P t T 2  µ | ( )| on panel (a) and the simulation for t 2  | ( )| on panel (b). In contrast with [35], the Hopf bifurcation is not followed by a period doubling cascade; the period doubling is not visible either in the cavity transmission data nor in the simulation results for 2  | | . We observe instead a regime, denoted here as regime 2, which is characterized by a repetition of double peaks whose periodicity is the same as for regime 1. This is shown in figure 7(c), for the experimental data and (d) for the simulation. Following the laser chaos denomination [36,37] this regime can be denoted as T 11 . The final regime 3 is shown in figures 7(f) and (g) for experimental and simulation, respectively. Here the single peaks are followed by a lethargic phase. If the detuning is reduced, progressing from positive to negative values, the optomechanics instabilities observed in the cavity transmission start with regime 1, evolve through regime 2 and 3, and finally decay. If instead the detuning is increased, progressing from negative to positive values, the instability starts from regime 1, evolves through regime 2 and 3, goes back to regime 1, and finally decay. The detuning values for which the different regimes are reached are indicated on the top of each panel in figure 7 for the case of decreasing detuning (i.e. negative sweep) corresponding to figures 5(e) and (f) for the experiment and the simulation, respectively. The experiment and the numerical simulations have very similar dynamical evolution when similar regimes are compared (close to the edges of the instable region, in the central part of the instable region).

Phase portraits
Applying a time delay approach to the temporal evolution of the P T cavity transmission, either experimental or from the simulation, we have reconstructed the attractor as presented in figure 8. The time delay τ utilized in the reconstruction of the phase space is estimated from the position of the first minimum in the average mutual information of signals P T (t) and P t T t + ( ), as a function of τ [38]. Regime 1 in figure 8(a) is characterized by a circular evolution within a 3D space, regime 2 in panel (b) by a torus-like structure, more developed for the experimental data, and less developed for the similar numerical data not reported here. Finally the lethargic regime 3 corresponding to the peaks in the cavity transmission reported in figures 7(f) and (g), is characterized by a large evolution within a plane, followed by an evolution within a perpendicular plane as shown in panel (c).

Conclusions
The present work investigates the instabilities occurring within an optomechanical system based on a three-layer membrane where the photothermal processes are more effective that the radiation force processes and therefore play the key role. These instabilities are analyzed on the basis of a model for photothermal forces previously introduced by Metzger and Karrai [5]. We have rewritten the equation system describing the optomechanics response in a convenient form for the instability analysis, and we have applied analytical and numerical solutions to characterize the unstable region. This model is compared to a large series of experimental results on the basis of the system parameters obtained from previous measurements of the membrane laser cooling, except for the photothermal process Λ parameter. The main outcome of our detailed study of the instabilities, in both experiment and simulation, is the determination of the Λ parameter that was measured with a low accuracy in our previous investigation of the membrane laser cooling [22,23]. We have used in our simulations 1250 L = -, slightly larger than the value 1000 L =derived in the previous laser cooling investigation, because it represents the minimum value where the simulation predicts the presence of a dynamical hysteresis for the instabilities at the large ±5 s −1 scanning rate, in units of the cavity linewidth. A larger Λ value leads to hysteresis cycles much larger than the experimental observations. For that value, the observed pulse sequence and their phase portraits are well represented by simulations. Instead the predicted ranges of dynamical hysteresis do not correspond to the experimental observations. We have verified that a 10 percent variation of all other parameters inserted into the model cannot Figure 7. Temporal dependence of the optomechanical instabilities observed on different theoretical and experimental variables versus time. From the top to the bottom rows, regimes from 1 to 3, as defined within the text. Experimental data for the cavity transmission P T in red color of the left column and simulation in blue color for 2  | | on the right column. All data are for P 0.5 mW in = and for the cavity detunings value reported on top of each plot.
explain the remaining discrepancies. The instability domains depend on both Λ and ph t in a complex manner. We have verified that a 10% variation of ph t (corresponding to the experimental uncertainty) around the value 0.1 ms modifies these instability regions in a comparable way as a variation of 10% of Λ (around the −1250 value) does. On the basis of our simulations we estimate the Λ range 1125, 1375 - -[ ] . In order to improve the experiment/simulation match, it will be important to verify either the presence of a small optical coupling of the input laser to additional cavity modes, or the oscillation of the membrane containing additional acoustic modes. A further important question is the experimental presence of additional photothermal processes specifically related to the investigated multilayer system, and not described by the applied model. In order to reach a better characterization of the membrane photothermal processes, and verification of other, possibly neglected processes, it will be interesting to use a bottom up approach and derive Λ values from the analysis of the photothermal processes contribution to the mirror displacement, as applied in [39,40].
The very large extent of the instability region associated to the case of negative sweep, in both the numerical simulations and the experimental observations, similar to the report by Marino and Marin [17], appears to be the result of nonlinearity and sweep. From preliminary investigations, the inclusion of thermal noise strongly affects the boundaries of the unstable region and the amplitude of the oscillations in transient regimes, leading to an overall shrinking of the region where bistability is observed.