# Spatial localization and pattern formation in discrete optomechanical cavities and arrays

J. Ruiz-Rivas, Giuseppe Patera, C. Navarrete-Benlloch, E. Roldán, G.J. de Valcárcel

► To cite this version:

J. Ruiz-Rivas, Giuseppe Patera, C. Navarrete-Benlloch, E. Roldán, G.J. de Valcárcel. Spatial localization and pattern formation in discrete optomechanical cavities and arrays. New Journal of Physics, 2020, 22 (9), pp.093076. 10.1088/1367-2630/abb73e. hal-02966280

HAL Id: hal-02966280

https://hal.science/hal-02966280v1

Submitted on 16 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Spatial localization and pattern formation in discrete optomechanical cavities and arrays

To cite this article: J Ruiz-Rivas et al 2020 *New J. Phys.* 22 093076

View the [article online](#) for updates and enhancements.

## You may also like

- Response of some Oat Cultivars to Spraying with NPK Nanofertilizer

Rasha Adil Al-Bakri and Ahmed Mohammed Lehmood

- Instabilities and turbulence in stellarators from the perspective of global codes

E. Sánchez, A. Bañón Navarro, F. Wilms et al.

- Ordinary least square and maximum likelihood estimation of VAR(1) model's parameters and its application on covid-19 in China 2020

Y Nalita, R Rahani, E R Tirayo et al.
<figure><img src="image_2.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

27 April 2020

REVISED

25 August 2020

ACCEPTED FOR PUBLICATION

10 September 2020

PUBLISHED

23 September 2020

Original content from
this work may be used
under the terms of the
Creative Commons
Attribution 4.0 licence.

Any further distribution
of this work must
maintain attribution to
the author(s) and the
title of the work, journal
citation and DOI.

<figure><img src="image_3.png" /></figure>

PAPER

# Spatial localization and pattern formation in discrete optomechanical cavities and arrays

J Ruiz-Rivas¹, G Patera², C Navarrete-Benlloch³,⁴,⁵,*, E Roldán¹,* and G J de Valcárcel¹,*

¹ Departament d'Óptica, Universitat de València, Dr. Moliner 50, 46100-Burjassot, Spain

² Univ. Lille, CNRS, UMR 8523-PhLAM-Physique des Lasers Atomes et Molécules, F-59000 Lille, France

³ Max-Planck-Institut für die Physik des Lichts, Staudtstrasse 2, 91058 Erlangen, Germany

⁴ Wilczek Quantum Center, School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai 200240, People's Republic of China

⁵ Shanghai Research Center for Quantum Sciences, Shanghai 201315, People's Republic of China

* Authors to whom any correspondence should be addressed.

E-mail: carlos.navarrete-benlloch@mpl.mpg.de, eugenio.roldan@uv.es and german.valcarcel@uv.es

**Keywords:** cavity optomechanics, optomechanical arrays, photon and phonon localization, spatial solitons

## Abstract

We investigate theoretically the generation of nonlinear dissipative structures in optomechanical (OM) systems containing discrete arrays of mechanical resonators. We consider both hybrid models in which the optical system is a continuous multimode field, as it would happen in an OM cavity containing an array of micro-mirrors, and also fully discrete models in which each mechanical resonator interacts with a single optical mode, making contact with Ludwig and Marquardt (2013 *Phys. Rev. Lett.* **101**, 073603). Also, we study the connections between both types of models and continuous OM models. While all three types of models merge naturally in the limit of a large number of densely distributed mechanical resonators, we show that the spatial localization and the pattern formation found in continuous OM models can still be observed for a small number of mechanical elements, even in the presence of finite-size effects, which we discuss. This opens new venues for experimental approaches to the subject.

## 1. Introduction

The emergence of patterns that spontaneously break some spatial symmetry is widespread in nonlinear optical systems, especially in large aspect-ratio cavities. Such patterns show up across the plane transverse to the light propagation direction—hence the name ‘transverse patterns’ to refer to them—and have been theoretically and experimentally investigated in many different nonlinear optical cavities [1, 2]. Their study constitutes a well developed discipline in modern nonlinear optics. In fact, these nonlinear patterns belong to the wider class of dissipative structures, which are structures that self-sustain out of thermal equilibrium by continuous interchange of energy with the environment [3]. Part of the interest of this research program in nonlinear optics, apart from its intrinsic physical relevance, lies in the potential for optical information storage and processing of a particular type of pattern, namely cavity solitons, which are localized structures that can be individually written, erased, and even moved without affecting neighboring structures [4–7].

Naturally, most of the studies on optical transverse patterns so far have considered cavities containing usual passive or active nonlinear materials, like two-level systems, Kerr media, or second-order nonlinear crystals, to cite a few. More recently, optomechanical (OM) systems have begun to be considered in this respect. In [8] we analyzed theoretically the possibility of using the nonlinear coupling between the cavity field and a deformable mechanical element to generate transverse patterns in OM cavities. These are conceptually simple systems, consisting of an optical resonator with mechanical degrees of freedom that couple to the light oscillating inside it [9]. The coupling appears either through radiation pressure (e.g., when the mechanical degree of freedom corresponds to the oscillation of a perfectly reflecting cavity mirror) or through dispersive effects (e.g., when the mechanical degrees of freedom correspond to the local
displacement of a partially transmitting membrane). Also, the formation of transverse patterns in OM cavities containing cold atoms has been considered [10–12]. These systems are receiving intense and continued attention mainly in the context of modern quantum technologies, where phenomena such as cooling [13–15], induced transparency [16, 17], squeezing [18–20], as well as quantum-coherent effects [21] have been demonstrated over the last decade. Generating dissipative structures in OM cavities could endow them with new capabilities. It could also lead to a pattern forming system in which quantum fluctuations play an important role: as OM cavities have demonstrated their ability to work within the quantum regime, they open the possibility of studying quantum dissipative structures under the strong influence of quantum fluctuations. Spontaneous dissipative structures indeed present quantum correlations among different spatial points within the structure [22]. These effects are stronger for parameters close to bifurcations, but some quantum properties remain even far from such points. For example, under certain conditions, the spontaneous spatial symmetry breaking induced by the structure gives rise to quantum correlations (squeezing, in particular) that are non-critical, in the sense that they do not depend on the distance to the pattern-forming instability [23–25]. Most of these exciting quantum phenomena have been predicted for optical parametric oscillators [23–29], but unfortunately the models are still far from realistic experimental implementations [30].

Contrarily, OM cavities have demonstrated the feasibility of the simplest models as well as a large versatility because of the variety of possible platforms, materials and designs [9, 31, 32]. As a first step toward understanding this fully quantum picture, it is important to characterize the conditions required for pattern formation in OM cavities.

Within the context of extended OM systems, Rakich and Marquardt [33] have recently formulated a quantum theory of continuum optomechanics that is closely related to our theory [8]. The theory is intended for treating OM interactions occurring along extended waveguides, which are naturally space-dependent problems in which Brillouin scattering is the essential coupling, thus connecting Brillouin physics with optomechanics. Below we address the formal connection between continuum optomechanics and pattern forming OM cavities.

In [8] we demonstrated the feasibility of pattern formation in OM cavities by modeling the mechanical element, be it an end mirror or an intracavity membrane, as a continuously deformable element. Interestingly, a certain necessary condition, alien to other nonlinear optical cavities, must be fulfilled for nonlinear pattern formation in OM cavities, namely that the mechanical array must possess a sufficiently homogeneous mode, i.e. a state that is invariant under translations in the plane transverse to the light propagation direction (cavity axis). This means that the equilibrium displacements of the mechanical elements under a spatially uniform illumination must equal each other, what is called a homogeneous spatial state. In [8] we suggested a way for implementing such condition through a quasi-one-dimensional membrane, that is, a membrane clamped by a large aspect-ratio frame, such that only one or a few transverse modes could be excited along the short direction. Under these conditions, the system naturally develops a one-dimensional homogeneous mode along the long direction, which becomes unstable through pattern forming instabilities under appropriate parameter settings. One of the limitations of that scheme is that only one-dimensional (1D) transverse-patterns can be generated.

In the present paper we investigate a conceptually different possibility of pattern formation to that of [8] by considering an OM cavity with an oscillating micro-structured end mirror consisting of an array of $N$ weakly-coupled micro-mirrors, see figure 1 for a sketch of the 1D case; of course, the system can also be implemented in 2D. This configuration, indeed, allows automatically a homogeneous mode. Below we show that in the limit of large $N$, one recovers the continuous linear-coupling model in [8], and hence all the predictions in that work apply when the number of elements in the array is large. The interesting point is that one can consider also the limit of small $N$ and study the transition from the continuous to the discrete limit. Among the interesting numerical results that we show below, it is remarkable the fact that a discrete analog of the continuous-limit cavity solitons can be observed with a relatively small number of coupled micro-mirrors, say $N \approx 10$.

Of course, the problem we are addressing is closely connected with the theory of OM arrays [34–37] as our model can be thought of as a special implementation of these systems in which all the OM elements are nonlinearly driven through the same multimode intracavity optical field. As the continuum optomechanics theory is recovered from the continuum limit of the OM array model, this opens the way for a formal connection between pattern formation in OM array cavities and continuum optomechanics.

After this introduction, we present our model for the OM cavity with microstructured end-mirror in section 2. In section 3 we establish the connection with our continuous model of [8] and in section 4 we present some numerical results showing, in particular, the transition from the quasi-continuum limit to a small number of mechanical elements. Then, in section 5, we establish a formal connection between the
<figure><img src="image_4.png" /><figcaption>Figure 1. Sketch of the OM cavity. A coherent field of slowly-varying complex amplitude $A_{\text{inj}}(z, \mathbf{r}, t)$ is injected into the cavity through a coupling mirror placed at the left end ($z=0$). The intracavity field is given by the superposition of two counter-propagating waves with slowly varying complex amplitudes $A_{\pm}(z, \mathbf{r}, t)$. The right end of the cavity ($z=L$) is formed by an array of N micro-mirrors that are weakly coupled by means of springs of stiffness $\kappa_{\perp}$. The array is one-dimensional in the figure for simplicity, but we consider two-dimensional ones along the text. Each mirror can oscillate back and forth along the longitudinal direction ($z$ axis) with angular frequency $\Omega_m$. We measure the displacement $q_j(t)$ of each micro-mirror with respect to the reference flat state at $z=L$ (gray dashed rectangle) when the cavity is not illuminated. For two adjacent micro-mirrors, $a$ is the distance between their centers and $b$ is the distance between their edges, when they are at rest.</figcaption></figure>

OM array models [36] and the continuum optomechanics, as well as we show some relevant numerical results. Finally, in section 6 we give our main conclusions.

## 2. Model

Consider an optical cavity with large-area mirrors, one of which is plane, partially transmitting, and immune to radiation pressure because of its stiffness and mass, while the other is a perfectly reflecting array of micro-mirrors. We assume that they can be described by N identical spring-mass systems, weakly coupled to each other by means of springs of stiffness $\kappa_{\perp}$ (see figure 1). Also, we assume that the tilt induced by the coupling springs is negligible, so that each micro-mirror oscillates along the z-axis with angular frequency $\Omega_m$.

The field injected in the cavity through the coupling mirror is assumed to be a paraxial, coherent beam

$$E_{\text{inj}}(z, \mathbf{r}, t) = i\mathcal{V}A_{\text{inj}}(z, \mathbf{r}, t) e^{i(k_L z - \omega_L t)} + \text{c.c.}, \quad (1)$$

where $\mathbf{r} = (x, y)$ denotes the position in the plane transverse to the cavity axis (z-axis), and $\mathcal{V}$ is a constant having the dimensions of voltage, which we choose as $\mathcal{V} = \sqrt{\hbar\omega_c/4\epsilon_0 L}$ in order to make contact with quantum optics (see appendix A), $L$ being the length of the unperturbed cavity and $\omega_c$ being the frequency of the longitudinal cavity mode closest to the injected frequency $\omega_L$, with corresponding wave vector $\mathbf{k}_L = \omega_L/c$.

The generic intracavity field $E(z, \mathbf{r}, t)$ can be written as

$$E(z, \mathbf{r}, t) = i\mathcal{V} (A_+ e^{ik_L z} + A_- e^{-ik_L z}) e^{-i\omega_L t} + \text{c.c.}, \quad (2)$$

which is the superposition of two waves with slowly varying complex amplitudes $A_{\pm}(z, \mathbf{r}, t)$, propagating along the positive ($A_+$) and negative ($A_-$) z direction. Note that we are considering a single longitudinal cavity mode (that closest to the injection's frequency), which we justify at the end of the section. With similar assumptions as those in [8, 38] (see appendix B for the derivation), the field $A_+(z=L, \mathbf{r}, t)$ at the surface of the mechanical array, which we denote by $A(\mathbf{r}, t)$, obeys the following evolution equation

$$\partial_t A = \gamma_c \left( -1 + i\Delta + i l_c^2 \nabla_{\perp}^2 + i \frac{G}{\gamma_c} Q \right) A + \gamma_c \mathcal{E}. \quad (3)$$

Here $G = \omega_L/L$ is the frequency pull parameter [9], $\gamma_c = cT/4L$ is the cavity damping rate (with $T$ the transmissivity of the fixed mirror), $\Delta = (\omega_L - \omega_c)/\gamma_c$ is the dimensionless detuning parameter, $\nabla_{\perp}^2 = \partial_x^2 + \partial_y^2$ is the transverse Laplacian, $\mathcal{E}(\mathbf{r}, t) = 2T^{-1/2}A_{\text{inj}}(L, \mathbf{r}, t + t_c)$ is a scaled version of the injected field amplitude ($t_c = 2L/c$ is the cavity round-trip time) and $l_c^2 = 2L/k_L T$ is the square of the diffraction length. The latter is the characteristic size for which diffraction is more effective, similarly to how a Gaussian beam diffracts more at its waist, so that details of length $l_c$ of a given image will suffer more diffraction than the rest of details. Furthermore, being any image a multimode coherent superposition of spatial modes, the number of disks having a radius $l_c$ contained in the image provides a rough estimate for the number of transverse modes required to describe it. In (3) we have introduced a field $Q(\mathbf{r}, t)$ that measures the local displacement of the mechanical array with respect to the reference state, say the 'flat
state', for which $Q = 0$ for any $\mathbf{r}$ (corresponding to the gray dashed rectangle in figure 1). Next we derive its equation of motion.

We describe the displacement of the mechanical array in terms of the individual displacements $\{q_j\}_{j \in \mathbb{Z}^2}$ of its constituent micro-mirrors, labeled by a double index $\mathbf{j} = (j_1, j_2)$ in a 2D configuration, as

$$Q(\mathbf{r}, t) = \sum_{j} q_j(t) w_j(\mathbf{r}), \qquad (4)$$

where $w_j(\mathbf{r})$ is the characteristic function of the j-th mirror which equals 1 when $\mathbf{r}$ is on the surface of micro-mirror $\mathbf{j}$ and is zero otherwise. In the following we assume for simplicity that the micro-mirrors are much larger than the separation between them, that is, $a \gg b$ in figure 1. Note that the latter condition implicitly assumes that the micro-mirrors are square-shaped, although their specific shape is unsubstantial, as long as all our previous assumptions hold (e.g., the space left between them is smaller than their size). Each of these displacements is assumed to satisfy the equation of motion of a damped and forced harmonic oscillator, the force acting on mirror $\mathbf{j}$ having two contributions, $F_j = F_j^{(\text{RP})} + F_j^{(\perp)}$, respectively coming from radiation pressure and from the coupling to neighboring mirrors. The first contribution is readily obtained by integrating the radiation pressure (A.5) over the surface $S_j$ of the corresponding micro-mirror

$$F_j^{(\text{RP})} = \hbar G \int_{S_j} d^2 \mathbf{r} |A(\mathbf{r}, t)|^2, \qquad (5)$$

where $k_c = \omega_c/c$. As for the force coming from the coupling to neighboring mirrors, we assume that it originates from a potential that harmonically couples neighbors as $V_j^{\perp} = \kappa_{\perp} \sum_{\langle l \rangle_j} (q_l - q_j)^2/2$, where $\langle l \rangle_j$ means that the sum is performed over nearest neighbors, hence the corresponding force is computed as $F_j^{(\perp)} = -\partial V_j^{\perp}/\partial q_j$, which reads

$$F_j^{\perp} = \kappa_{\perp} \mathcal{L}_{\nu}[q_j]. \qquad (6)$$

where

$$\mathcal{L}_{\nu}[q_j] \equiv \sum_{\langle l \rangle_j} (q_l - q_j), \qquad (7)$$

being $\nu$ the coordination number. For instance, for an inner point in a two-dimensional square lattice ($\nu = 4$),

$$\mathcal{L}_{\nu}[q_j] = q_{j+u_x} + q_{j+u_y} + q_{j-u_x} + q_{j-u_y} - \nu q_j, \qquad (8)$$

where $u_x = (1, 0)$ and $u_y = (0, 1)$. For points at the boundary trivial modifications can be done to this expression. Note that $\mathcal{L}_{\nu}[q_j]$ is the finite-difference approximation to the Laplacian of a continuous field $Q(\mathbf{r}, t)$ of which $\{q_j\}$ is a sampling.

Putting everything together we get

$$\ddot{q}_j + \gamma_m \dot{q}_j + \Omega_m^2 q_j - \frac{\kappa_{\perp}}{m} \mathcal{L}_{\nu}[q_j] = \frac{\hbar}{m} G \int_{S_j} d^2 \mathbf{r} |A(\mathbf{r}, t)|^2, \qquad (9)$$

with $\gamma_m$, $\Omega_m$ and $m$ the damping rate, oscillation frequency, and mass of the micro-mirrors, respectively. Equation (9) together with the optical field (3) and definition (4) form the equations of our model.

For the sake of completeness, we give in the remaining of this section a rough estimate for some relevant parameters of the model. We stress that we found pattern formation and spatial localization in a very broad region of the parameter space (see below), which means that the following estimates are illustrative, and not a narrow condition that needs to be matched experimentally. Indeed, parameters may vary significantly depending to the technology chosen for the experimental implementation (see for example section 5). The most relevant parameter is the diffraction length $l_c$ in relation to the characteristic size $a$ of a single micro-mirror. $l_c$ increases with the cavity length $L$, and decreases with the frequency of the injected field or the transmissivity $T$ of the coupling mirror. Hence, for example, for $L \sim 1 \text{ mm} - 1 \text{ cm}$, and green-light injection ($\lambda \sim 0.5 \text{ µm}$), one gets $l_c \sim 50 \text{ µm} - 1 \text{ mm}$, evidencing the tunability of this parameter. Considering $a \approx l_c$, then, an array of about $10 \times 10$ micro-mirrors would have a total size of about $0.5 \text{ mm} \times 0.5 \text{ mm}$ to $1 \text{ cm} \times 1 \text{ cm}$.

We recall that our model considers just one cavity longitudinal mode, which is obviously justified as far as no other longitudinal mode becomes excited. This is granted whenever the frequency separation between neighboring longitudinal modes, i.e. the so-called cavity free spectral range (FSR), is large enough as compared to the relevant mechanical frequencies, which is the case we consider. Otherwise a three-mode parametric instability [39–42] that involves a mechanical mode and two optical modes with different longitudinal order can settle in the system. That kind of instability can limit the performance of
gravitational wave interferometers, whose very long arms (km sized) lead to FSR compatible with mechanical frequencies existing in the system [39, 40]. Our relatively small cavities (5 orders of magnitude shorter) have FSR exceeding the GHz, therefore well above typical mechanical frequencies. Our model does contain however three-mode parametric instabilities between transverse modes with the same longitudinal order (which are assumed quasi-degenerate in frequency, as required for a pattern-forming cavity with translational invariance in the transverse plane). Nevertheless we are considering here a dynamical regime where the homogeneous steady-state solution is unstable with respect to static pattern-forming instabilities only, which roughly requires working with red-detuned injection, as shown by a standard linear stability analysis such as the one we performed in our previous work of reference [8]. Note that this also excludes the usual OM self-pulsing or Hopf instability [9], which occurs when the optical injection is blue detuned with respect to the cavity mode frequency.

### 3. Continuous limit

In order to make contact with our previous work [8] and with [33], we now consider the limit in which the number of micro-mirrors per unit length tends to infinity while keeping a finite mass density and sound speed. The following reasoning is valid for a generic 2D lattice with lattice constant $a$ and generated by a pair of linearly-independent vectors $\{\mathbf{u}_1, \mathbf{u}_2\}$. Each mirror of the array is centered at the nodes of the lattice, which are univocally identified by the indices $j = (j_1, j_2) \in \mathbb{Z}^2$ (see figure 2 for the example of a triangular lattice). We first write the displacements as a function of the mechanical field $Q$ as

$$q_j = \int_{\mathbb{R}^2} \frac{d^2\mathbf{r}}{a^2} Q(\mathbf{r}) w_j(\mathbf{r}). \qquad (10)$$

Next, using the immediate properties

$$\int_{S_j} d^2\mathbf{r} |A(\mathbf{r}, t)|^2 = \int_{\mathbb{R}^2} d^2\mathbf{r} |A(\mathbf{r}, t)|^2 w_j(\mathbf{r}), \qquad (11a)$$

$$\int_{\mathbb{R}^2} d^2\mathbf{r} Q(\mathbf{r}) w_{j+m}(\mathbf{r}) = \int_{\mathbb{R}^2} d^2\mathbf{r} Q(\mathbf{r} - a[m_1\mathbf{u}_1 + m_2\mathbf{u}_2]) w_j(\mathbf{r}), \qquad (11b)$$

with $\mathbf{m} = (m_1, m_2) \in \mathbb{Z}^2$, the equation (9) for the motion of a generic displacement $q_j$ is turned into

$$\partial_t^2 Q(\mathbf{r}) + \gamma_m \partial_t Q(\mathbf{r}) + \Omega_m^2 Q(\mathbf{r}) - \Omega_\perp^2 \sum_{n=1}^{\nu} (Q(\mathbf{r} - a\mathbf{v}_n) - Q(\mathbf{r})) = \frac{\hbar a^2}{m} G |A(\mathbf{r})|^2 \qquad (12)$$

with $\Omega_\perp = \sqrt{\kappa_\perp/m}$, and where $\{\mathbf{v}_n\}_{n=1,2,...,\nu}$ are the vectors that connect a given lattice node with its $\nu$ nearest neighbors. Note that we are considering a micro-mirror that is not at the boundary of the flexible mirror (as in the continuous limit the fields extend up to infinity), and hence have used (6). As specific examples, let us consider the square and triangular lattices. The first one is generated by $\mathbf{u}_1 = \mathbf{u}_x$ and $\mathbf{u}_2 = \mathbf{u}_y$, with nearest neighbors connected through the vectors $\mathbf{v}_1 = \mathbf{u}_x = -\mathbf{v}_2$ and $\mathbf{v}_3 = \mathbf{u}_y = -\mathbf{v}_4$. In the case of the triangular lattice, which is represented in figure 2, we have $\mathbf{u}_1 = \mathbf{u}_x$ and $\mathbf{u}_2 = (1, \sqrt{3})/2$, while nodes are connected with their neighbors through $\{\mathbf{v}_n = (\cos \frac{\pi n}{3}, \sin \frac{\pi n}{3})\}_{n=1,...,6}$.

For any lattice, we next consider the continuous limit by assuming that the lattice constant $a$ is much smaller the scale in which spatial variations in the field $Q$ take place (long-wavelength approximation). In this limit we can approximate

$$Q(\mathbf{r} - a\mathbf{v}) \approx Q(\mathbf{r}) - a\partial_{\nu}Q(\mathbf{r}) + a^2\partial_{\nu}^2Q(\mathbf{r})/2, \qquad (13)$$

where $\partial_\nu = \mathbf{v} \cdot \nabla_\perp$, with $\nabla_\perp = (\partial_x, \partial_y)$. Using this expansion on equation (12) for a square lattice, while keeping finite both the speed at which transverse perturbations propagate in the flexible mirror $v = a\Omega_\perp$ and its surface mass density $\sigma = m/a^2$, we obtain

$$\partial_t^2 Q + \gamma_m \partial_t Q + (\Omega_m^2 - v^2 \nabla_\perp^2) Q = \frac{\hbar}{\sigma} G |A|^2. \qquad (14)$$

For the triangular lattice, one obtains the same form for the equation, but with a modified speed $v$, specifically enhanced by a factor $\sqrt{3}/2$. Remarkably, this is the same equation of motion we derived in the model of [8] (for the linear coupling case). For sake of later use, we note that in that work it was evidenced
<figure><img src="image_5.png" /><figcaption>Figure 2. Triangular lattice as an example of a possible geometry for our array of micro-mirrors (shown in gray). The two-dimensional lattice is generated by two linearly-independent unit vectors $\mathbf{u}_1$ and $\mathbf{u}_2$, and each micro-mirror is univocally determined by an index $j \in \mathbb{Z}^2$.</figcaption></figure>

the importance of the rigidity parameter $\rho = v/\Omega_m l_c$, which can be expressed in terms of our model parameters as

$$ \rho = \frac{a\sqrt{\kappa_{\perp}/m}}{\Omega_m l_c}. \qquad (15) $$

In the continuous limit, pattern formation requires $\rho$ smaller than a certain critical value that is on the order of 1 or above, depending on the light's detuning [8]. Considering then $a \approx l_c \approx 126 \ \mu\text{m}$ as in the previous section, common mechanical oscillation frequencies $\Omega_m$ on the order of MHz, and a sound speed on the broad range of $v \sim 10^2-10^3 \ \text{m s}^{-1}$, we obtain rigidity parameters $\rho \sim 0.8-8$, as required.

## 4. Numerical simulations

From the previous derivation of the continuous model we can conclude that, at least in the limit of a very large number of micro-mirrors, the system we are proposing is equivalent to that studied in [8]. This implies that all the results that we obtained in [8] for the linear coupling model apply in this limit, both the analytical (concerning homogeneous steady states and their stability properties) and the numerical ones (types of patterns, generalized bistability, temporal dynamics, etc). Of course, one must wonder how large must the density of micro-mirrors be for the results of the continuous model to still apply, as well as how do the results change when departing from such continuous limit. We have performed extensive numerical simulations of both the discrete (equations (3) and (9)) and continuous (equations (3) and (14)) models, and figures 3 and 4 summarize our main findings. For simplicity, we have restricted the simulations to one dimension, but similar conclusions are drawn in 2D. We have numerically simulated the continuous model by using the usual split-step method, which at any time step provides an approximation of the fields at certain space points. The same method can be applied to the discrete model, and in particular, we take M spatial points for the optical field at every micro-mirror, denoting by $(j, l)$ point $l$ of mirror $j$, so that the field amplitude $A(x)$ is represented by the array $\{A_{j,l}\}_{j=1,2,...,N}^{l=1,2,...,M}$, giving a total of $N \times M$ points. The next step consists in choosing a finite-differences form of the integral appearing in the mechanical equations (9). We have found that, for stability purposes, an integration rule of the type

$$ \int_{S_j} dx |A(x, t)|^2 \approx a \sum_{l=0}^{M+1} d_l |A_{j,l}(t)|^2, \qquad (16) $$

where $A_{j,0} = A_{j-1,M}$ and $A_{j,M+1} = A_{j+1,1}$, is what works best, that is, we use a discrete representation of the integral over mirror $j$ that includes the last point of the previous mirror and the first point of the next one. The weights satisfy the constrain $\sum_{j=0}^{M+1} d_j = 1$, and we have chosen a second order integration rule $\{d_l\}_{l=0,1,...,M+1} = \{1, 23, 24, 24, ..., 24, 23, 1\}/24M$ which seems to provide very good convergence properties.

Following our previous work [8] and for the sake of convenience, we express equations (3) and (9) in terms of dimensionless variables. To this purpose we define the following dimensionless versions of the
<figure><img src="image_6.png" /><figcaption>Figure 3. Normalized field amplitude (squared) |F(x)|² and mechanical field Z(x) in the steady-state as a function of the position, for a 1D system of finite size x ∈ [-40l_c, 40l_c]. The solutions have been obtained by numerical resolution of the discrete (solid blue) and continuous (dashed red) models described in the text, and under a top-hat illumination. In all cases M = 11, γ = 0.1, Ω = 10, Δ = -2.2, and ρ = 1.13. The injection values E_0^2 have been chosen in the region where solitons ((a)–(d)) or periodic patterns ((e)–(h)) are expected from the continuous model. For the discrete model we consider N micro-mirrors as specified in the figure, so that their individual size is a = 80l_c/N.</figcaption></figure>

mechanical displacement and optical field,

$$z_j = \frac{G}{\gamma_c} q_j, \qquad F = \frac{G}{\Omega_m} \sqrt{\frac{\hbar a}{m \gamma_c}} A. \qquad (17)$$

We also define the dimensionless injection E = (G/Ω_m)(ħa/mγ_c)¹/²E and write κ⊥/m = v²/a² in terms of the effective rigidity parameter ρ = v/Ω_m l_c which together with the detuning was shown to control the appearance of dissipative structures in the continuous model [8]. We finally introduce normalized versions of other parameters, namely γ = γ_m/γ_c, Ω = Ω_m/γ_c.

Combining this normalization with the discrete form (16) of the integral, and introducing into equations (3) and (9) dimensionless versions of time and space, τ = γ_c t and ¯x = x/l_c, respectively, we obtain the normalized equations

$$\frac{d^2 z_j}{d\tau^2} + \gamma \frac{dz_j}{d\tau} + \Omega^2 z_j = \rho^2 \Omega^2 \frac{l_c^2}{a^2} \mathcal{L}_2[z_j] + \Omega^2 \sum_{l=0}^{M+1} d_l |F_{j,l}|^2, \qquad (18a)$$

$$\partial_{\tau} F = (-1 + i\Delta + i\partial_{\bar{x}}^2 + iZ) F + E, \qquad (18b)$$

where we define the normalized mechanical field Z(x, t) = ∑_j z_j(t)w_j(x) and remind that L_2[z_j] refers to the finite-difference version of the Laplacian, see (7). Note that after rescaling the mechanical and optical amplitudes, the equations have become independent of the OM coupling represented by G. This means that, at the classical level, G plays no role in the dynamics of the system, other than as a scaling factor of the amplitudes and the injection. In contrast, G is expected to have a bigger impact at the quantum level, where the presence of quantum noise does not allow for the same kind of trivial elimination of such parameter.

In actual experiments a plane-wave injection is realized by a Gaussian beam having a waist much larger than the (finite) transverse dimensions of the OM cavity. However, for the numerical implementation, the use of the split-step algorithm requires either periodic boundaries or an infinitely extended system. For this reason, in order to simulate the finite size of the system, we use the well-known approach that considers a flat injection presenting a finite transverse size, while we keep periodic boundary conditions in the transverse plane. In particular, we take a top-hat injection profile with finite width, modeled as a super-Gaussian E(¯x) = E₀ exp(-x²/2σ_x²). In figure 3, we show in solid blue the stationary structures we
<figure><img src="image_7.png" /><figcaption>Figure 4. Normalized field amplitude (squared) |F(x)|² and mechanical field Z(x) in the steady-state as a function of the position, for a 1D system of finite size x ∈ [-40l_c, 40l_c]. Figures 4(a) and (b) show a localized structure that has been written around position x̄ = 12; an additional localized structure has been written around x̄ = -12 in figures 4(c) and (d), which clearly does not disturb the previous structure. The basic shape of the injection (E²) is represented in dashed red, on top of which a Gaussian profile with the proper width and position is initially fed in order to write the localized structures. The parameters in this simulation are N = 7, M = 11, γ = 0.1, Ω = 10, Δ = -2.2, ρ = 1.13, E₀ = √1.5, and σₓ = 23.</figcaption></figure>

have found in a spatial window x ∈ [-40l_c, +40l_c] for different number of micro-mirrors N (whose individual size is then a = 80l_c/N), taking M = 11 field points per micro-mirror. We are showing results for the special case γ = 0.1, Ω = 10, Δ = -2.2, and ρ = 1.13, and studied the spatial structures for two values of the injection, E₀² = 2.25 and 2.7 (with σₓ = 40). For these parameters the continuous limit predicts the appearance of, respectively, cavity solitons and periodic patterns [8]. In the figure we show in dashed red the corresponding structures found in the continuous limit for these same parameters. We see that cavity solitons are better captured with a small number of micro-mirrors than periodic patterns. This is so because we need enough space to hold such an extended structure, but they can still be observed with not so many micro-mirrors, as we show below.

In figure 4 we illustrate what happens when the number of micro-mirrors is small, a limit that is easier to implement experimentally. We take N = 7 micro-mirrors and show how localized structures can be supported by a single micro-mirror. We are able to write and erase such structures by an additional Gaussian optical injection at the desired position, and hence they are equivalent to the cavity solitons present in the continuous case. We find it remarkable that with such a small number of micro-mirrors the cavity solitons of the continuous model can still be recovered.

Once we have demonstrated that the predictions of the continuous model are robust for the discrete model, we pass to study a factor that may severely affect pattern formation, namely boundary effects. When a homogeneous solution is stable, any perturbation of it will be damped and will tend to disappear. However, the decay of the perturbation may be accompanied by spatial oscillations with characteristic spatial decay and a wavenumber [43], so that depending on the relation between these quantities and the extension of the array and the individual micro-mirror size, resonances may occur. As the boundaries of the mechanical array constitute a perturbation of the homogeneous solution, we can expect a strong influence of their presence in the patterns developed by the system.

We first illustrate the spatial oscillations in figure 5. For the parameters chosen (the system is driven by E₀² = 3.5), the system has a stable homogeneous solution Z̄ = |F̄|² = 2.73 for the continuous model [8], however the numerical simulation shows damped spatial oscillations starting from the two edges. We compare three different situations where the spatial extension of the mirror array is 40l_c, 80l_c and 160l_c (see figures 5(a)–(c) respectively). For all the other parameters fixed, the characteristic decay length is the same for the three situations. In particular, we obtain λrelax ≈ 8.5l_c by fitting the envelope of the damped oscillations of |F(x)|² to the function

$$y(x) = |\bar{F}|^2 + K e^{(|x|-x_{\text{max}})/\lambda_{\text{relax}}} \quad (19)$$

where K is an additional fitting parameter and x_max is the largest transverse position that we consider, that is, x ∈ [-x_max, x_max]. This fit is shown as a green-dashed line in figure 5. When the size of the system is comparable to this length, we observe the formation of a stable pattern, see figure 5(a). On the contrary when the system size becomes larger, then spatial oscillating perturbation is completely damped at the center where the homogeneous solution is re-established, see figure 5(c).
<figure><img src="image_8.png" /></figure>

As stated, because of the discreteness of the mirror array, the perturbation induced by the boundaries gives origin to another effect when it resonates with the spatial period of the mirror array. In this case we observe the emergence of a pattern of spatial period $\lambda_c = 2\pi / \text{Re}\{k_c\}$ when the size $a$ of the micro-mirrors is such that

$$a \approx \lambda_c/2 \qquad (20)$$

where $k_c$ is the critical wave vector (in general a complex number) for which one of the eigenvalues of the linear stability analysis is equal to zero. In practical terms, we find it by setting to zero equation (C12d) in [8]. In figure 6 we present the case of an array of mirrors of total size equal to 160l_c, when the system is driven in the upper homogeneous branch by a pump field $E_0^2 = 3.5$, for three different values of mirror size. In figure 6(a) we choose $N = 40$ so that $a = 4$, in figure 6(b) $N = 55$ so that $a \approx 2.9$, and in figure 6(c) $N = 80$ so that $a = 2$. The critical wavevector for the case considered is such that $\text{Re}[k_c] \approx 1.1$. As discussed in the previous paragraph, the boundary conditions induce a perturbation that is spatially modulated relaxing toward the homogeneous solution with some characteristic length. In the case we are considering there is enough room for these oscillations to relax on the homogeneous solution. Indeed in figures 6(a) and (c) the spatial period of the discrete array is very different from $\lambda_c/2 \approx 2.85$, while for (b) expression (20) is matched so that a periodical pattern emerges.

## 5. Connection with the standard OM array theory

In the limit where the intracavity optical field varies slowly enough with respect to the size $a$ of the micro-mirrors and the diffraction length $l_c$ is larger than $a$, we have about one optical mode locally interacting with each micro-mirror. In this case our model defined by equations (3) and (9) can be mapped onto the OM array suggested in [36] where ‘[...] a localized mechanical mode interacts with one laser-driven cavity (optical) mode [...] and where both photons and phonons can hop between neighboring sites’.

Then, we consider the optical field $A_j = A(\mathbf{r}_j)$ at each node $\mathbf{r}_j$ of the lattice defined by our mechanical array, and approximate the Laplacian in (3) by its finite-differences expression of (7),

$$\nabla_{\perp}^2 A \approx \frac{1}{a^2} \mathcal{L}_{\nu}[A_j]. \qquad (21)$$

As a consequence, equation (3) can be mapped onto the classical part of the quantum Langevin equation for the optical field $\alpha_j$ in [36]

$$\frac{d\alpha_j}{dt} = \left[-\frac{\bar{\kappa}}{2} + i\bar{\Delta} + ig_0\bar{q}_j\right]\alpha_j + i\frac{J}{\nu}\sum_{\langle l \rangle_j}\alpha_l - i\alpha_L \qquad (22)$$

and (9) can be mapped onto the classical part of the quantum Langevin equation for the mechanical field $\beta_j$ in [36]

$$\frac{d\bar{q}_j}{dt} = \bar{\Omega}\bar{p}_j, \qquad (23a)$$

$$\frac{d\bar{p}_j}{dt} = -\frac{\Gamma}{2}\bar{p}_j - \bar{\Omega}\bar{q}_j + 2g_0|\alpha_j|^2 + \frac{K}{\nu}\sum_{\langle l \rangle_j}\bar{q}_l \qquad (23b)$$
<figure><img src="image_9.png" /><figcaption>Figure 6. Normalized field intensity $|F(x)|^2$ (solid blue) in the steady-state as a function of the position, in a 1D system of finite size with $x \in [-80l_c, 80l_c]$ for (a) $N = 40$, (b) $N = 55$ and (c) $N = 80$. The parameters in this simulation are $M = 11$, $\gamma = 0.1$, $\Omega = 10$, $\Delta = -2.2$, $\rho = 1.13$ and $E_0 = \sqrt{3.5}$. The red horizontal dashed curve represent the homogeneous solution expected from the continuous model.</figcaption></figure>

where $\bar{q}_j = \beta_j + \beta_j^*$, $\bar{p}_j = -i(\beta_j - \beta_j^*)$ and

$$ \alpha_j = aA_j, \qquad (24a) $$

$$ \bar{q}_j = q_j/\Delta q_0, \qquad (24b) $$

$$ \bar{p}_j = m\dot{q}_j/\Delta p_0, \qquad (24c) $$

with $\Delta q_0 = \sqrt{\hbar/2m\Omega}$ and $\Delta p_0 = \hbar/\Delta q_0$ the zero point fluctuations. Then, we obtain a one-to-one correspondence with the dynamics of the OM arrays of [36], established by the following correspondence between the parameters of the two models.

$$ g_0 = G\Delta q_0, \qquad (25a) $$

$$ \bar{\kappa} = 2\gamma_c, \qquad (25b) $$

$$ \bar{\Gamma} = 2\gamma_m, \qquad (25c) $$

$$ \alpha_L = -i\gamma_c\mathcal{E}, \qquad (25d) $$

$$ \frac{K}{\nu} = \frac{\kappa_{\perp}\Delta q_0}{\Delta p_0}, \qquad (25e) $$

$$ \frac{J}{\nu} = \frac{\gamma_c l_c^2}{a^2}, \qquad (25f) $$

$$ \bar{\Delta} = \gamma_c \left( \Delta - \nu \frac{l_c^2}{a^2} \right), \qquad (25g) $$

$$ \bar{\Omega} = \sqrt{\Omega_m^2 + \nu\kappa_{\perp}/m}. \qquad (25h) $$

The above derivation, thus, allows the connection between our model and the standard OM array model [36] and, given that both discrete models are discrete versions of continuous models, the derivation also establishes the relation between our pattern forming model and the continuum optomechanics theory of [33].

Next we proceed to numerically illustrate the type of patterns appearing in this standard OM array limit. We have chosen a square array with $32 \times 32$ mirrors and a total area of $(32l_c)^2$ ($x, y \in [-16l_c, 16l_c]$), such that the area of each mirror is $a^2 = l_c^2$. This means that we are just at the limit where our model can be mapped onto the discrete model. For improving the validity of the correspondence between the two models (i.e. $a < l_c$) either one can increase the number of mirrors or reduce the size of the array. While the first option has, in principle, no limitation apart from having enough computational power or—experimentally—having a large array, the second option is physically limited by the fact that the total size of the array should provide enough room for the emergence of structures that have a characteristic linear size on the order of $l_c$. For an array area smaller than $(32l_c)^2$ we observed a reduced capacity of sustaining structures. Finally we have chosen the remaining parameters as in the previous section: $\gamma = 0.1$, $\Omega = 10$, $\Delta = -2.2$ and $\rho = 1.13$. The corresponding parameters of the OM array model, given by equations (22) and (5), are $K \approx 0.82\bar{\kappa}$, $J \approx 3.67\bar{\kappa}$, $\bar{\Delta} \approx -5.87\bar{\kappa}$, $\bar{\Omega} \approx 23.85\bar{\kappa}$, and $\bar{\Gamma} = 0.1\bar{\kappa}$. Similarly to our previous equations, the bare OM coupling $g_0$ can be removed from the equations by rescaling the variables, and hence again plays no role apart from a scaling factor.

By virtue of the model correspondence, the OM array presents the same curve of steady-state homogeneous solutions as the model of equations (3) and (9) (also of [8]), which is bistable for the chosen
<figure><img src="image_10.png" /></figure>

parameters. Interestingly, when the system is driven by an injection field $E_0^2 = 3.0$, the system spontaneously generates a 4 × 4 square pattern both in the optical and mechanical fields (the numerical simulations are again performed by using periodic boundary conditions and a top-hat super-Gaussian injection in two dimensions, with $\sigma_x = \sigma_y = 0.95 \times 16l_c$). This can be appreciated in figure 7(a), where we show the normalized mechanical displacement $z_j$ as a function of the OM-cell index $j$. When the injection field is reduced to $E_0^2 = 2.35$, the optical and mechanical fields spatially localize into a cavity soliton as the steady state. This is shown in figure 7(b), where it is appreciated that the soliton involves only 5 elementary cells. This solution coexists with the homogeneous solution of the lower branch of the steady state curve. As a consequence it is possible to write several localized structures for $E_0^2 = 2.35$. This is shown in figure 7(c) where two solitons have been written at $j = (9, 9)$ and $j = (23, 23)$. Differently from examples shown in figure 4, in the case of OM arrays it is not possible to write a soliton for a too small number of cells (for example $N = 7$). This is because in the situation of figure 4 we are in the $l_c < a$ limit, so that structures smaller than one mirror can be sustained by the optical field, allowing for a small size of the array. In contrast, in the case of OM array of this section, we are in the $l_c \ge a$ limit, for which the characteristic size of self-localized structures is larger than the size of the micro-mirrors, as clearly appreciated in figure 7. In terms of the model equations (23a) and (23b) this is equivalent to asking $J/\nu > \gamma_c$. This limitation could be circumvented by considering a generalization of the model in [36] where each elementary cell is an OM cavity which is degenerate for at least two optical modes.

## 6. Conclusions

In this paper we have proposed an architecture for an OM cavity that allows for the generation of dissipative structures. The device consists of an OM cavity with an oscillating end-mirror formed by an array of weakly-coupled micro-mirrors. This configuration fulfills the basic requirement necessary for OM dissipative structures: the existence of a homogeneous mechanical transverse mode [8]. This proposal offers then an alternative to our previous one [8] that consisted in mounting a flexible mirror on a large aspect ratio frame (with dimensions $L_x \gg L_y$).

The model we are proposing coincides mathematically with that studied in [8] in the limit of large density of micro-mirrors, making all its analytical and numerical results applicable in such limit. We have numerically shown that this is also the case, more qualitatively, when the number of micro-mirrors is not very large, and have found that with a relatively small number of elements there exist solutions reminiscent of the continuous-case cavity solitons. More concretely, we have found that a discrete model consisting of $N \approx 10$ micro-mirrors with size $a \approx 4l_c$ is enough to observe localized structures exactly as predicted by the continuous limit model. Periodic patterns may require a larger number of micro-mirrors depending on their periodicity, but in any case they should still be well captured with a reasonable number of these (say $N < 100$). In the second part of the paper we have also connected our model with the discrete OM arrays model recently put forward in [36]. The connection appears through the discretized version in [36] meaning that our OM cavity model (with micro-structured mirror) is equivalent to an OM array when the diffraction length is larger than the individual micro-mirror size. Hence, implementing pattern formation in a microstructured OM-cavity is a way of implementing OM arrays and, with more generality, pattern formation is a robust implementation of continuum optomechanics. We hope that our work be useful in the experimental search of dissipative structures in OM devices.
Acknowledgments

We thank Chiara Molinelli for useful discussions in the initial brainstorming phase. This work was funded by Spanish Ministerio de Ciencia, Innovación y Universidades—Agencia Estatal de Investigación, and European Union FEDER (projects FIS2014-60715-P and FIS2017-89988-P). CNB acknowledges additional support from a Shanghai talent program and from the Shanghai Municipal Science and Technology Major Project (Grant No. 2019SHZDZX01). GP acknowledges support from Lille University within the framework ‘Internationalisation de la recherche 2019 - collaboration bilatérales’.

Appendix A. Interpretation of the optical field amplitude and radiation pressure

In the main text we wrote the electric field propagating to the right as

$$
\mathbf{E}_{+}(z, \mathbf{r}, t) = i\mathbf{x}\mathcal{V}A_{+}(z, \mathbf{r}, t) e^{ik_{L}z - i\omega_{L}t} + \text{c.c.}, \qquad (A.1)
$$

where here we include the polarization of the electric field (defining the x direction with corresponding unit vector **x**), which was omitted in the main text for simplicity. The aim of this section is to explain how the choice **V** = √ħω<sub>c</sub>/4ε<sub>0</sub>L characteristic of quantum optics allows us to give a simple interpretation to the amplitude A<sub>+</sub> (and similarly for A<sub>-</sub> and A<sub>inj</sub>), as well as writing the expression for the radiation pressure exerted by **E**<sub>+</sub> in terms of this amplitude.

Let us first remind that, within the paraxial approximation, the magnetic field associated to (A.1) can be
written as (y is the unit vector in the y direction)

$$
\mathbf{B}_{+}(z, \mathbf{r}, t) = i\mathbf{y}c^{-1}\mathcal{V}A_{+}(z, \mathbf{r}, t) e^{ik_{L}z - i\omega_{L}t} + \text{c.c.}; \quad (A.2)
$$

the corresponding Poynting vector is then written as (z is the unit vector in the z direction)

$$
\mathbf{S}_{+} = \frac{1}{\mu_0} \mathbf{E}_{+} \times \mathbf{B}_{+} = - \frac{\mathcal{V}^2 \mathbf{z}}{\mu_0 c} (A_+ e^{ik_L z - i\omega_L t} - \text{c.c.})^2, \quad (A.3)
$$

whose magnitude averaged over an optical cycle

$$
\langle S_+ \rangle|_{z=L} = \frac{2\pi}{\omega_L} \int_{t-\pi/\omega_L}^{t+\pi/\omega_L} d\tau |S_+(L, \mathbf{r}, \tau)| \approx \frac{2\mathcal{V}^2}{\mu_0 c} |A(\mathbf{r}, t)|^2,
$$

provides the instantaneous measurable power impinging point **r** of the mirror located at *z* = *L* per unit area (irradiance). Note that we have made use of the slowly time-varying nature of the amplitude, and remember that we defined *A* (**r**, *t*) = *A*<sub>+</sub>(*L*, **r**, *t*) in the main text. Now it is customary in quantum optics to take *V* = √ħω<sub>c</sub>/4ε<sub>0</sub>*L* so that

$$
|A(\mathbf{r}, t)|^2 = \frac{t_c \langle S_+ \rangle|_{z=L}}{\hbar \omega_c}, \qquad (A.4)
$$

can be interpreted as the number of photons per unit area which arrive to point **r** of the mirror during a
round-trip (*t*<sub>c</sub> = 2*L*/*c* is the cavity round-trip time). With this choice, the theory is quantized by
interpreting this amplitude as an operator satisfying equal-time commutation relations

$$
[\hat{A}(\mathbf{r}, t), \hat{A}^{\dagger}(\mathbf{r}', t)] = \delta(\mathbf{r} - \mathbf{r}') \text{ and } [\hat{A}(\mathbf{r}, t), \hat{A}(\mathbf{r}', t)] = 0.
$$

From the Poynting vector, we can get the radiation pressure exerted onto a point **r** of the flexible mirror as $P(\mathbf{r}, t) = \langle S_+ \rangle|_{z=L}/c$; this is a quantity of fundamental relevance to our work, and in our case takes the particular expression

$$
P(\mathbf{r}, t) = \frac{\hbar k_c}{t_c} |A(\mathbf{r}, t)|^2; \qquad (A.5)
$$

given our interpretation of $|A(\mathbf{r}, t)|^2$, this coincides precisely with the total momentum (momentum per photon × number of photons) hitting point **r** of the flexible mirror per unit time and area.

**Appendix B. Derivation of the light field equation**

Here we derive equation (3) of the main text. To this aim we use the approach of references [8, 38], which
consists in propagating the complex amplitudes A_± (z, **r**, t) along a full cavity round-trip. Assuming that
they are slowly varying in space and time, they satisfy the paraxial wave equation

$$
(\partial_z \pm c^{-1} \partial_t) A_{\pm} = \pm \frac{i}{2k_L} \nabla_{\perp}^2 A_{\pm}. \qquad (B.1)
$$
Given the amplitude $A_+(z = L, \mathbf{r}, t)$, after reflection on the flexible mirror we get

$$A_{-}(L, \mathbf{r}, t) e^{-ik_{L}L} = -A_{+}(L, \mathbf{r}, t) e^{ik_{L}[L+2Q(\mathbf{r},t)]}, \quad (B.2)$$

where $Q(\mathbf{r}, t)$ represents the displacement of the mirror from its rest position ($Q = 0$ at rest). The amplitude $A_{-}(L, \mathbf{r}, t)$ propagates from $z = L$ to $z = 0$ giving rise to a new amplitude

$$A_{-}(0, \mathbf{r}, t + \frac{1}{2}t_c) = U_{L}A_{-}(L, \mathbf{r}, t), \quad (B.3)$$

where

$$U_{L} = \exp \left[ i(L/2k_{L})\nabla^{2} \right], \quad (B.4)$$

is the paraxial propagation operator in free space. After reflection onto the coupling mirror, a new amplitude

$$A_{+}(0, \mathbf{r}, t + \frac{1}{2}t_c) = -\sqrt{R}A_{-}(0, \mathbf{r}, t + \frac{1}{2}t_c) + \sqrt{T}A_{\text{inj}}(0, \mathbf{r}, t + \frac{1}{2}t_c), \quad (B.5)$$

is got, with $R$ and $T$ the reflectivity and transmissivity factors of the coupling mirror, respectively ($R + T = 1$ is assumed: lossless mirror). Finally, propagation from $z = 0$ to $z = L$ yields $A_{+}(L, \mathbf{r}, t + t_c) = U_{L}A_{+}(0, \mathbf{r}, t + \frac{1}{2}t_c)$. Adding all parts together one gets

$$A(\mathbf{r}, t + t_c) = \sqrt{R} e^{2ik_L L} U_L^2 \exp[2ik_L Q(\mathbf{r}, t)] A(\mathbf{r}, t) + \sqrt{T} A_{\text{inj}}(L, \mathbf{r}, t + t_c), \quad (B.6)$$

where we used $U_{L}A_{\text{inj}}(0, \mathbf{r}, t + \frac{1}{2}t_c) = A_{\text{inj}}(L, \mathbf{r}, t + t_c)$. We now take into account that $R \to 1$ (equivalently, $T \to 0$) so that $\sqrt{R} = \sqrt{1 - T} \to 1 - T/2$. Next we assume that light is almost resonant with the cavity, specifically we impose that $2(\omega_L - \omega_c) L/c = \delta$ is of order $T$, where $\omega_c$ is the cavity longitudinal mode frequency (hence $\omega_c = m\pi c/L$, $m \in \mathbb{N}$) closest to $\omega_L$, what allows approximating $\exp(2ik_L L) = \exp(2i\omega_L L/c) \approx 1 + i\delta$. We assume as well that $k_L Q(\mathbf{r}, t)$ is of order $T$ (the mirror displacement/deformations are much smaller than the optical wavelength), so that $\exp[2ik_L Q(\mathbf{r}, t)] \approx 1 + 2ik_L Q(\mathbf{r}, t)$. Similarly we assume that the effect of diffraction is small (this implies that both mirrors must be sufficiently close each other, either physically or by means of lenses) so that we can expand $U_L^2 \approx 1 + i(L/k_L)\nabla_{\perp}^2$. All these assumptions imply that the overall variation of $A$ between consecutive round-trips is very small and then one can approximate $\partial_t A$ by $[A(\mathbf{r}, t + t_c) - A(\mathbf{r}, t)] t_c^{-1}$. With all these approximations we get, to the lowest nontrivial order,

$$\partial_t A(\mathbf{r}, t) = \gamma_c \left( -1 + i\Delta + i l_c^2 \nabla_{\perp}^2 + i \frac{4k_L}{T} Q \right) A + \gamma_c \mathcal{E}, \quad (B.7)$$

where all the parameters are defined in the main text; this is precisely (3), and it is the same light-field equation we derived in the linear-coupling model of [8].

## References

[1] Staliunas K and Sánchez-Morcillo V J 2003 *Transverse Patterns in Nonlinear Optical Resonators* (Berlin: Springer)

[2] Mandel P 2005 *Theoretical Problems in Cavity Nonlinear Optics* (Cambridge: Cambridge University Press)

[3] Cross M C and Hohenberg P C 1993 *Rev. Mod. Phys.* **65** 851

[4] Firth W J and Weiss C O 2002 *Opt. Photonics News* **13** 54

[5] Barland S et al 2002 *Nature* **419** 699

[6] Esteban-Martín A, Taranenko V B, Roldán E and de Valcárcel G J 2005 *Opt. Exp.* **13** 3631

[7] Ackemann T, Firth W and Oppo G-L 2009 *Adv. At. Mol. Opt. Phys.* **57** 323

[8] Ruiz-Rivas J, Navarrete-Benlloch C, Patera G, Roldán E and de Valcárcel G J 2015 *Phys. Rev. A* **93** 033850

[9] Aspelmeyer M, Kippenberg T J and Marquardt F 2014 *Rev. Mod. Phys.* **86** 1391

[10] Tesio E, Robb G R M, Ackemann T, Firth W J and Oppo G-L 2014 *Phys. Rev. Lett.* **112** 043901

[11] Labeyrie G, Tesio E, Gomes P M, Oppo G-L, Firth W J, Robb G R M, Arnold A S, Kaiser R and Ackemann T 2014 *Nat. Photonics* **8** 321

[12] Baio G, Robb G R M, Yao A M and Oppo G-L 2020 *Phys. Rev. Res.* **2** 023126

[13] Arcizet O, Cohandon P-F, Briant T and Heidemann A 2006 *Nature* **444** 71

[14] Gigan S et al 2006 *Nature* **444** 67

[15] Kleckner D and Bouwmeester D 2006 *Nature* **444** 75

[16] Weis S, Rivière R, Deléglise S, Gavartin E, Arcizet O, Schliesser A and Kippenberg T J 2010 *Science* **330** 1520

[17] Safavi-Naeini A H, Mayer Alegre T P, Chan J, Eichenfield M, Winger M, Lin Q, Hill J T, Chang D E and Painter O 2011 *Nature* **472** 69

[18] Brooks D W C, Botter T, Schreppler S, Purdy T P, Brahms N and Stamper-Kurn D M 2010 *Nature* **488** 476

[19] Safavi-Naeini A H, Gröblacher S, Hill J T, Chan J, Aspelmeyer M and Painter O 2013 *Nature* **500** 185

[20] Purdy T P, Yu P-L, Peterson R W, Kampel N S and Regal C A 2013 *Phys. Rev. X* **3** 031012

[21] Verhagen E, Deléglise S, Weis S, Schliesser A and Kippenberg T J 2012 *Nature* **482** 63
[22] Lugiato L A, Brambilla M and Gatti A 1999 *Adv. At. Mol. Opt. Phys.* **40** 229

[23] Pérez-Arjona I, Roldán E and de Valcárcel G J 2006 *Europhys. Lett.* **74** 247

[24] Pérez-Arjona I, Roldán E and de Valcárcel G J 2007 *Phys. Rev. A* **75** 063802

[25] Navarrete-Benlloch C, Roldán E and de Valcárcel G J 2008 *Phys. Rev. Lett.* **100** 203601

[26] Lugiato L A and Grynberg G 1995 *Europhys. Lett.* **29** 675

[27] Lugiato L A and Gatti A 1993 *Phys. Rev. Lett.* **70** 3868

[28] Gatti A and Lugiato L A 1995 *Phys. Rev. A* **52** 1675

[29] Santagiustina M, Colet P, San Miguel M and Walgraef D 1997 *Phys. Rev. Lett.* **79** 3633

[30] Vaupel M, Maitre A and Fabre C 1999 *Phys. Rev. Lett.* **83** 5278

[31] Safavi-Naeini A H, Van Thourhout D, Baets R and Van Laer R 2019 *Optica* **6** 213

[32] Liu Q, Li H and Li M 2019 *Optica* **6** 778

[33] Rakich P and Marquardt F 2018 *New J. Phys.* **20** 045005

[34] Heinrich G, Ludwig M, Qian J, Kubala B and Marquardt F 2011 *Phys. Rev. Lett.* **107** 043603

[35] Holmes C A, Meaney C P and Milburn G J 2012 *Phys. Rev. E* **85** 066203

[36] Ludwig M and Marquardt F 2013 *Phys. Rev. Lett.* **111** 073603

[37] Cernotík O, Mahmoodian S and Hammerer K 2018 *Phys. Rev. Lett.* **121** 110506

[38] Kolpakov S, Esteban-Martín A, Silva F, García J, Staliunas K and de Valcárcel G J 2008 *Phys. Rev. Lett.* **101** 254101

[39] Braginsky V B, Strigin S E and Vyatchanin S P 2001 *Phys. Lett. A* **287** 331

[40] Braginsky V B, Strigin S E and Vyatchanin S P 2002 *Phys. Lett. A* **305** 111

[41] Chen X *et al* 2015 *Phys. Rev. A* **91** 033832

[42] Evans M *et al* 2015 *Phys. Rev. Lett.* **114** 161102

[43] Sánchez-Morcillo V J and Staliunas K 1999 *Phys. Rev. E* **60** 6153