|Home | About | Journals | Submit | Contact Us | Français|
The simultaneous control of optical and mechanical waves has enabled a range of fundamental and technological breakthroughs, from the demonstration of ultra-stable frequency reference devices, to the exploration of the quantum-classical boundaries in optomechanical laser-cooling experiments. More recently, such an optomechanical interaction has been observed in integrated nano-waveguides and microcavities in the Brillouin regime, where short-wavelength mechanical modes scatter light at several GHz. Here we engineer coupled optical microcavities to enable a low threshold excitation of mechanical travelling-wave modes through backward stimulated Brillouin scattering. Exploring the backward scattering we propose silicon microcavity designs based on laterally coupled single and double-layer cavities, the proposed structures enable optomechanical coupling with very high frequency modes (11 to 25GHz) and large optomechanical coupling rates (g0/2π) from 50kHz to 90kHz.
Brillouin scattering occurs due to the interaction of optical and mechanical waves and it leads to the inelastic scattering of pump photons to Doppler red-shifted (Stokes) or blue-shifted (anti-Stokes) photons. In optical waveguides and microcavities this interaction occurs due to a combination of the photo-elastic effect1, induced by strain, and moving-boundary effect caused by the mechanical mode distortion of the optical boundaries2. These two scattering processes are strongly influenced by optical and mechanical properties of the confining structure and can be tailored for various applications. For instance, the generation of anti-Stokes photons, which is accompanied by destruction of phonon quanta, can be used to cool down mechanical modes in optical optomechanical cavities3,4; whereas the generation of stokes photons, which create phonons (heating), may foster the development of high-coherence lasers, ultra-stable radio frequency (RF) synthesizers5,6,7,8, and broadband tuning of RF filters9. Such confinement-enhanced optomechanical interaction has been observed as Stimulated Raman-like10 and Brillouin scattering in a range of photonic waveguides11,12,13,14,15,16,17,18 — where both energy and momentum conservation are directly fulfilled. In microcavities, however, the short roundtrip length and narrow linewidth further constrain the conservation laws, requiring both pump and scattered waves to be resonant with the optical cavity modes in order to ensure efficient Brillouin scattering. These constraints have limited the current cavity demonstrations of Brillouin scattering either to mm-scale cavities5,19,20,21, whose optical free-spectral range matches the mechanical resonant frequency; or heavily multimode micro-cavities, where the frequency difference of the distinct transverse optical modes accidentally matches the mechanical frequency, both at the cost of reduced optomechanical coupling.
Here we explore a compound microcavity system based on silicon microdisk cavities and demonstrate its potential to drastically enhances backward Brillouin scattering (BBS) for mechanical modes at tens of GHz. The compound microcavity scheme is illustrated in Fig. 1a and can ensure a doubly-resonant condition for the pump and stokes wave, yet preserving the small footprint necessary to achieve large optomechanical coupling. The individual cavities forming the compound system could be of various types, however, they must support mechanical modes whose frequencies can match the optical frequency splitting of the coupled modes. By engineering the mechanical modes, of single-disk (sd, Fig. 1b) and double-disk (dd, Fig. 1c) optical microcavities, to avoid cancellation between the photo-elastic (pe) and moving-boundary (mb) effect22, we demonstrate that both cavity designs could be exploited as components of the compound cavity scheme, thus offering a promising route towards the demonstration of low threshold backward stimulated Brillouin lasing in a compact silicon platform.
In backward Brillouin scattering the optical pump and the scattered Stokes waves propagate in opposite directions, resulting in a large wavevector mismatch that favors the interaction between light and short-wavelength propagating mechanical modes1. In disk microcavities, the optical and mechanical modes are azimuthally traveling waves with dependence exp(±imϕ) (here m is an integer and ϕ the azimuthal angle). Hence, a pump laser exciting an optical cavity mode with frequency and azimuthal number (ωp, mp) may be scattered into another optical mode (ωs, ms) through the interaction with a mechanical mode (Ω, M), provided that both energy and momentum (phase-matching) are conserved, i.e ωs(ms)=ωp(mp)±Ω(M) and ms=mp±M. While in forward Brillouin scattering the phase-matching condition favors mechanical modes close to their cut-off condition M=0 (ms=mp), in backward Brillouin scattering (BBS) the scattered light frequency shift is proportional to the optical wavevector mismatch and can easily reach tens of GHz in solids. In order to enhance BBS, such a large frequency shift would require the pump wave to be detuned from the optical resonance by tens or even hundreds of linewidth — in a single-resonance cavity, such a large detuning would drop the benefits of the resonant cavity build-up for the pump wave. In the proposed compound cavity system, illustrated in Fig. 1a, the interaction between the optical modes (through their evanescent fields) leads to a frequency splitting that can match the mechanical mode frequency. The outcome of this scheme is illustrated in Fig. 1d with the pump wave tuned to the higher frequency coupled mode at ωp, while the scattered wave is resonant with the lower frequency coupled mode at ωs, thus ensuring a high photonic density of states (PDOS) both at the pump and scattered frequencies (see Fig. 1e).
The large azimuthal numbers involved in BBS imply that the phase-matched mechanical modes are localized near the cavity edge23,24,25, compared to low azimuthal number modes that are spread throughout the cavity; such an edge localization effectively increases the overlap between the optical and mechanical modes26. The mechanical mode induces strain and boundary deformation at the cavity edge, which efficiently Bragg scatter light and couples forward and backward propagating optical modes. The resulting energy exchange between the optical pump and Stokes wave can be modeled using coupled mode theory27,28, which leads to a set of coupled equations for the amplitudes of the pump, stokes, and mechanical waves (see Methods). We use the coupled mode theory to derive the threshold power necessary to achieve stimulated Brillouin lasing, which occurs when the Stokes photons gain, induced by the optical pump, suppresses the Stokes cavity mode loss. By assuming an undepleted pump the following expression can be derived for the input threshold power27 (see Supplementary Information),
where is the so-called single-photon cooperativity; is the vacuum optomechanical coupling rate for the compound cavity, Γ is the mechanical loss rate and are the pump and stokes detunings relative to the cavity resonance doublets (), as illustrated in Fig. 1f; κs and κp are the corresponding total (intrinsic and extrinsic) loss rates, κe is the extrinsic pump loss rate due to coupling to the driving mode of the bus waveguide.
The lowest input threshold power is reached when both pump and Stokes waves are resonant with the mode doublet, (Δs, Δp)=0. Since the Stokes photons are initially generated by spontaneous Brillouin scattering (due to thermally driven phonons), their frequencies are Doppler shifted ωs=ωp−Ω. When the pump and Stokes optical modes are split by a frequency difference J, their detuning is given by Δs=Δp+(J−Ω) (see Fig. 1f). Consequently, the threshold power scales as , thus the minimum threshold occurs when the optical splitting precisely matches the mechanical frequency, i.e. J=Ω. Such a threshold power scaling reveals the importance of the doubly-resonant condition ensured by the compound cavity scheme. For instance, the minimum threshold power achievable using a standard single-resonance cavity occurs in the so-called resolved sideband limit (Ωmκ), at an optimum pump detuning (Δp=Ω); this limit can be obtained from eq. 1: threshold assuming a large cavity separation (J→0, then κs→κp and ). Hence a single-resonance cavity has a threshold power that is (Ωm/κp)2 larger than the proposed compound cavity doubly-resonant approach. For a typical 5µm radius microdisk (Ωm/κp)2≈100, a roughly two-orders of magnitude higher threshold; where we assume an intrinsic quality factor of 2×105 (κp,s/2π≈960MHz) and Brillouin frequency Ω/2π=VlM/r≈22GHz — where M/r≈4π(neff/λ) and neff=1.73 for the transverse-magnetic (TM) optical mode phase index (λ=1.55m), Vl is the Si bulk dilatational wave velocity.
Such high mechanical frequencies at tens of GHz can also be readily matched to the optical resonance splitting accessible by employing either single or double-disk silicon cavities. For example, we show in Fig. 1g the numerically calculated frequency splitting curves (solid lines) for a compound microcavity based on silicon single disks (see Methods), demonstrating splitting at tens of GHz for a lateral gap of 100nm. This contrasts with larger and lower refractive index coupled microcavities whose frequency splitting lies in the sub-GHz range29 and could not enable resonant backward Brillouin scattering.
We demonstrate the feasibility of our compound cavity scheme by investigating two designs that can achieve high optomechanical coupling rates and mechanical frequencies at tens of GHz, the sd and dd cavities shown in Fig. 1b,c. Both designs could be readily fabricated30 and enhance the optomechanical interaction for distinct mechanical mode families. The frequency dispersion is the starting point to infer general characteristics of the phase-matched mechanical modes leading to backward Brillouin scattering. Many aspects of the mechanical modes dispersion of sd and dd cavities can be regarded as mixtures among whispering gallery mechanical modes of an infinite cylinder, and Lamb-wave modes of a free-standing silicon slab23,31. The mechanical mode dispersion of even and odd modes (with respect to z=0 plane — see Fig. 1b) of a single disk are shown in Fig. 2. The dispersion curves were calculated using an axisymmetric finite element method (see Methods). In the dd-cavities, the mechanical modes are approximately symmetric/anti-symmetric combinations of the even and odd parity sd-cavity modes, therefore, the key characteristics of both designs may be inferred by inspecting only the sd-cavity mode structure.
The mechanical modes of sd cavities that may efficiently interact with optical modes can be divided in four groups: whispering gallery (w-modes), Rayleigh (r-modes), dilatational (d-modes), and flexural modes (f-modes). Their dispersion curves are shown in Fig. 2(a–c), and signaled by markers in Fig. 2(b,c); the corresponding displacement profiles are shown in Fig. 2(d–g). The whispering gallery group (w-modes) modes are remarkably similar to the modes of an infinite cylinder, as suggested by the excellent agreement between their dispersion curves (red-dashed lines), shown in the zoomed dispersion diagram of Fig. 2b. Despite the very small thickness/radius ratio of our disk (t/r=0.05), the mechanical displacement components of w-modes are also similar to infinite cylinder modes, essentially in the radial-azimuthal (rϕ) directions. The peak displacement of increasing frequency w-modes shifts radially inward, which suggests a varying overlap with the optical mode.
The Rayleigh (r-modes), dilatational (d-modes), and flexural (f-modes) mode groups are signatures of the thin disk; their dominant radial-vertical (rz) displacement are noticeable in Fig. 2e–g. The r-mode is a singleton group and has the lowest frequency dispersion branch in Fig. 2a, which is characterized by a phase velocity lower than both the longitudinal (Vl) and transverse (Vt) bulk velocities23,31. Such a surface wave localization, noticeable in the displacement profile of Fig. 2e, compromises its overlap with the optical mode.
The dilatation (d-modes) originate from the tight vertical confinement of our thin disks. These characteristics are evidenced not only by their displacement profiles, shown in Fig. 2g, but also through their matching dispersion characteristics, highlighted in Fig. 2b (blue-dashed curve). The onset of the distinct disk mechanical mode families in Fig. 2a is also well matched by slab-mode cutoff frequency. Based on their displacement profile, the d-mode group is likely a strong candidate to large overlap with optical modes.
On the other hand, the f-group modes are cantilever-like modes and exhibit and odd symmetry relative to the z=0 plane (see Fig. 1b), resulting in a negligible interaction with sd-cavity optical modes. In the dd-cavities, however, the symmetric combination of upper and lower-disk f-modes strongly modulate the air-gap between the disks, thus having potential to strongly interact with the dd-cavity optical modes. Indeed, these symmetric f-modes are similar to those explored in previous double-disk devices32,33, but due to their large azimuthal number they can readily vibrate in the 10GHz frequency range.
The spatial overlap between optical modes and mechanical modes is necessary but not sufficient for a large optomechanical coupling. The optomechanical interaction in these high refractive index structures occurs due to a combination of the photo-elastic effect1 and deformation of the cavity boundaries2. The calculated optomechanical coupling rate must consider both effects, g0=gpe+gmb, where gpe stands for the volumetric photo-elastic contribution (pe-effect), and gmb for the cavity moving-boundaries contribution (mb-effect). In Fig. 3a,b we show the calculated photo-elastic (gpe, red), the moving-boundary (gmb, green) and total coupling rate g0 (bars) for the phase-matched mechanical modes in the single (Fig. 3a) and double-disk (Fig. 3b) structures. We focus on the mechanical modes that are phase-matched with fundamental TM (transverse magnetic) optical mode (Fig. 3c) since the TM-modes exhibits the largest coupling rate and potentially higher optical quality factors34.
For both sd and dd-structures, the whispering, Rayleigh, and dilatational mode groups can be identified in Fig. 3a,b. The relative contributions from the pe and mb-effects, however, varies significantly for each structure and mode group. To understand this in detail, we analyze the weighting function role played by the optical field in the mb and pe-effects (see Methods). For the mb-effect, the normal mechanical displacement (u) is weighted by the optical field as (see Supplementary Information),
where Er, Eϕ and Ez are the energy-normalized electric field components, δεmb=ε1−ε2 and with and being the permittivities of silicon and air, respectively. The spatial dependence of the three terms in Eq. 2: gmb_terms are shown in Fig. 3d,e for both sd and dd-cavities. It is evident that the mb-effect is dominated by the azimuthal component (−δεmb|Eϕ|2) in both structures. The opposite sign of the azimuthal term relative to the radial contribution — due to the π phase difference between forward and backward azimuthal field components — substantially distinguishes the backward from the forward Brillouin optomechanical interaction. The weighting terms peaks around r=4.6µm also hint which mechanical modes should benefit from the mb-effect.
The photo-elastic (pe-effect) contribution is a bit more intricate due to the anisotropic nature of the strain-induced permittivity perturbation. Each component is calculated from the permittivity perturbation tensor, defined as , in which p is the photoelastic tensor of the isotropic silicon and S is the strain tensor induced by the mechanical waves (see Methods). For TM optical mode, however, the main permittivity perturbations are the vertical () and azimuthal () components since they are weighted by the dominant fields Ez and Eϕ. Despite the complexity of the pe-effect, the dominant photo-elastic coefficient in silicon is p11, hence an insight about which modes will lead to a strong pe-effect can be obtained analyzing only the azimuthal () and vertical () perturbations.
The w-mode group has the largest optomechanical coupling rate and give rise to several peaks in the high frequency range (20~28GHz), for both sd (Fig. 3a) and dd cavities (Fig. 3b), which are unique to BBS due to the large mechanical azimuthal numbers imposed by the phase-matching condition. Due to their dominant displacement components (ur, uϕ), the largest strain components is azimuthal, Sϕϕ≈Muϕ/r. Such a large azimuthal strain leads to a pe-effect dominated optomechanical coupling, reaching the highest coupling rate for the single-disk, g0/2π≈61kHz at 24.34GHz for the 16th radial order w-mode (w16); the strain profile of selected w-modes are shown in Fig. 3f. In this mode group, accounts for about 70% of the total coupling coefficient (see Supplementary Information, Table S2). The double-disk counterpart of this mode family (wdd) behaves similarly, yet with a reduced strength (g0/2π≈35kHz for the ) due to the larger optical and mechanical mode volumes. In either designs, the tiny in-plane displacement of this mode group always ensures a negligible contribution from the mb-effect.
Another interesting feature of this mode group is the peaked g0 distribution in Fig. 3(a,b), which can be understood by inspecting the overlap between the dominant azimuthal strain component Sϕϕ (Fig. 3f) and the optical mode profile (Fig. 3c). Despite the strain oscillations along the radial direction, there is a net tensile strain (Sϕϕ>0) region indicated by the dashed arrow in Fig. 3f; as the mechanical frequency increases, the net strain region shifts radially inwards and sweeps the spatial matching between the strain and optical fields. The fine mechanical frequency spacing between adjacent w-modes and such a sweeping overlap lead to the observed peaked g0 distribution in Fig. 3(a,b). The physical origin of this net positive strain region relies on the hybrid longitudinal-transverse nature of the w-group, which can be precisely traced using the analytical solution of an infinite cylinder: the fast radial oscillation periods seen in Fig. 3f arise from the transverse-wave contribution to this mode, whereas the more slowly varying net positive strain regions are caused by longitudinal-wave contribution (see Supplementary Information, Fig. S4).
The Rayleigh mode (r1), which has the lowest resonant frequency (at 11.12GHz, Fig. 3a), has a dominant radial displacement (ur) in the sd-cavity, which results in the edge-localized radial strain shown in Fig. 3g. On the other hand, its flexural odd counterpart in the dd-cavity, , has also a large vertical (uz) component. Combined with the large azimuthal component, this leads to a dominant shear Sϕz strain, which is shown in Fig. 3g. Due to its edge localization, away from the optical mode maximum, a minor role from the mb-effect is expected for the r1 mode in the sd-cavity. However, the boundary contribution for in the dd-cavity is outstanding, indeed this mode (11.15GHz) has the second largest coupling rate among all the dd-cavity modes, reaching g0/2π=31kHz. Furthermore, the dd-cavity has a unique feature as its mb-effect strength can be tailored by adjusting the slot height between the two disks; in Fig. 3i we show that the total coupling rate (g0) of the mode can be improved by 300% by reducing the slot height from 150nm to 50nm.
Finally, a very high optomechanical coupling rate could be expected for the d-mode group (16GHz <Ω/2π<18GHz). These modes are striking examples that intuition based on mode shapes may fail when predicting the BBS interaction strength in microstructures. The d-modes display not only a large vertical strain but also a large deformation of the cavity boundaries, as shown in Fig. 3h. Indeed, these two effects are very strong individually but their opposite sign lead to a cancellation effect, a clear competition between the mb (gmb) and pe-effects (gmb)22. For the double-disk structure, whose optical weighting function is shown in Fig. 3e, the slot effect enhancement does not readily improve the optomechanical coupling with the dilatational modes, this is due to a balanced contribution from the azimuthal (dashed black line) and vertical field (dashed blue line) components along the slot boundaries, which oppositely contribute to the mb-effect and almost cancel it.
The calculated BBS coupling rates for both sd and dd-cavities improve over a few orders of magnitude compared to coupling rates reported for larger silica microspheres (g0/2π≈20Hz)35 and crystalline resonators (g0/2π≈61Hz)20. Such high coupling rates ensure that the proposed designs will have a low input power threshold to achieve stimulated Brillouin lasing, despite the lower quality factors often achievable with silicon cavities. For instance, we assume an optical mode resonant around the wavelength of 1550nm, with an intrinsic optical quality factor of 2×105 (κp/2π=κs/2π≈1.2GHz), and mechanical quality factor of 103, which are conservative optical and mechanical mode parameters. If these cavities designs are employed in the proposed compound cavity scheme, and the mechanical frequency matches the optical frequency splitting, i.e., (Δs=Δp=0) in Eq. 1: threshold. The predicted threshold power for SBS lasing is only Pth≈8mW for the single-disk w16 mode; Pth≈31mW for the double-disk mode.
For the cantilever-like flexural mode of the dd-cavity (), the threshold power is Pth≈17mW (assuming the same mechanical quality factor). This threshold power for the -mode can be reduced even further for smaller double-disk separation. For instance, if tSiO=50nm, it is possible to achieve g0/2π≈75kHz (Fig. 3i), leading to a threshold power of only Pth≈3mW (considering the same optical and mechanical losses).
Experimentally, in order to ensure the simultaneous resonant condition, a set of coupled optical cavities with varying coupling gaps could be fabricated. The importance of the compound cavity scheme also becomes clear if we compare the results above with a single cavity (single or double-disk) input threshold power, which is predicted by eq:threshold assuming a resonant stokes signal Δs=0 and the optimal pump-detuning (Δp=Ωm). Using the same optical and mechanical linewidth above, and a single-photon optomechanical coupling rate (, see Methods), the threshold for Brillouin lasing would be as high as Pth≈3.2W for the w16, and Pth≈12.3W for the mode. Such a high input powers are rather impractical due to strong detrimental effects they induce, such as strong nonlinear light absorption or even thermal damage of the bus waveguides.
Our results provide a clear guideline towards the observation of backward Brillouin scattering lasing in integrated CMOS-compatible silicon cavities. Both the single-disk and double-disk designs were shown to enable large optomechanical coupling rates, a few orders of magnitude larger than previously reported resonators, allowing for mW-level input power lasing thresholds. The single-disk device, despite its simplicity, was shown to enable mechanical modes at very high frequencies (25GHz) with large optomechanical coupling rates. The double-disk device, although exhibiting a lower optomechanical coupling for high frequency modes, was shown to drastically enhance the interaction with flexural modes around 11GHz. Although we concentrated our discussion on silicon-based devices, our results could be adapted to similar structures fabricated from other high index materials, such as III-V, Si3N4 and SiO2.
The frequency splitting simulation was performed using a two-dimensional approximation to the actual sd-cavity. In this approximation, the modes of the coupled infinite cylinders are calculated while constraining the out-of-plane wave number,
where m/r and zm,1/r are the azimuthal and radial components of the wave vector with norm k0n, k0 is the free-space wave number for λ=1.55µm, n=3.5 is the silicon refractive index, r=5µm, m=35 is the optical azimuthal number for the TM-mode of the sd-cavity and zm,1 is the first zero of the Bessel function Jm(z). This is equivalent to the Marcatilli effective index method and we verified that it leads to an electric field envelope that agrees well with the numerical mode obtained from the axisymmetric calculation.
We obtain the dispersion relation Ω(M) by solving the eigenfrequency problem derived from the full-vectorial elastic wave equation through the finite-element method (see Supplementary Information), due to the highly multimode character of the mechanical dispersion, we show in Fig. 2a–c a grey color shading proportional to the mechanical density of modes instead of the calculated pairs (ΩM, M). The mechanical density of modes is calculated as , where f and g are Gaussian weighting functions with a normalized product, full-width-half-maximum (FWHM) wavenumber σM=0.1, FWHM-frequency σΩ/2π=117.5MHz and given one azimuthal acoustic number are calculated each of the frequencies from the elastic wave equation.
The resulting energy exchange between the optical pump and Stokes wave can be modeled using coupled mode theory27,28 (see Supplementary Information), which leads to a set of coupled equations for the amplitudes of the pump wave (ap), stokes wave (as) and mechanical wave (b),
where ai’s are normalized such that |ai|2 is the intra-cavity photon number and b is normalized such that |b|2 is the phonon number. is the frequency detuning between the pump (i=p) or stokes wave (i=s) relative to the optical cavity mode frequencies . κp,s represents the optical loss rate for each mode, (Ωm,Γm) are the mechanical mode frequency and damping rate, respectively; the optomechanical coupling rate is and represents the photon coupling rate between the stokes and the pump wave induced by the zero-point fluctuation of the mechanical mode, which will be calculated shortly using electromagnetic perturbation theory2,36. Note that is related to the more usual single-cavity coupling rate as . The factor 1/2 arises because the coupled optical mode is distributed within two optical cavities whereas the mechanical mode is localized within a single cavity due to the air gap between the cavities. κe is extrinsic coupling rate to the feeding waveguide carrying a photon-flux |sp|2. Fth is a white-noise random thermal force responsible for driving the mechanical motion.
where the permittivity differences are given by δεmb=ε1−ε2 and , in which and are the permittivities of the silicon and air, respectively. is the surface-normal component of the displacement vector u (normalized to unit); is zero-point fluctuation of the mechanical mode with effective mass meff; the fields and Dj, are boundary-tangential electric field and boundary-normal electric displacement field to the cavity surface S of the pump (j=p) or scattered (j=s) optical mode (energy-normalized). The pe-effect contribution is given by ref. 36 (see Supplementary Information),
where is the photo-elastic perturbation in the permittivity inside the cavity volume V, p is the photoelastic tensor of silicon, and S=xzpfsu is the strain tensor induced by the mechanical waves. The optical and elastic mode profiles are numerically calculated using the finite-element method (see Supplementary Information).
Si refractive index nSi=3.5, silica refractive index , air refractive index nair=1.0, wavelength of interest λ=1550nm, Si Young’s modulus ESi=170GPa, silica Young’s modulus GPa, Si Poisson’s ratio νSi=0.28, silica Poisson’s ratio , Si density mass ρSi=2329kg/m3, silica density mass kg/m3 and Si photo-elastic coefficients37 p11=−0.09, p12=0.017 and p44=−0.0535. Here we neglect silicon anisotropy and assume the values of ESi and νSi along to principal crystal axes38.
How to cite this article: Espinel, Y. A. V. et al. Brillouin Optomechanics in Coupled Silicon Microcavities. Sci. Rep. 7, 43423; doi: 10.1038/srep43423 (2017).
Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors would like to acknowledge Omar Florez and Paulo Dainese for fruitful discussions. This research was funded by the Sao Paulo State Research Foundation (FAPESP) (Grants 2012/17765-7 and 2012/17610-3), the National Counsel of Technological and Scientific Development (CNPQ), and the Coordination for the Improvement of Higher Education Personnel (CAPES).
The authors declare no competing financial interests.
Author Contributions Y.E. performed the numerical simulations and conceived the idea with help from G.W. and T.A. G.O.L and F.S helped Y.E. with the numerical simulation and analytical analysis. All authors discussed the results and their implications and contributed to writing this manuscript.