Observation of a phononic Mollow triplet in a multimode hybrid spin-nanomechanical system

Reminiscent of the bound character of a qubit's dynamics confined on the Bloch sphere, the observation of a Mollow triplet in the resonantly driven qubit fluorescence spectrum represents one of the founding signatures of quantum electrodynamics. Here we report on its observation in a hybrid spin-nanomechanical system, where a nitrogen-vacancy spin qubit is magnetically coupled to the vibrations of a silicon carbide nanowire. A resonant microwave field turns the originally parametric hybrid interaction into a resonant process, where acoustic phonons are now able to induce transitions between the dressed qubit states, leading to synchronized spin-oscillator dynamics. We further explore the vectorial character of the hybrid coupling to the bidimensional deformations of the nanowire. The demonstrated microwave assisted synchronization of the spin-oscillator dynamics opens novel perspectives for the exploration of spin-dependent forces, the key ingredient for quantum state transfer.

Supplementary Figure 2: Schematics of the energetic structure of the NV levels in absence (left) and presence (right) of an external magnetic field. In presence of an external magnetic field, the new eigenstates |Ψ i become mixtures of the initial eigenstates Ψ 0 i . As a consequence the transition rates k ij are also modified. This analysis permits to estimate the variations of the mean fluorescence rates and the evolution of the ESR contrast (see text).

SUPPLEMENTARY NOTE 1: EXPERIMENTAL DETAILS
The experimental setup combines in a volume of a few microns, the NV functionalized nanowire, the magnetic field gradient source, the microwave antenna and the electrostatic actuation. An extremely high stability was required to operate the system: in a magnetic field gradient of 10 5 T/m, moving the nanowire by 1 µm shifts the static magnetic field experienced by the NV spin by 100 mT, and represents a 3 GHz frequency shift. Requiring a drift of less than 3 MHz imposes a stability of 1 nm over the measurement duration, which can last for hours. The strategy adopted was to keep the NV-functionalized nanowire fixed in space and to organize the experiment around it. We have developed an overhanging MW waveguide, supporting a magnetic microstructure providing the strong magnetic field gradient. It is also used to apply an electrostatic force to drive the spin-functionalized nanowire into motion. The optical microscope objective used for injection and for collecting the reflected pump light as well as the NV fluorescence was mounted on piezo actuators to track the NV defect in space. The whole experiment was mounted in a homemade vacuum chamber. Note that the results reported here, except for determination of the vibration orientations, were obtained at ambient pressure. However, the chamber nonetheless played a crucial role in ensuring a stable operation of the system.

Optical setup
The overall layout of the optical table is sketched in Supplementary Figure 1. The laser source employed for fluorescence measurements and optomechanical readout is a 532 nm diode pumped solid-state frequency doubled Nd:YaG laser (Laser Quantum, Torus, 200 mW maximum output power). A confocal fluorescence microscope arrangement is implemented, combined with the possibility to measure the reflected and transmitted intensities on quadrant photodiodes to read out the vibrations of the NV-functionalized nanowire. Prior to injection in the cahmber, the laser beam diameter and polarization are adjusted and its intensity modulated with an AOM in double pass followed by a fiber cleaning of the optical profile. The AOM is controlled by an arbitrary TTL pulse generator (SpinCore, PulseBlaster 500 card) in order to allow ns timing of optical pulses. The laser enters the vacuum chamber through an optical window (covered with a broadband anti-reflection coating) and is reflected by a dichroic mirror located inside the cavity. Preventing the pump beam and fluorescence channel to traverse the same optical window ensures that no additional parasitic fluorescence is detected. The laser is focused onto the functionalized nanowire and the emitted NV fluorescence is collected through a high numerical aperture microscope objective (model EC plan -NEOFLUAR 100 ×0.75 from Zeiss), with a 4 mm working distance. The fluorescence exiting the chamber by an optical port is then filtered to select the NV emission band, subsequently split into two paths by a polarizing beam splitter and detected by two APDs. An active electronic splitter based on logic gates duplicates the photon counts with variable delays. The different outputs are interfaced with a Nanonis scan controller, which controls the laser position, to a PCIe 6323 board (National Instruments) for ESR measurements and to a FastComTec MCS6A card for time resolved fluorescence detection, needed for the measurement of Rabi oscillations. The two APDs also permit to realize autocorrelations measurements of the emitted fluorescence. A white light source and a CCD camera (Watec 101N+) can be used for monitoring and permit a coarse alignment of the nanowire, waveguide and magnetic gradient source. The reflected intensity is collected through a 90/10 non-polarizing beam splitter which is placed on the injection path, and focused on a low noise quadrant photodiode amplifier. Similarly a second microscope objective can be installed facing the injection objective to collect the transmitted light. The DC output of both detectors can be used for fine positioning of the experiment, while the differential HF output serves for readout of the nanowire vibrations [1].

Microwave -Radiofrequency apparatus
Spin manipulation is enabled by a MW antenna which also holds the magnetic gradient source. It is brought in the vicinity of the functionalized nanowire by means of a closed loop Physik Instrument Nanocube, moving by 100×100×100 µm. The MW source is a SMB100A Signal Generator from Rohde & Schwarz. Its intensity is externally modulated by a fast microwave switch from Minicircuit (ZASWA-2-50DR++, coaxial High Isolation Switch) which is TTL controlled by the PulseBlaster card that governs the time sequence of dynamical spin measurements. When needed, a 25 W MW amplifier from Amplifier Research was used. Two vacuum feed-throughs with SMA connections are employed to bring MW and RF signal in and out of the vacuum chamber. An external 50Ω 2W load was placed on the output port, to avoid power dissipation inside the chamber. To implement the electrostatic actuation, two bias tees are placed upstream and downstream of the chamber MW/RF ports. When needed, a low noise RF amplifier was used (Falco Systems), the source being either a network analyzer (Agilent 5061B) for response measurements (such as in Fig. 1f) or an arbitrary signal generator. The RF output of the second bias tee was left open in order to impose an electrostatic voltage on the nanowire and reduce RF currents at the nanowire position.

Micro-positioning and position control
As already mentioned, the functionalized nanowire serves as a spatial reference in the experiment. It is mounted on a coarse XYZ translation stage (Microcontrole, 562) equipped with piezo motors (Newport, Picomotors) that serves for coarse positioning of the nanowire on the optical axis. The injection objective is mounted on a XYZ piezostage (Physik Instrument, 733.DD) with a 30 × 30 × 10 µm scan range. A Nanonis scan engine serves for imaging the system and tracking the NV fluorescence in time to compensate the thermal drifts with respect to the optical axis. The piezo stage holding MW waveguide and magnetic structure is mounted on a second XYZ translation stage, that also supports the translational stage holding the nanowire. As such, after a preliminary alignment phase, where the motorized stages are involved, the nanowire ends within the scan range of the optical microscope, while the magnetic gradient source can be micro-positioned in the vicinity of the functionalized nanowire. Once this configuration is reached, all of the positioning parameters can be derived from physical measurements and optical monitoring is not needed anymore. The piezo stage holding the magnetic field gradient is controlled by a National Instrument (PCIe6323) board interfaced by a homemade scanning software.

SiC nanowires functionalized with a NV spin qubit
Long SiC nanowires with large aspect ratios and oscillation frequencies in the 100 kHz-10 MHz range exhibit large zero point fluctuations. Here we seek to demonstrate the Mollow triplet phenomena of the driven NV spin dynamics in the hybrid setup and to expose aspects of vectorial coupling in the system. Hence, we are interested in working in the spin dynamical resolved sideband regime, where the mechanical oscillation frequency surpasses the spin decay rates. Consequently, we have chosen to work with rather short nanowires. In order to attach a single NV defect to the tip of the SiC nanowire [2], we use the arrangement of the position controllers for the nanowire and waveguide of the hybrid setup. The nanowire is mounted in a copper sample holder which is fixed to the translational stage. A drop of the nanodiamond solution (Microdiamant, 50 nm average diameter) is deposited on a Silicon wafer which is put on the same piezo stage that otherwise holds the waveguide. White light illumination of the detection volume of the confocal microscope and the CCD camera facilitate navigation of the drop of NV solution, which is steered upward until its surface is penetrated by the nanowire. We then direct green laser light onto the nanowire tip. This creates heat convection favoring a flow of nanodiamonds near the resonator. In order to limit the number of nanodiamonds attaching themselves to the wire by a surface effect, the process is interrupted after a few seconds by moving the wire out of the solution. We then use the confocal microscope and the ability to perform fluorescence scans to inspect the nanowire . If NV defects have attached to the nanowire, they become visible as bright spots in these measurements. By measuring an autocorrelation function of the light emitted from the fluorescence spots, their single photon source character can be verified. Ultimately, an ESR measurement and detection of the spin resonance near 2.87 GHz verifies that the emitters are indeed negatively charged NV defects. If no NV defect is found, the process is repeated. Despite this approach being completely heuristic, the NV yield is usually quite high and after a few iteration a single NV defect can be attached to the nanowire extremity.

Nanowire actuation
Two different actuations were used in this work: piezoelectric and electrostatic. The piezoelectric actuation was performed by driving the nanowire through a thin piezo element (200 µm thick) glued on the nanowire support. It can be used even in absence of the RF/MW electrode, however the orientation of the corresponding force vector can neither be anticipated nor tuned. A drawback is the presence of parasitic resonances in the system response, which is not flat over the frequency range required to drive both mechanical eigenmodes (in particular in air where the damping rate is rather large). The electrostatic actuation is exerted by applying a potential between the nanowire metallic support and the waveguide, through the bias-tee inputs. The orientation of the electrostatic force vector (e F ) can be adjusted by moving the waveguide position with respect to the nanowire. Its frequency response is extremely flat over the range of interest. This is why this strategy has been employed in the measurements illustrated in Fig. 4.
However it is responsible for a parasitic modulated magnetic stray field which partially contributes to the observations of Fig. 4. It has been taken into account in the fit of the experimental data as an additional term independent of the drive frequency Ω d , inserted in the argument of the absolute value of equation (3), with a magnitude of 1.86 MHz and a phase of -81 degrees. This parasitic contribution is not present in case of piezo-electric actuation which does not create parasitic RF currents at the NV location.

SUPPLEMENTARY NOTE 2: FITTING OF THE MULTIMODAL AND 2D MOLLOW TRIPLET
The data of Figure 4a were fitted using the half the norm of the following expression: where the last complex term contains the magnetic stray field contribution visible when using the electrostatic actuation. The RF field that is employed to mechanically drive the nanowire (combined with a DC offset), is flat over the excitation frequency range, so that we can assume and verify that its phase and amplitude are frequency independent (the suspended RF-MW waveguide employed has a flat response in this frequency range). Its magnitude is δω bckgnd 0 /2π = 1.86 MHz, with a dephasing φ bckgnd of −81 o . The argument of the above expression can then be used to measure the dephasing between the parametric modulation undergone by the spin and the applied electrostatic voltage modulation (the electrostatic force follows instantaneously the actuation voltage). The other fitting parameters, in particular the orientations of the force and hybrid coupling vectors are reported in Figure 4b (220 • and 50 • respectively). They are in agreement with the spatial map shown in Figure 2f at the measurement point, the force vector pointing towards the microbead as expected for electrostatic actuation. This orientation is also in agreement with the individual phase contributions from each of the modes, shown in Figure  4c, which are identical and equal to −180 • at low frequencies. In this situation, when the actuation force increases, the nanowire gets attracted towards the microbead. As a consequence both mode see the same reduction of the spin qubit energy. The situation would have been different if the e F , e λ vectors would have lied in different quadrants with respect to the eigenmodes orientations ( e 1 , e 2 ).

Calibration of the parametric hybrid coupling strength
This section aims at describing the static calibration of the dynamical parametric coupling strength in further detail. As outlined in the article, the ability to measure the ESR measurements at different positions in space allows to map the spatial dependance of the NV spin resonance ω 0 (r) in the field gradient of the magnetic particle. We introduce the theoretical formalization supporting our method, and present the numerical magnetic simulations which permit modeling the spatial dependence of the vectorial parametric coupling strength used to generate Fig. 2d.

Static spatial dependence of the spin qubit energy
The static Zeeman HamiltonianĤ Z = −gµ Bσz B(r) can be linearized around the chosen working point r 0 to obtain the hybrid spin-mechanical Hamiltonian: with: Whereω 0 is the NV spin zero field splitting (2.87 GHz). The spatial variation of the magnetic field render the spin qubit eigenstates position dependent. They are labeled |Ψ i (r) , with i = −1, 0, 1 (so that in absence of magnetic field they converge towards the usual spin 1 eigenstates of the NV). Following the alignment procedure described in the main text (preserving the maximum florescence rate) permits working at a position (r 0 ) where the field generated by the magnetic particle at the rest position is aligned with the NV quantization axis B(r 0 ) = B z .e z (the following calculations are however independent of exact orientation of the magnetic field at the rest position). In that case, the eigenstates of the spin qubit at the rest position are identical to the inital states. The position dependent spin qubit energy corresponds to which can be numerically estimated for a given magnetic source (see the following section). In practice the ordering of the energetic level structure does not vary in the region of interest, so that we can continuously follow the ESR resonance having the lowest frequency.

Vectorial parametric coupling strength
When the spin qubit moves in space, the qubit eigenstates |Ψ i (r 0 + δr) will vary and result in a mixture of the eigenstates obtained at the rest position |Ψ i (r 0 ) . We verify here that the intuitive static calibration of the parametric coupling strength (based on the gradient of the measured ω 0 (r) map) works as expected when one takes this intermixing into account. The parametric modulation strength δω 0 (r 0 ) is defined as: The first term can be split as Z , but at first order in δr we have: We have: Reproducing the same reasoning with the |0 state, this validates that we are only left with the interaction term: The first term is proportional to Ψ 1 (r 0 + δr)| σ i δr j ∂ j B i |Ψ 1 (r 0 + δr) , which can also be written as δr · . The bracketed expression is nothing but the gradient of the total static Hamil-tonianĤ 0 +Ĥ Z . Since the eigenfunctions |Ψ 1 (r 0 + δr) are eigenfunctions of H 0 + H Z , which does not depends upon time, the Hellmann-Feynmann theorem can be used, which states for an Hamiltonian implicitly dependant on µ (the position in our case): Therefore one writes the parametric coupling strength as: and we formally recognize here the spatial gradient of ω 0 (r) (see equation (6)). As such, the parametric interaction hamiltonian can be expressed as: using the vectorial coupling constant λ: as employed in the article.

Magnetic simulations
We have performed numerical simulations of the spatial dependence of the spin energy E(r 0 ) to model the energy maps obtained in figure 2d, 2e of the article. The NV fluorescence was experimentally recorded as a function of the spin position in the magnetic field gradient in presence of a fixed microwave tone. To be able to reproduce both the position of resonance slices and the spatial dependance of the fluorescence, a detailed model of the NV center must be considered. Following [3], the fluorescence properties are fully understood using a seven level model of the NV center, see Supplementary Figure 2: a spin triplet for both ground and excited states and a dark metastable singlet state labeled |Ψ 0 1→7 .
The spin ground states |m S = −1, 0, 1 = |Ψ 0 1,2,3 and excited states |m S = −1, 0, 1 = |Ψ 0 4,5,6 are respectively eigenstates of the Hamiltonian H gs = ω gs .e z is aligned along the quantization axis e z , the energetic structure of the spin is preserved, as well as the transition rates between eigenstates. They are defined by k 0 ij (r 0 ). The NdFeB particle creating the magnetic field gradients is an homogeneous sphere. Its magnetic stray field is therefore exactly modeled as a dipole of magnetic moment M = (M V ).e z , where M is the zero field magnetization and V the volume of the particle: In the general case when a transverse magnetic field is present the previous Hamiltonians are no longer diagonal which modify the eigenstates as follow: The coefficients α ij (r 0 ) and the new eigenfrequencies ω i→j (r 0 ) are numerically computed with the assumption that |Ψ 7 = |Ψ 0 7 . As a consequence the transition rates between the new eigenstates are given by: The field dependent NV fluorescence rate can therefore be calculated at each position: where the averaged population of each statesn i is retrieved by solving the rate equation of the system in the steady state: The microwave field with a fixed frequency ω M W is only able to drive the |Ψ 1 → |Ψ 2 transition at the position where the resonant condition ω 0 (r 0 ) ≡ (E 1 (r 0 ) − E 1 (r 0 )/ = ω M W is met. At this position the fluorescence rate is reduced by the ESR contrast: At each spatial position of the NV center in the magnetic field B(r 0 ) the optically detected magnetic resonance signal is modeled with a Lorentzian curve L(ω, r 0 ) of amplitude F (r 0 ).C(r 0 ) and linewidth Γ S centered at ω = ω 0 (r 0 ). The resulting fluorescence rate at position r 0 is given by F (r 0 ) − L(ω M W , r 0 ), which allows to reconstruct the fluorescence map of figure 2e.
The typical zero field transition rates used in this simulation, extracted from the literature [3], are given in the following table: This simulation is able to reproduce and understand all the experimental features which appear in Figure 2d and validates the static estimation of the parametric coupling strength calibration procedure presented in the article. The distance to the bead being fixed, the only free parameter being the magnetic moment of the NdFeB particle, which is evaluated to be M 1.5 10 −9 A.m 2 . This value is compatible with the magnetization specified for this material (M = 5.5 10 6 A.m −1 ) which gives for a radius of R = 9 µm: M = M 4 3 πR 3 = 1.7 10 −8 A.m 2 .

Bloch equations
The goal of this section is to establish a set of Bloch equations that we will use to describe the evolution of the coupled spin-oscillator system.
We will consider the magnetic field coupling of a spin 1 system to a single mechanical mode in a most generic approach. Assuming that the mechanical oscillator is classical, the 2D vibrations of its oscillating extremity are denoted by δr. We consider the quantum master equation for the density matrix ρ: where the Hamiltonian H is written as which includes the coupling term The evolution of σ can be written: where we have separated the contribution of the coherences originating from the |−1 state: The last step allows us to restrict our analysis to the {|0, |1 } subspace for the evolution of the spin qubit instead of a spin 1 three level system. Such a simplification permits a straightforward interpretation of the evolution of the spin populations and coherences on the Bloch sphere. This restriction to the {|0, |1 } subspace is experimentally justified in the conditions described in the main text. The Hamiltonian can then be reduced to: while the damping contribution iṡ where Γ p conveys the role of the light induced polarization of the spin qubit through optical pumping [4]. The time evolution of the density matrix elements ρ ij = i| ρ |j , given byρ = 1 i [H tot , ρ] +ρ relax can be expressed as: We have used the time dependence of the resonant MW field Ω R (t) = Ω R cos ωt (it is assumed to be aligned along the x direction, but the physics does not change in case of a y orientation). We have introduced the parametric modulation strength δω 0 (t) as well as the resonant contribution δω R (t) due to mechanical oscillation in the magnetic field gradient which both vary at the timescale of the oscillator. By going into the MW rotating frame, we define: We introduce the coordinates of the Bloch vector: and obtain the Bloch equations: (34)

Numerical simulations
The previous set of equations are used to numerically simulate the dynamics of the spin qubit in presence of mechanical motion which generates parametric energy modulation. The time evolution of the population σ z (t) can then be simulated, and its Fourier transform is then calculated. The presence of light can be turned on (for ESR measurements) or off (for the Rabi dynamics) by turning Γ p on or off. To realize the numerical simulations displayed in figure 4d of the manuscript we proceed as follow. The parametric modulation strength δω 0 [Ω d ] is deduced for each drive frequency Ω d as explained in the manuscript by analyzing the driven trajectories in the parametric coupling landscape. At t=0 the spin qubit is supposed to be optically initialized in its ground state, (u = 0, v = 0, w = −1/2), the mechanically induced dynamical change of the MW-spin qubit interaction strength is neglected here (δω R =0), the light is turned off (Γ p = 0) and we emloy a time modulated parametric modulation strength: δω 0 (t) = δω 0 sin(Ω d t + φ), φ being an arbitrary dephasing parameterizing the state of the parametric modulation at the beginning of the Rabi measurement protocol. Then the 3 coupled Bloch equations are solved with a Runge-Kutta algorithm of fourth order and the evolution of σ z (t) is recorded as a function of time. A sufficiently small time step is employed (100 points per oscillation period) and the duration of the simulation is chosen to ensure a proper sampling of the patterns appearing in the subsequent Fourier Transform which is realized on the temporal traces. In order to account for the absence of synchronization of our measurement protocol with the driving force, we realize numerical simulations of σ z (t) for varying time delays (4 different values of φ are used, sampling 2π) and average the obtained results. This procedure has been used to simulate Figure 4d and Supplementary Figure 4. In our experiment the intra-multiplicity transition at δω 0 [Ω d ]/2 can be seen, traducing the presence of a slight detuning between the MW field and the spin qubit resonance, which can change due to thermal drifts. Each measurement of Rabi oscillations is typically accumulated for one hour and paused every 5 minutes approximately to realize an ESR measurement. This permits to adjust the microwave frequency to the qubit resonance, which slightly shifts due to spatial drifts in the strong magnetic gradient. A 50 kHz detuning is an average value of the qubit resonance shit measured between two relocking steps of the microwave frequency onto the qubit resonance. It is a good estimation of the average detuning one can expect in our measurements. It has thus been taken into account in our numerical simulations.

Interpretation of the Mollow triplet in terms of a double dressing
We formalize here the two-tone interaction (MW and phonon fields) with the NV spin qubit in terms of a double dressing of the spin. We begin by considering the effect of the quasi-resonant MW field on the spin enabling a first dressing. Then, we investigate the role of the mechanical phonon field which is resonantly coupled to the MW dressed qubit. We explain the triple peak features and anti-crossing in the FFT spectra. We focus here on the case of a single mechanical mode and investigate its contribution to the phononic Mollow triplet. The generalization to two eigenmodes, driven out of resonance is straightforward and detailed in the manuscript. The total Hamiltonian is with where we have included the interaction of the microwave and phonon fields (described by theâ,â † andb,b † annihilation and creation operators respectively) on the spin qubit, using the matrices Dressing the spin qubit with the MW field We first consider the action of H MW , which does not affect the phonon field. We use the rotating wave approximation, which corresponds to dropping the off-resonant terms(b † σ + andbσ − ) in H MW . This is justified for intense coherent MW fields with mean phonon numbersN as long as Ω v The coupling term of the so called Jaynes-Cummings model [5] corresponds to interactions where a spin excitation is traded for a MW field excitation. It only couples between sets of states {|0, N , |1, N − 1 }. We therefore introduce the multiplicities E(0) ≡ {|0, 0 } and for N ≥ 1 The Hamiltonian describing the spin-MW subsystem in the RWA approximation does not couple different multiplicities and can then be written as an infinite tensor product with H (0) ≡ 0. Introducing the MW detuning δ ≡ ω − ω 0 , the Hamiltonian H (N ) acting in the E(N ) subspace reads Diagonalization of the Hamiltonian leads to new eigenstates |± N with energies featuring a splitting ∆ N with Using and the new eigenstates can be written: They are called the dressed states. Importantly, the interaction between MW and spin causes an energy splitting of ∆ N in E N , giving rise to the Jaynes-Cummings energy ladder with sets of rungs |− N , |+ N . If the MW field is described by a large amplitude coherent field |α = e −|α| 2 /2 ∞ n=0 α n √ n! |n , the photon distribution is given by P (N ) = e −|α| 2 α 2N N ! , which is centered aroundN = |α| 2 with width ∆N =

√N
. The variation of the rung spacing ∆ N does not vary significantly over the spreading of the photon distribution ifN 1. Then, for a large coherent drive strength the dominant splitting can be described by the classical field amplitude, Ω v R

√N
= Ω R , so that ∆ N ≈ Ω 2 R + δ 2 for N −N ∆N . Since the MW dressing does not alter the phonon state, the generic eigenstates at that step are |± N , M .
Before investigating the role of the phonon field, it is important to calculate the effect of the σ z operator, describing the parametric interaction, on the dressed states. Within one multiplicity E(N ), it can be expressed in the basis The σ z operator couples the MW dressed eigenstates: the phonon field described by H int will then be able to induce transitions among the MW dressed states if its frequency Ω m matches the splitting ∆ N . Notice that σ z does not couple dressed states belonging to different multiplicities ( ± N | σ z |± N = 0 if N = N ).

Dressing the MW dressed states with the phonon field
In the picture of dressed states the role of longitudinal and transversal magnetic field components are interchanged: for near resonant pumping, the dressed spin is polarized along an axis in the equatorial plane of the Bloch sphere. This is a consequence to having rotated the eigenvectors onto the equator of the original Bloch sphere, see equation (49). Hence, when a mechanical oscillation moves the NV in space, the magnetic field experienced by the spin qubit is modulated on timescales on the order of the oscillation period. Both parallel and transverse components of the B field (with respect to the NV quantification axis, z) are modulated by the nanomotion but we only restrict our study to the parametric coupling case. As we will see, our measurement is based on a spectroscopy of the dressed qubit, and the transition frequencies are not affected at first order by perpendicular coupling terms (σ x (â +â † )). However these terms will have to be taken into account when inspecting the coherence of the dressed qubit (which is beyond the scope of that paper). Thus the mechanical motion will cause a modulation of the magnetic field along the z-direction experienced by the spin oscillating at frequency Ω m /2π. Thus, it can induce transitions between the dressed states if the field is resonant with their energy splitting ∆ N = Ω R . We then lead to develop a second dressing of the dressed qubit with the phonon field. The hybrid interaction Hamiltonian, does not couple multiplicities with different photon number N but permits raising or lowering phonon number M by one unit. Its effect can thus be studied between dressed states combined with adjacent phonon number states (|± N , M ). The following table represents the matrix elements of the interaction Hamiltonian (within an overall g z factor) Having realized that the action of σ z can induce state flips in the dressed state basis (contribution of non-diagonal terms in (50), it signifies that the phonon field can now dress the MW dressed basis |± N . It is natural to inspect its effect in term of a dressing of the dressed states with the phonon field and work on the following set of phononic multiplicities: The role of the interaction Hamiltonian is further investigated in Supplementary Figure 3. In addition to the intramultiplicity coupling terms, it also highlights that the interaction Hamiltonian moreover couples different multiplicities (F(N, M ) ↔ F(N, M ± 1)). However, the inter-multiplicity transitions are due to fast rotating coherences and we can neglect their contribution by means of a second rotating wave approximation concerning the phonon field, oscillating at Ω m /2π. Again such an approximation is valid as long as g z √ M Ω R . As explained in Supplementary Figure 3, we are then only left with the intra-multiplicity transitions, |+ N , M ⇔ |− N , M + 1 , which preserve the global dressed state and phononic excitation number: Then, the interaction Hamiltonian restricted to each MW dressed multiplicity can be simplified to Its action in the F(N, M ) basis, having uncoupled energies where an energy offset of (N ω − δ/2 − ∆ N /2 + M Ω m ) has been subtracted, can be written as where g N,M is defined as The new eigenenergies E |± N,M and eigenstates |± N,M can thus be described by where For a large MW drive, Similarly, for large oscillation amplitude √ M g z → δω 0 , δω 0 being the amplitude of the parametric energy modulation caused by the oscillation in the magnetic field gradient. The splitting becomes In case of resonant pumping δ = 0, it simplifies to: We note that the mechanical contribution is fully encoded in the Ω m , δω 0 terms. Its generalization to a mechanical oscillator driven out of resonance can simply be obtained by replacing Ω m by the classical driving frequency Ω d and δω 0 by the corresponding parametric coupling strength δω 0 [Ω d ] whose derivation has been explained in the main text and which depends on the spin qubit trajectory in space. Also the extension to a driven multimode resonator is described in the main text.

Frequencies of Rabi oscillations
The frequencies present in the Rabi oscillations correspond to the energy differences between the doubly dressed states of neighbouring multiplicities coupled by the σ z operator.
The frequencies indicated above indeed well reproduce the observed components in the Fourier spectra of Rabi oscillations. We point out that the observation of the intramultiplicity transition at frequency ∆ N,M is only possible when the MW detuning is non-zero (δ = 0). Dependency on Ω R and δω 0 is verified in Supplementary Figure 4.

Interpretation of Rabi oscillation spectra in terms of dressed spin fluorescence
In this section we further analyze the experimental observations. It is important to recall that at resonance (δ = 0) the MW dressed qubit is aligned along the x-axis of the Bloch sphere representing the original spin qubit. As we have explained, the roles of the σ x and σ z operator for the dressed spin are interchanged and the mechanical motion along the z-axis resonantly couples to the dressed states. The MW dressed qubit therefore experiences irradiation from the coherent phonon field oscillating at Ω m . Here, we describe the evolution of the dressed spin by a quantum master equation. We introduce Bloch equations for the MW dressed qubit coupled to the phonon field. This permits justifying the above mentionned rotation of the point of view in the Bloch sphere operated by the MW dressing. This also allows analyzing the spin locking signatures: in particular, we here find the dependence of the weights of the Mollow triplet peaks on the oscillation amplitude and phase. Finally, this theoretical derivation permits to justify that the measurement of the Rabi oscillation of the bare qubit permits to measure the temporal evolution of the dipole of the MW dressed qubit. Fourier Transforms of such traces thus correspond precisely to the scenario of observing a Mollow triplet: when the fluorescence of a two level system is measured in presence of an intense pumping field which resonantly drives the TLS.
The reduced density matrix for the dressed qubit Experimentally, we measure the temporal evolution of σ z : (65) Considering an example for A: we see that all parts of the sum with N = M vanish. Thus, where ρ (N ) is the density matrix reduced to the N th multiplicity, with ρ N ++ = + N | ρ |+ N . This leads to: which reduces in case of resonant pumping to: Here we have defined the reduced populations and coherences: which illustrates that we can introduce an effective two level system, which entirely describes the MW dressed qubit in case of large MW strength. The full derivation is given in [6]. Equation (68) signifies that Rabi oscillations can be understood as a measurement of the temporal evolution of the dipole or coherence (∝ σ +− + σ −+ ) of the MW dressed qubit. The spectrum of the Rabi oscillation thus represents the equivalent of a fluorescence spectrum of the MW dressed qubit. Note that in contrast to the original Mollow experiment, the Mollow triplet is not observed in a measurement of the side fluorescence but directly in the time resolved evolution of the dressed qubit dipole. We can therefore expect to be able to interpret spectra of Rabi oscillations in terms of phenomena similar to spectra of photons emitted by a single atom irradiated by a light field.

Master equation for the dressed qubit
The dressed states {|± N } forms a base of the spin MW subsystem. We describe here the dynamics of the pseudo qubit in terms of Bloch equations in presence of a parametric modulation originating from mechanical motion. The global density matrix ρ follows the evolution equatioṅ with If a large enough MW field is applied (N ∆N 1), which is the case in most experimental conditions when operating with a NV spin qubit, the variation of θ N over the dominantly populated multiplicities is negligible. Following the approximation explained previously, we can restrict our analysis to one single multiplicity with N ≈N .
For a coherent mechanical oscillation, generating a coherent parametric modulation strength δω 0 , we can rewrite g z a † + a σ z as δω 0 cos(Ω m t + α)σ z (where α is now the initial phase of the RF field). In that situation, we havė Bloch equations describing the dressed qubit evolution The damping mechanisms acting on the dressed qubit are not the same as the one involved in the original qubit, see for example [7]. The changes in the dressed qubit frequency and its orientation on the Bloch sphere make the dressed qubit sensitive to environmental noises at different frequencies and along different orientations. Here, we will introduce the decoherence rate Γ spin of the dressed qubit and and simplify the dressed state energy decay rate to 2Γ spin , which is assuming that T 1 processes of the dressed spin qubit are negligible compared to T 2 processes. We thus arrive at a set of Bloch equations for the dressed stateṡ which can be summed over N in order to yield the Bloch equations of the pseudo spin: These equations are commonly solved by means of a rotating wave approximation. In that case we express (73) in terms of the coherences in the frame rotating at the mechanical oscillator frequency, defining for which the Bloch equations in the rotating frame arė We then define the Bloch vector b α as: so that which is the quantity measured through the Rabi oscillations of the original qubit (see equation (68)). The equation of evolution of the Bloch vector is then Notice that the dependence on the initial phase α is hidden in the initial conditions at t = 0. Since in our protocol, the spin is always prepared in the ground state |0 by optical pumping, we have and u α (t = 0) = 1 2 w α (t = 0) = 0.
We will restrict the following considerations to the resonant condition, where Ω m = Ω R . We then look for solutions like and find solutions for λ by solving the eigenvalue problem det(M − λI) = 0. The eigenvalues are for which the Bloch vector evolves as We note that for large enough amplitudes δω 0 , the character of the evolution changes from an exponential decay to exponentially decaying oscillations. When δω 0 > Γ spin , The measured Rabi oscillation can be described by: where it is clear that u α is the component oscillating in phase with the mechanical motion and v α in quadrature: which we can reexpress as: With our definition of the σ z operator, the population P α 0 = 1 − σ α z (t) is evolving as P α 0 (t) = 1 2 1 + e −Γspint cos α cos (Ω m t + α) Comparison with the experimental observations Experimentally, our measurement protocols are not synchronized to the oscillatory phase, so that we have to average over α in order to recover the observed spin dynamics: For sufficiently large amplitudes, δω 0 > Γ spin , when the Mollow triplet appears, the expression simplifies to P 0 (t) = 1 2 + 1 4 e −Γspint cos (Ω m t) The three peaks of the Mollow triplet are present in the above equation, as well as their weights, widths and relative phases (which is preserved despite averaging over the mechanical phase α). Notice that the constant term reflects the mean value of Rabi oscillations: P 0 (t → ∞) = P 1 (t → ∞) = 1 2 . By transforming expression (89) into its Fourier space we can write down the analytic form of the phonon Mollow triplet: Equation (90) fully captures amplitude, frequency and relative phase of each component of the Mollow triplet at the spin locking condition. We note that the frequencies of each peak are in agreement with the expressions obtained from the double dressing approach. Expression (90) was used to fit the experimentally measured phononic Mollow triplets with adjustable weights.
Weak and strong coupling of the parametric interaction Expression (90) describes Rabi spectra taken at the resonance or spin locking condition Ω m − Ω R = 0. This corresponds to the experiment shown in Figure 3 of the article in which the oscillation strength δω 0 is gradually increased. The emergence of the triplet structure with augmenting oscillation amplitude is visible. Since the separation of the sidebands is directly related to the coupling strength, we expect that the individual components are distinguishable once the parametric modulation strength is larger than the peak linewidth: δω 0 Γ spin . Such a condition is usually called a strong coupling condition, which applies here to the resonant interaction of the phonon field with respect to the MW dressed spin. It is possible to obtain a spatial interpretation of this criterion: it suggests that the oscillation amplitude projected along the parametric coupling vector e λ is larger than the Mollow length (δr Mollow ≡ Γ spin /|∇ω 0 |, see text): δr · e λ > δr Mollow (91)

Sideband asymmetry
The measured Mollow triplets present a systematic asymmetry in the sideband heights, which is well reproduced by numerical simulations (see Figure 4d and Supplementary Figure 4). This asymmetry is absent in the numerical simulations when the rotating wave approximation is realized, as can be seen in the analytical expressions exposed above. It is a consequence of the Bloch-Siegert mechanism, a "phononic light shift" due fast counter rotating terms in the dressed qubit-phonon interaction, which renormalizes the pseudo qubit energy splitting Ω R when it is strongly driven by the phonon field [8,9]. In the experiment the observed height asymmetry may also occur due to other experimental imperfections. For example, the relative sideband amplitude depends strongly on the detuning (Ω d − Ω d ) between the mechanical drive and Rabi frequencies as evidenced in the article figure 3d and in the σ z matrix elements in the doubly dressed base. In the experiment presented by figure 4a, the Rabi frequency has been adjusted each time to the mechanical drive frequency by tuning the microwave power. Any experimental errors originating from the calibration of the power dependence of the Rabi frequency could lead to a small but finite and uncontrolled detuning between Ω R and Ω d which could cause a sideband height asymmetry. However the observed sideband amplitude ratio can be as large as 50%, which cannot be fully explained by this simple artefact since it would require a too large detuning. Indeed our experimental tuning procedure doesn't take into account the Bloch-Siegert shift of the dressed spin qubit frequency, as observed in [8] with a RF driving field. The Bloch-Siegert shift starts to play a role as soon as the coupling strength becomes comparable to the spin energy splitting. In our system, coupling strengths on the order of the dressed qubit energy (δω 0 ≈ Ω R ) are experimentally obtained. It can be shown, see [10] that the maximum in the magnitude of the central peak is displaced towards larger microwave powers by Ω BS R = Ω m + δω 2 0 16Ωm which can lead to a detuning as large as 750 kHz for a coupling strength of 8 MHz. As a consequence, this creates a sideband assymetry that increases with the coupling strength and can reach up to 50% in our case. When adding an inherent experimental noise floor to these theoretical spectra, the upper sideband will be vanishing at large separation, i.e. around the upper mechanical resonance. It is important to mention that the correction due to Bloch-Siegert effects on the measured parametric modulation strength δω 0 occurs at the second order in δω 0 /Ω m since we lie (see Figure 3d) in the parabolic section of the Mollow sidebands (numerical simulations put an upper bound of 8% on the relative error made). As a consequence, it is justified to neglect the Bloch-Siegert shift on the estimation of δω 0 , even if it can play an important role on the sidebands amplitude asymmetry. In the experiment, we have chosen the experimental procedure consisting in tuning the Rabi frequency to the mechanical drive frequency at each step, without compensating from the potential Bloch-Siegert shift. This experiment is aiming at measuring experimentally the coupling strength δω 0 for each mechanical drive frequency. In that sense the inherent Bloch-Siegert shift of the Rabi frequency cannot be anticipated in advance without knowing the actual value of δω 0 .