Intermodal Brillouin scattering in solid-core photonic crystal fibers

Paulo F. Jarschel,1, 2, a) Erick Lamilla,1, 2 Yovanny A. V. Espinel,1, 2 Ivan Aldaya,3, 1 Julian L. Pita,1, 2, 4 Andres Gil-Molina,2, 4 Gustavo S. Wiederhecker,1, 2 and Paulo Dainese1, 2, 5 1)Photonics Research Center, University of Campinas, Campinas 13083-859, SP, Brazil 2)Gleb Wataghin Physics Institute, University of Campinas, Campinas 13083-859, SP, Brazil 3)Campus of São João da Boa Vista, State University of São Paulo, São João da Boa Vista 13876-750, Brazil 4)School of Electrical Engineering, University of Campinas, SP, Campinas 13083-852, Brazil 5)Corning Research & Development Corporation, One Science Drive, Corning, New York 14830, USA


I. INTRODUCTION:
Applications exploring optical waveguides and cavities supporting multiple spatial modes have greatly expanded in recent years. Albeit previously considered as an impairment for optical communications, multimode systems regained attention as the basis of mode division multiplexing (MDM) 1,2 , and can substantially enhance the capabilities in many other applications such as sensing 3-6 , particle manipulation [7][8][9] , and nonlinear optical devices as for example in frequency comb generation in multimode ring resonators 10 and nonreciprocal devices based on Brillouin scattering in multimode waveguides 11,12 . Excitation of different spatial modes can be performed using various approaches such as electronicallyaddressable spatial light modulators (SLMs) [13][14][15][16] , photonic lanterns 17,18 , and integrated mode combiners/multiplexers 19 . However, interaction between propagating modes is more difficult to achieve, despite being a crucial functionality for alloptical systems, as in switching, mode conversion, and optical isolators. Nonlinear optical effects offer a path to enable and control intermodal interactions, for example using Kerrinduced long period gratings to perform all-optical mode conversion 20 . Stimulated Brillouin scattering (SBS) is another nonlinear mechanism that can be explored for intermodal interaction, with unique properties [21][22][23][24][25] . Besides enabling direct power exchange between different spatial modes, SBS can be explored as a mode-selective isolator or mode-selective variable optical attenuator (VOA) 26,27 , and has recently been explored as the basis for non-reciprocal devices 11,12 .
Interaction through Brillouin scattering can occur between modes that are co-propagating or counter-propagating, respectively referred as forward Brillouin scattering (FBS) or backward Brillouin scattering (BBS) 21 . In integrated silicon waveguides, high intermodal Brillouin gain has been recently reported 28 . In fibers, the first demonstration of stimulated ina) Electronic mail: jarschel@ifi.unicamp.br Paulo F. Jarschel and Erick Lamilla contributed equally to this work.
termodal FBS was obtained in all-solid fibers 22,29 , with relatively low gain. Efficient experimental demonstrations were only obtained in BBS configuration, using few-mode fibers (FMFs) [23][24][25] . As an alternative, solid-core PCFs offer enhanced optomechanical interactions due to the greater flexibility to engineer both optical and mechanical modal properties [30][31][32][33] , and is the basis for the results presented here. In this paper, we provide a comprehensive experimental and theoretical analysis of intermodal FBS in a PCF. By exploring different supported mechanical modes, we show that both polarization-independent and polarization-dependent interaction can be implemented, and demonstrate forward intermodal Brillouin gain comparable to backward Brillouin-based mode conversion in FMFs. Finally, we investigate the fundamental limitations to conversion efficiency imposed by optomechanical coupling and by different forms of mechanical dissipation.

II. PHASE-MATCHING AND SYMMETRY CONSIDERATIONS
In intermodal FBS, depicted in Fig. 1(a), optical forces generated by the beating of two optical modes selectively excite mechanical modes in the fiber, which in turn induce power exchange between the optical modes as they propagate. In particular, we investigate the interaction between the fundamental HE 11 mode and one of the high-order HE 21 modes, employed as pump and probe (Stokes line), respectively. Efficient interaction occurs when both energy and momentum are conserved, i.e. ω p = ω s + Ω and β p = β s + q, where Ω, ω p , and ω s are respectively the mechanical, pump, and Stokes angular frequencies, and q, β p , and β s are the corresponding propagation constants ( Fig. 1(b)). In practice, these conditions yield a process that is highly frequency selective, as the pump-Stokes frequency detuning must match the frequency of a given mechanical mode within the mechanical resonance linewidth, typically in the order of MHz. This property enables Brillouin-based devices to be highly wavelength selective. Another necessary condition is imposed by spatial symmetry selection rules. Following conventional mode coupling notation, the optomechanical coupling coefficient is propor- tional to the spatial overlap of the interacting fields 21 : where S is the fiber cross-section, E * p and E s are the pump and Stokes electric field profiles, and δε is the tensor of the permittivity perturbation induced by the mechanical mode. In our particular case, E p and E S exhibit even and odd reflection symmetry, respectively. As a result, the mechanical modes must induce a perturbation with odd symmetry to yield a nonzero overlap integral. This is satisfied for example by a flexural mechanical mode, as illustrated in Fig. 1(c). In the particular case of PCFs, several mechanical modes satisfy both phase-matching and spatial overlap selection rule, opening possibilities to explore novel mode-conversion functionalities in a wide range of discrete frequencies with or without polarization selectivity.

III. EXPERIMENTAL SETUP
The setup used in our experiments is illustrated in Fig. 2. On the launch side, both pump and Stokes lasers propagate in free-space. The Stokes beam reflects off a SLM, which can dynamically transform the phase profile of an incident beam, resulting in the excitation of one of the supported high-order modes (See Supplementary Material S2.A). A half-wave plate is used to adjust the Stokes polarization to match that of one of the HE 21 modes, which are roughly linearly polarized. The pump polarization is adjusted to be either parallel or orthogonal to Stokes (referred to as x and y axis from here on). Both beams are then combined and launched into the fiber.
Low repetition rate pulses are used for both pump and Stokes (generated using external modulators) in order avoid power damages to the SLM. Such pulsed scheme also helps identifying the optical modes excited through their group delay differences (see Supplementary Material section S2.A for details). Initially, in order to optimize the SLM phase scanning Fabry-Perot interferometer. L3 and L4 are 50x microscope objective lenses. (b) Time series showing the differential delay corresponding to the HE 11 and HE 21a modes. The dashed line corresponds to the position of the differential delay for the high-order mode, obtained by finite element method simulations. Experimental (top) and simulated (bottom) beam profiles for each scenario and SEM image of the PCF under study are shown in the insets. masks and characterize the excitation of various guided optical modes, we blocked off the pump beam and used short 35 ps pulses (at a repetition rate of 155 MHz) for the Stokes signal. The output of the PCF was split to enable simultaneous imaging of the beam via an InGaAs camera, and detect the output pulses with a 20 GHz photodiode. Fig. 2(b) shows the output pulses, alongside with the output beam profiles for a blank SLM phase mask (red) and a phase profile optimized to excite one of the HE 21 modes (blue). In each case, the output trace shows a single pulse, either a fast pulse corresponding to the fundamental mode or a delayed pulse, which in this case corresponds to the HE 21a mode. The output beam shows clean and well defined profiles, corresponding to the respective simulated mode shapes. A micrograph of the PCF used in our experiments is shown in Fig. 2(b), from which we measured the following parameters: 3.35±0.05 µm solid-core diameter, 4.0±0.2 µm cladding pitch and 4.0±0.1 µm cladding hole diameter. Both the differential delay and mode shape results compare well with simulations performed in Comsol Multiphysics (dashed lines and bottom inset images in Fig. 2b), employing the actual cross-section extracted from the SEM image. Using the same short-pulse modulation, we also verified that the pump arm excites the fundamental mode with a high purity (no other modes are detected in time domain traces).
Once the mode excitation is optimized, we characterized the propagation and coupling losses for both HE 11 and HE 21a modes using the cutback method. For a fiber length of 30 m we obtained a propagation loss of 0.04 dB/m and coupling loss of 1.5 dB for HE 11 , while for HE 21a , propagation and coupling losses were determined to be 0.15 dB/m and 4.5 dB, respectively. Although it is certainly possible to optimize the modelaunching scheme to minimize coupling losses in a practical device, it is not critical here for the purpose of demonstrating proof of principle.
To perform the FBS experiments, we switched to longer 80-ns square pulses at a repetition rate of 250 kHz, and each signal was individually amplified by erbium-doped fiber amplifiers, reaching peak powers of up to 2 W while preserving a relative low average power. The Stokes laser was kept at a constant operating wavelength, while the pump frequency was swept so that different mechanical modes could be excited. The pump-Stokes frequency detuning was continuously monitored using an electrical spectrum analyzer. At the output of the fiber, the beam was split and one arm is monitored with an InGaAs camera while the other is directed to two possible detection schemes, selected by a flip mirror. In the first scheme, we perform a broadband mechanical spectroscopy by sweeping the pump-Stokes frequency detuning over a wide range. Here, a novel detection scheme using a lock-in amplifier (LIA) is implemented to detect the mechanical resonances signatures. In small gain-regime, the magnitude of the LIA signal is linearly proportional to the Brillouin gain spectrum (details regarding this approach can be found in Section S2.B of the Supplementary Material). While this technique is fast and does not require tunable narrow-band optical filters, it cannot discriminate between the Stokes gain or pump depletion, since the chopper modulation is transferred from Stokes to pump due to the Brillouin interaction along the fiber. We therefore use a second detection scheme once a certain mechanical resonance is identified in the LIA spectrum. In this second, a scanning Fabry-Perot interferometer (FPI) (7.5 MHz linewidth, 1.5 GHz free-spectral range) is employed to individually measure the pump and Stokes powers at the output of the fiber, and the energy transfer between the two signals can be quantitatively characterized.

A. Intermodal FBS spectrum
Using the LIA detection method, we first obtained the intermodal Brillouin spectrum shown in Fig. 3 for parallel and orthogonal pump-signal polarizations. Multiple mechanical resonances are observed, with two clearly dominant peaks centered at 1.30 GHz and 2.13 GHz, and weaker resonances at 90 MHz, 750 MHz, and around 1.0 GHz and 1.75 GHz. No other significant resonances were observed between 2.5 GHz and 10.0GHz. To better understand the features in the experimental spectrum and identify the mechanical modes responsible for each peak, we simulated the Brillouin gain for all modes supported by the PCF structure in the frequency range of interest. Details on the simulation can be found in Section S1 of the Supplementary Material. For the calculated Brillouin gain spectrum, the vertical axis represents Brillouin gain (normalized to the highest peak). Overall, the simulated spectrum semi-quantitatively explains the most important features observed in the experiment. Particularly, it is possible to identify the mechanical modes with highest gain and obtain the relative strengths of the strongest peaks. In addition, it gives insight on the observed polarization dependence. We now discuss these features in more detail.
In order to identify the mechanical modes related to the dominant peaks, we include simulated displacement profiles as insets in Fig. 3. These modes were calculated for a full PCF structure (obtained from the SEM profile) and compared to the corresponding modes in a suspended rod with the same core diameter. As expected, the flexural nature of these modes can be clearly observed, and moreover, the field profiles remarkably resemble those in a simple rod: all peaks correlate to flexural modes, being A the fundamental, and B, C, and D higher-order modes, with displacement nodes within the core region. Another important aspect is the 2-fold and 3-fold symmetry of these modes, which is critical to understand the polarization dependence observed both in the experiments and simulations. Interestingly, the strongest mechanical resonances (peak B at 1.30 GHz and C/D at 2.13 GHz) present significantly different polarization dependency. On one hand, at 2.13 GHz we observe strong power exchange for both parallel and orthogonal polarization, while at 1.30 GHz only when pump and Stokes have orthogonal polarization. In other words, by selecting the frequency detuning, it is possible to create a device that is either dependent or virtually independent of the polarization.
The physical mechanism behind this observation is quite unique, and not expected based on the simple rod analogy (where the peak at 1.3 GHz is polarization independent, contrary to our observation in a PCF). In a rod, there are two degenerate modes with a 6-fold symmetry displacement profile at the 1.3 GHz region, rotated by 30 degrees relative to each other (one orientation is shown as mode B in Fig. 3). For parallel or orthogonal pump-probe polarization, only one or the other of these two orientations is excited, with however identical Brillouin gain (explaining why in a rod this peak is polarization independent). In a PCF, the picture changes. Even though both orientations of this mechanical mode are still present, one of them couples with the photonic crystal cladding much more strongly than the other. This hybridization with the cladding modes simply means that the mechanical energy for that particular orientation is no longer concentrated in the core, which leads to reduced its overlap with the optical mode (thus lower Brillouin gain). The orientation that hybridizes with the cladding is exactly the one that would in principle couple orthogonal pump-probe polarizations, but, as observed experimentally and confirmed in the simulation, it does so very weakly. Physically, it is quite clear why one orientation couples to the cladding more strongly than the other. The 6-fold symmetry of these two modes coincides with the 6-fold symmetry of the photonic crystal cladding. This means that the mode whose maximum displacement lobes coincide with the glass webs hybridizes strongly, while the other, whose displacement nodes are aligned to the glass webs, does not.
We now turn to the mechanical modes giving rise to the polarization independent peak at 2.13 GHz. At about this same frequency, a rod supports two degenerate modes with a 2fold symmetry, one rotated by 90 degrees relative to the other. Again, one mode couples parallel polarization, and the other orthogonal. Differently than the previous 6-fold mechanical modes, none of the 2-fold orientations hybridizes strongly with the cladding. Physically, the mismatch in symmetry (2fold vs. 6 fold) reduces the hybridization of the core modes with the photonic crystal cladding, and one can see in Fig. 3, that the mechanical energy is mostly concentrated in the core for both modes C and D. As a result, in the PCF, the orientation in C strongly couples parallel polarization, while the rotated mode in D couples orthogonal polarization, resembling the behaviour in a rod. This is also confirmed by numerical simulations (note however that these modes are no longer degenerate, due to slight mechanical distortion in the structure).
It is worth pointing out that in some of the peaks, one can observe in both experiment and simulated spectra that their shapes deviate from a typical Lorentzian. This is because these peaks are formed from a cluster of several mechanical modes around their central frequency, due to coupling between core and cladding vibrations. Obviously, the pulsed excitation character convolved with the averaging in the LIA signal (which integrates the signal as the pump frequency is swept) effectively broadens the observed spectrum. Having said that, convolving the simulated spectrum with a filter equivalent to the experimental averaging (∼ 15 MHz bandwidth), is not sufficient to fully explain the experiment (see for example a single peak at 2.13 GHz in the experiment while the simulation even after convolution shows two separate peaks). We point out that other possible reasons for this disagreement might be a mismatch between the actual and simulated fiber structure or fluctuations along the fiber length. Another aspect that upon careful observation might seem surprising is that the fundamental flexural mode at 90 MHz is not the dominant Brillouin peak. In a rod, one can show that the equivalent mode indeed exhibits the strongest Brillouin gain, and that is not the case in the PCF studied here, as observed experimentally and confirmed by the simulations. We will return to this point when we investigate in more detail the physical process dominating the strength of the optomechanical coupling as well as the damping mechanisms of the mechanical modes.

B. Intermodal FBS energy transfer
To quantitatively analyze the intermodal energy transfer at the most prominent Brillouin peaks, we performed a fine frequency scan and used the FPI detection scheme that alllows discrimination between pump and Stokes signals. Initially, both have a fixed frequency separated by 2.13 GHz, corresponding to one of the Brillouin peaks in Fig. 3, and the FPI is scanned to measure the output power of each signal separately. The results are presented in Fig. 4(a) for three conditions: (i) both pump and Stokes on, (ii) only Stokes on, and (iii) only pump on. We used 1 W of input peak power for the Stokes and 2 W for the pump. Clearly, when both signals are present, a strong depletion of the pump is accompanied by a Stokes gain, as a result of the intermodal FBS process. A convenient way to quantify the energy transfer is the on/off gain, defined as the ratio of the Stokes output power levels with pump on and off. Analogously, the pump depletion is the ratio of the pump output power measured when the Stokes signal is turned on and off. The curves from Figs. 4(b, c) show the measured on/off Stokes gain and pump depletion for a narrow sweep around the 2.13 GHz frequency detuning. As the frequency separation approaches the Brillouin resonance, gain and depletion are maximized, reaching about 3 dB peak values. The solid lines represent numerical fitting that resulted in a Lorentzian curve with a 45 MHz FWHM. Performing a deconvolution of the signal with an effective Lorentzian accounting for the pulsed source and the FPI linewidth, we estimate a Brillouin linewidth of 42 MHz. To extract the gain coefficient, we can numerically solve the coupled equations in FBS process and fit the experimental data (See section S2.C of the Supplementary Material for details). We applied this procedure for three relevant Brillouin peaks by setting the corresponding pump-Stokes frequency detuning, for parallel and orthogonal pump/Stokes polarization conditions. The results are summarized in Table I.

V. DISCUSSION AND CONCLUSIONS
From the experimental results listed on Table I, the highest gain coefficients are 21.4 W −1 km −1 for a detuning of 2.13 GHz, and 19.9 W −1 km −1 for 1.31 GHz. In fact, this gain coefficient obtained in forward configuration in PCF is on the same order as observed in backward intermodal Brillouin scattering in FMFs 24 . This is quite remarkable because again, in all-solid fibers, forward Brillouin scattering involve mechanical modes that are distributed throughout the cladding with very little overlap with the optical modes. In PCF, transverse mechanical confinement enables high gain in forward configuration. Having said that, a more complete picture is provided here to better understand the physical limits and possible enhancement directions.
Fundamentally, the maximum Brillouin gain for a given resonance depends on the strength of the optomechanical coupling and on the linewidth of that particular resonance. It is therefore useful to separate the discussion of the optomechanical coupling strength from the damping mechanisms of the mechanical modes. In Fig. 5, for all mechanical modes we plot separately g/Q, which represents the optomechanical coupling strength, and the inverse of Q, which represents the damping strength. It is more convenient to look at Q −1 as we can directly add contributions from different damping mechanisms (i.e. Q −1 = ∑ i Q −1 i ). The number of mechanical modes in the complete structure is quite large, and thus leads to very dense curves of the calculated coefficients. To help visualize the general trend, the results shown in this figure represent the envelop curves (highest g/Q and lowest 1/Q), while the raw data for all modes can be found in the supplementary materials.
From the g/Q curves in Fig.5(a), we can see that the lower frequency flexural modes indeed exhibit larger coupling strength (0.3 W −1 km −1 for 90 MHz) than higher frequencies modes (0.04 W −1 km −1 for 2.13 GHz), as expected from the suspended rod model. The results in Fig. 5(a) include contributions from moving-boundary and elasto-optic effects, which can reinforce or counter-act each other 21,34 . We can see that, for the fundamental flexural modes, large optomechanical coupling is due to the moving-boundary contribution at the glass-air interface, a mechanism that is irrelevant in solid fibers given their small core-clad index constrast. For higher-order flexural modes, the elasto-optic mechanism dominates with overall lower g/Q values. Despite stronger optomechanical coupling, our experiment shows that these low frequency flexural modes do not dominate the Brillouin spectrum. This can be explained by evaluating the damping mechanisms shown in Fig. 5(b). Here, the lower frequency modes suffer the strongest damping (poorer quality factors), overriding their high optomechanical coupling. In this analysis, we considered three loss mechanisms for the mechanical modes: viscosity, leakage, and squeezed film air damping. Viscous damping scales quadratically with frequency (Q ∝ 1/ω m ), which means that the quality factor would be inversely proportional to frequency if this was the dominant loss, again favoring stronger Brillouin interaction at lower frequencies. However, it is clear from our modeling that viscosity is not the limiting factor, as shown in Fig. 5(b). Leakage through the cladding varies from mode to mode and the profiles in Fig. 3 illustrate that some modes are more confined to the core than others. In the modeling, this was accounted for using a mechanical perfect matched layer at the outer silica cladding surface (see section S1.B of the Supplementary Material for details). Leakage is dominant above 400 MHz, as shown in Fig. 5(b).
The main limitation for the low-frequency flexural modes arises from the effect of squeezed film air damping phenomenon, by which a micro-vibrating membrane transfers part of its energy to surrounding gas molecules 35 . Since there are many membranes vibrating in phase inside the PCF structure, this mechanism can be significant 36 . We used Bao's model 35 to estimate this contribution. From the results in Fig.  5(b), it is clear that squeezed film air damping is strongest for low order flexural modes, limiting their quality factors. For the 90 MHz peak, this mechanism lowers the quality factor from 1800 to 33, or in terms of linewidth broadening, it is increased from tens of kHz up to 3 MHz, consistent with the experimental value of approximately 5 MHz. At high frequencies, the calculated broadening due to squeezed film damping is negligible, and the obtained value of 35 MHz for the 2.13 GHz peak is consistent with the experimental value of 42 MHz. Another potential linewidth broadening factor not considered in the analysis here is the geometry non-uniformity along the fiber length. According to our simulations, a 5% geometry scaling results in 18% and 8% Brillouin frequency variation for 90 MHz and 2.13 GHz peaks, respectively. In other words, this form of inhomogeneous broadening due to geometry fluctuation might affect low frequency flexural modes more strongly than higher ones. When all damping effects are considered, the simulated gain values drop to 6.9 W −1 km −1 , 20.3 W −1 km −1 and 19.6 W −1 km −1 , consistent with the experimental values from Table I. From a practical point of view, one could envision engineering structures that better confine mechanical modes to the core (i.e. reduce leakage), and rely on larger structures suspended-core fibers 37 to minimize squeezed film damping. In this ultimate scenario, the Brillouin gain for the 90 MHz flexural mode would be limited by viscosity as a fundamental mechanism, and could potentially reach values higher than 1000 W −1 km −1 .
In conclusion, we demonstrated intermodal forward Brillouin scattering in PCF between the HE 11 and HE 21 modes over a wide range of frequencies (90 MHz -2.5 GHz), with a maximum gain coefficient of 21.4 W −1 km −1 . Symmetry arguments support the observation of polarization dependent and polarization independent intermodal Brillouin interaction, and different damping mechanisms dictate the relative strength of FBS involving mechanical modes at low and high frequencies. This work opens the path to future engineering of PCF structures to enhance the interaction of optical modes of high order through the control of flexural mechanical modes and their dissipation mechanisms, leading ultimately to novel physical phenomena and highly efficient devices.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. In this section, we detail the most important aspects of modeling forward stimulated Brillouin scattering (FSBS) in a solid core solid-core photonic crystal fiber (SC-PCF). These results were used to generate Fig. 3 of the main text. In order to find both optical and mechanical modes, the full-vectorial wave equations were independently solved by the finite elements method (FEM). Specifically, we employed the Electromagnetic Waves and Solid Mechanics modules of COMSOL Multiphysics ® software (version 5.4). The computational domains used in our simulations are shown in fig. S1. All simulations were carried out at a fixed optical wavelength (λ = 1.55 µm), and we used the material properties listed in Table I TABLE II. Simulated group index, and differential group delays for various modes (relative to the HE 11x mode).

A. Optical modes
All simulations shown in the main text were carried out using the realistic PCF geometry shown in fig. S1, which was obtained from the SEM cross-section image. However, in order to precisely identify all the optical modes observed in the experiment, we compared the realistic geometry solutions with two idealized ones: (i) an ideal triangular-lattice PCF structure, and (ii) a stretched version of this ideal PCF to account for the asymmetry present in the actual structure (6% stretching along the x-axis). While in the ideal PCF, the higher order modes shown in fig. S2(a) display the usual donut-shaped intensity profile, this is not the case for the actual fiber modes, shown in fig. S2(b). Despite the subtle differences between the ideal and actual fiber geometries, the optical modes of the stretched-ideal geometry, shown in fig. S2(c), clearly capture most features of the actual fiber. Notably, the donut-like shape of the higher-order modes is lost as a result of lifting the degeneracy of the two polarization states of the HE 21 modes. Therefore, as far as the optical mode is concerned, the impact of geometry distortion is rather predictable. Mode shapes provide a first guidance to help identify which modes are excited experimentally, along with modal group delay. In fig. S2(d) we show the effective index (n eff ) for various modes, from which we can calculate the corresponding group index (n g = n eff − λ ∂ λ n eff ) and the differential group delay τ. The results for the realistic structure are shown in Table II, where the differential group delays are calculated relative to the fundamental mode HE 11x (τ = ∆n g L/c, using the fiber length, L = 30 m). In section S2 A it is shown how these distinct group-indexes and mode profiles enables the identification of the optical modes excited in our experiment.

B. Mechanical modes
Similarly to the optical modes, although the actual PCF geometry of fig. S1 was used in all simulations, exploring an idealized structure is helpful in identifying the mechanical mode families observed in the experiment. Due to specular symmetry with respect to both Cartesian axes, only a quarter of the full structure is used, which not only spares computational resources but also allows using boundary conditions that filters specific mode families. We applied Dirichlet Boundary Conditions (DBCs) in the last ring, u = 0, and symmetric (S, u ·n = 0), or antisymmetric (AS, u ·t = 0) boundary conditions along the x-axis (C 2 ) and the yaxis (C 1 ) to model the deformation of the flexural modes, as shown in fig. S3(a). In fig. S3(b)  spectral density of mechanical modes, we choose to represent them as a density plot in fig. S3(b). The color scale represents the fraction of mechanical energy that is concentrated in the glass core (black 100%, white 0%), indicating which modes might interact more strongly with the optical modes. We can see that the dispersion relation curves do not vary significantly in this range. These characteristics are also reflected in the realistic PCF. In the same plot, we show the dispersion relation for the flexural modes in a circular rod waveguide with equivalent dimensions (same core radius, r = 1.7 µm). Several modes in the PCF have dispersion and profile similar to those in the rod. For example, mechanical displacement profiles at β m = 1.88 × 10 5 m −1 are illustrated in fig. S3(b) for three modes A, B and C in the ideal PCF, and closely resemble their corresponding modes in the rod. Furthermore, these modes also have their correspondents in the actual PCF. For example, fig. S3(c) shows the flexural-like modes at ∼90 MHz, ∼1.3 GHz and ∼2.1 GHz. Note that for the low frequency mode, its mechanical displacement field is distributed throughout the core and cladding, whereas the high frequency mode is mainly concentrated in the core region. These results suggest a stronger influence of low-frequency mechanical modes upon the optomechanical coupling process.

Mechanical dissipation
a. Leakage loss mechanism: Leakage loss accounts for the energy escaping towards the solid silica region surrounding the PCF structure. To quantify its contribution, we included a perfectly matched layer (PML) 1,2 wrapping the photonic crystal, as shown in Fig. S1, and we calculated the energy flowing on it. To implement a PML in COMSOL ® , the complex coordinate stretching function used is f (ρ) = λ t sρ P (1 − i) 2 , where λ t is the typical wavelength of the transverse component of the mechanical wave deformation (u), s is the scaling factor that tunes the wave amplitude arriving at the edge of the PML, and P is the scaling curvature parameter that controls how strong is the amplitude decay inside the PML. The PML parameters used were: s = 0.25 and P = 2.0. We employed λ t = 2π[(Ω m /V s ) 2 − β 2 m ] −1/2 , with V s = 3766 m/s as the shear bulk velocity in silica, and β m = 2π(n eff,p − n eff,S )/λ = 1.86 × 10 5 1/m, from momentum conservation. The effective indexes n eff,p = 1.413067 and n eff,S = 1.367107 are the calculated at λ = 1.55 µm for the pump (HE 11x ) and Stokes (HE 21a ), respectively. Finally, Ω m is a mechanical frequency initial guess near the eigenfrequencies being sought for.
b. Viscous damping Viscous damping (VD) was included in the PCF model as an additional stress term proportional to the elastic strain rate 3 , σ VD = η : ∂ S ∂t , where η is known as the fourth-rank viscosity tensor. The mechanical wave equation is then written as:∇ with c = c + iΩη and c being the stiffness tensor. The tilde in the divergence operator above indicates the coordinates transformation that arises from the PML. As fused silica is an isotropic material, the stiffness tensor has only two independent coefficients, c 11 = 78.5 GPa and c 44 = 31.2 GPa 3 , as well as the viscosity tensor, η 11 = 0.00279 Pa·s and η 44 = 0.00104 Pa·s 3,4 . c. Squeeze-air film damping The presence of air in the micro-metric channels in the PCF structure leads to a loss mechanisms known as squeeze-air film damping (SFD) 6 . The energy transfer model has been successfully applied to different micro-electro-mechanical systems (MEMS), including slab-air-substrate systems and microstructured optical fibers with doublemembrane 6,7 . According to this model, the SFD linewidth in a slab-air-substrate system is given by: in which, S, ρ, and H are the slab's length, mass density, and thickness, respectively. h 0 is the gap between the slab and the substrate, M m is the molar mass of the air, R is the ideal gas constant, T is the temperature, and p denotes the pressure. We heuristically extended this model to a bore by mapping each circular air-hole to the vibrating slab model described by eq. (S2). The slab length is considered as the hole's perimeter (S = 2πr c ), while the slab-substrate gap h 0 is mapped to an effective bore diameter, given by the area to perimeter ratio, h 0 = πr 2 c /S. Such mapping is feasible for a purely radial motion within each bore, however, in a general mechanical mode of the PCF structure, not all holes are radially deformed (e.g., motion can be tangential to the boundary). Since we are considering only the squeeze film effect, we take this into account by introducing an effective number of radially deformed holes given by: in which, u t is the deformation in the PCF cross section of the mode, obtained by solving eq. (S1), and r c is the PCF core radius. The integral is performed across all silica-air boundaries defining each air hole. With this, the SFD contribution to the linewidth of a given mechanical mode, γ SFD , can be calculated combining eq. (S2) and eq. (S3): d. Total mechanical Q-factor and linewidth We first solved eq. (S1) in order to find the Q-factor from the combined leaky and viscous damping mechanisms, Q L+VD = ℜ(Ω m )/[2ℑ(Ω m )], where ℜ(Ω m ) and ℑ(Ω m ) are respectively the real and imaginary parts of the mechanical eigenfrequency, Ω m . Then, we extracted the SFD Q-factor from eq. (S4) as, Q SFD = ℜ(Ω m )/γ SFD . Finally, we obtain the total mechanical Q-factor, Q m = Q SFD Q L+VD /(Q SFD + Q L+VD ) and γ m = ℜ(Ω m )/Q m .

C. Coupled-mode description of Brillouin interaction
Using the formalism adopted by   9 , we consider the interaction between two optical modes corresponding to pump and Stokes waves oscillating at ω p and ω s with complex amplitudesã p andã s , respectively, and K mechanical modes with complex amplitudesb 1 ,b 2 , ...,b K . Both optical and mechanical amplitudes are defined in such a way that |ã p | 2 and |ã s | 2 represent the number of photons per unit length, and |b k | 2 the number of the phonons per unit length of the k-th mechanical mode. In the presence of multiple mechanical modes, an effective complex gain coefficient can be derived. Assuming large mechanical damping (v m /γ m 1), steady-state regime (∂ t → 0), and neglect the Stokes and pump losses (α s , α p ≈ 0), the coupled equations become: where ∂ z is the partial derivative with respect to the spatial coordinate, and v p , v s , and v m,k represent the group velocity of the pump, Stokes and the k-th mechanical mode, respectively. (·) * denotes the complex conjugate operation, ∆ m,k = Ω − ℜ(Ω m,k ), in which Ω = ω p − ω s .g 0,k is the waveguide optomechanical coupling rate [ √ m/s], whose modulus is related to the Brillouin gain coefficient, G B,k 9,10 . In the undepleted regime (∂ zãp ≈ 0), the evolution of the Stokes signal is given by: Using the pump power defined as P p =hω p v p ã p 2 , and Q-factor as Q m,k = ℜ(Ω m,k )/γ m,k , we can express equation S6 in terms of an effective complex Brillouin gain coefficient G c B : Explicitly, the effective gain coefficient is can be written in terms of the moving boundary ( f mb ) and photoelastic force densities ( f pe ) 9 : The MB integral must be evaluated across the perimeters of all holes in the first PCF ring (where the optical mode is large), whereas the PE integral must account for the whole cross-section area involving mechanical and optical modes, as is the grey region in fig. S1. Finally, the Stokes power defined as P s =hω s v s |ã s | 2 is: in which P s,0 =hω s v s |ã s,0 | 2 is the initial Stokes optical power and L is the fiber length. According to eq. (S9) the small signal gain is given by the real part of the complex function G c B .

D. Numerical results
Fig . S4, we plot separately g/Q, which represents the optomechanical coupling strength, and the inverse of Q, which represents the damping strength for all mechanical modes. The dashed curves represent envelop curves (highest g/Q and lowest 1/Q) as discussed in the main manuscript.

S2. EXPERIMENTAL CONSIDERATIONS
Here we give a more in-depth description of how the distinct optical modes were excited and characterized (section S2 A) and how the lock-in amplifier was employed to identify the active mechanical modes (section S2 B). The detailed version of our experimental setup is shown in fig. S5.

A. Modal excitation and characterization
In order to excite the HE 21 a mode, we used the phase mask shown in the inset (ii) of Fig. S6. This mask was built by multiplying a blazed grating and a π-difference phase mask. By doing so, we were able to use the first-order diffraction of the grating to convert the incident Gaussian beam into the desired higher order mode with improved beam quality. The conversion efficiency of this approach, however, depends on the incident polarization, so a fiber polarization controller was included at the output of the fiber before conversion to free space. Additionally, with the aim of maximizing the matching between the output beam and the higher order mode of the fiber, the phase mask was rotated and shifted, whereas a λ /2-plate was included to minimize the polarization mismatch. It is important to note that during the realization of the experiments, the purity of the mode was periodically monitored by short pulse measurements as the ones shown in Fig. S6, and by imaging the output beam with the InGaAs CCD camera. Note that for a phase mask with a π along the horizontal axis (iii), the beam profile resembles well the HE 21 a mode, and the measured group delay is 140 ps/m, which is in good agreement to the simulated value of 126 ps/m.

B. FBS frequency identification
As indicated in fig. S5, the probe laser (Stokes) is modulated using a chopper in order to enable a lock-in amplifier (LIA) detection scheme. This set the basis to build Figure 3 in the main text. It is not obvious, however, that signal resulting from this detection scheme should be related to the Stokes probe gain when the detuning matches a mechanical resonance. For example, for lossless optical modes, the Manley-Rowe relations 9 imply that the photons lost from the pump mode are transferred to the probe wave. Since the mechanical frequency is much smaller than the optical frequency, this is equivalent to state that the pump power depletion must match the probe power gain. Since the signal generated at photodetector PD2 is proportional to the total power (pump + Stokes), no change would be observed on the LIA signal whether the lasers are tuned to a mechanical resonance or not. In other words, the net LIA signal in this lossless scenario would not be sensitive to the Brillouin gain.
To investigate the connection between the LIA signal and the actual Stokes gain in a realistic lossy fiber, considering pump (fundamental mode) and Stokes (higher-order mode) signals with loss coefficients of 0.04 dB/m and 0.15 dB/m, we solved the steady-state eq. (S5) with a square-wave modulated probe signal. Since the bandwidth of the photodetectors PD2 or PM are much smaller than the pump-probe frequency detuning, their time-dependent output signal is simply given by the sum of pump and probe powers. At the fiber input, as it would be detected by photodetector PM, these signal are depicted in fig. S7(a). The simulated output pulses are shown in fig. S7(b). At resonance, energy transfer occurs in the time intervals when the Stokes probe is on. However, the Stokes gain is no longer equal to the pump depletion, simply as a consequence of different attenuation rates (α p = α S ). In other words, power flows from a low-loss channel (fundamental mode) into a high-loss channel (HE 11a ). As a result, when the lasers are tuned to a mechanical resonance, the pump power variation (δ P p ) is clearly different from the Stokes power variation(δ P S ). To confirm the importance of the loss imbalance for the lock-in detection scheme, we show in fig. S7(c) how these power variations would depend on the detuning (∆) between the pump-Stokes beat-note (Ω = ω p − ω S ) and the mechanical frequency of a given mechanical mode, i.e., (∆ = Ω − Ω m (normalized by the mechanical linewidth, γ m ), also, in fig. S7(c) we show how different different values of pump and Stokes mode optical losses impact the LIA signal. For instance, if the two modes present exactly the same loss, the LIA signal would have a flat response. Also this model explains whether dips or peaks are to be observed in the LIA signal: if α p > α S , a peak would be observed, meaning that the Stokes gain contributes more to the signal, while if α p < α S , a dip is measured, indicating that the pump depletion is more significant at the output signal. In order to obtain the traces in Fig. 3 of the main text, the frequency detuning between the pump and Stokes was swept using the piezo-scan available in the probe-laser, so a temporal trace of the output power with well-defined dips could be obtained. It is worth noting that in order to isolate the power reduction due to intermodal power transfer, it is useful to measure the baseline in absence of FBS by blocking the pump wave. The frequency resolution of the employed method is limited by two factors. On the one hand, the use of pulses for both Stokes and pump poses a lower resolution boundary because the effective linewidth is broadened. As it is well known, longer pulses result in narrower linewidth but, unfortunately, it is not possible to arbitrarily increase the pulse length because longer pulses lead to lower peak power. However, the incident power limitation of the spatial light modulator (SLM) (200 mW) should also be considered. Therefore, 80 ns pulses were adopted in our experiments since this value simultaneously allowed a relatively high peak input power of 2 W and an acceptable frequency resolution of 12.5 MHz. On the other hand, we should consider the intrinsic resolution of the LIA, which further depends on the pump-Stokes frequency sweeping speed and LIA integration time. These parameters should be chosen to optimize the trade-off between frequency resolution, and noise bandwidth. In this way, for a fixed sweeping speed, a larger integration time leads to improved signal-to-noise ratio but poorer resolution and vice versa. To get the curves shown in Fig. 3(a) of the main manuscript, we configured the lock-in amplifier to have a sweeping speed of 300 MHz/s and an integration time of 30 ms, which results in a ∼9 MHz resolution. Therefore, to find the combined frequency resolution, we performed the convolution of a sinc function with a full-width at half-maximum (FWHM) of 12.5 MHz and a Lorentzian with a FWHM linewidth of 9 MHz, resulting in a total frequency resolution of 15.5 MHz. To obtain the Brillouin gain coefficient, we measured several on/off gain (and depletion) points as a function of both Stokes and pump input powers, and the results are shown in Fig. S8(a) and Fig. S8(b), for the 2.13 GHz peak. From the data points, we can see that the Stokes gain value reaches 2.6 dB for the Stokes and pump respective input powers of 0.95 W and 0.15 W, while the pump depletion reaches 2.8 dB for 0.95 W and 1.07 W. By introducing loss into eq. (S5), the power exchange is governed by the following coupled rate equations: dP p (z) dz = −α p P p (z) − gP p (z)P s (z) (S10) dP s (z) dz = −α s P s (z) + gP p (z)P s (z), where P p and P s correspond to the pump and Stokes powers, respectively, z is the position along the fiber, α corresponds to the attenuation coefficient for the different optical modes, and g is the Brillouin gain coefficient. Equations (S10) and (S11) were solved numerically using the measured attenuation coefficients for HE 11 (0.04 dB/km) and HE 21 (0.15 dB/km) modes, and leaving g as the fitting parameter of the Least Squares method.