Optical metamaterials and metasurfaces, which emerged in the course of the last few decades, have revolutionized our understanding of light and light–matter interaction. While solid materials are naturally employed as key building elements for construction of optical metamaterials mainly due to their structural stability, practically no attention was given to study of liquid-made optical two-dimensional (2-D) metasurfaces and the underlying interaction regimes between surface optical modes and liquids. We theoretically demonstrate that surface plasmon polaritons and slab waveguide modes that propagate within a thin liquid dielectric film trigger optical self-induced interaction facilitated by surface tension effects, which leads to the formation of 2-D optical liquid-made lattices/metasurfaces with tunable symmetry and can be leveraged for tuning of lasing modes. Furthermore, we show that the symmetry breaking of the 2-D optical liquid lattice leads to phase transition and tuning of its topological properties, which allows the formation, destruction, and movement of Dirac-points in the k-space. Our results indicate that optical liquid lattices support extremely low lasing threshold relative to solid dielectric films and have the potential to serve as configurable analogous computation platform.

Metamaterials are composite man-made or natural materials that possess emergent optical properties, which stem from specific spatial arrangements of the constituent subwavelength basic units^{1}^{,}^{2} and lead to numerous effects, such as formation of bandgap in dielectric photonic crystals,^{3}^{,}^{4} suppression of plasmon frequency in metallic metamaterials,^{5} control over the radiation dynamics of embedded active materials^{6}^{,}^{7} (see also Ref. 8 and references therein), wave front control using thin elements,^{9} polarization and phase control in both transmission^{10} and reflection^{11} modes, and control of light properties in low-loss and high-index dielectric resonant Mie nanoparticles.^{12} Dynamical tuning of optical metamaterial properties is particularly appealing as it allows one to study new regimes and effects of light–matter interaction, holds promise for future metamaterial-based devices with functionalities achieved by structuring matter on the subwavelength scale,^{13} and may be also of interest as platforms to simulate many body quantum effects, which are challenging to realize in real quantum systems.^{14} In particular, numerous studies, partly fueled by seminal advances in condensed matter physics, such as the discovery of graphene,^{15} have revealed analogies between propagation of light in photonic crystals to dynamics of relativistic Dirac fermions near Dirac points in crystals and dynamics of electrons in topological insulators. In the case of two-dimensional (2-D) photonic structures, Dirac cones have attracted significant interest because of the existence of robust surface states due to the breaking of parity and time-reversal symmetry,^{16} and intriguing transport properties such as pseudodiffusive transmittance,^{17} persistence of the Klein effect,^{18} and breakdown of conical diffraction due to symmetry breaking of hexagonal symmetry^{19} or due to nonlinear interactions.^{20}^{,}^{21} Furthermore, 2-D periodic structures play an important role in realizing lasing effects in so-called distributed feedback (DFB) structures, which provide continuous coherent backscattering from the periodic structures without mirrors; these were originally proposed in one-dimensional (1-D)^{22}^{,}^{23} and later were extended to 2-D systems enabling lasing amplification of waveguide modes in photonic crystals^{24} and of surface plasmon polariton (SPP) oscillations in metal structures,^{25} and more recently also in thin polymer membranes.^{26} Since metasurfaces in general and 2-D periodic structures in particular are conventionally constructed using solid metals and high index dielectrics, their tuning properties are constrained by physical properties, such as carrier^{27} and material density, which can be manipulated by external electrical, magnetic, acoustic, and temperature fields,^{28} and also by mechanical stretching of the elastic DFB structure, which affects the resonant lasing frequencies.^{29}

Liquids on the other hand provide an attractive platform to introduce significant changes to the optical properties of metamaterials and metasurfaces, due to the liquids’ capability to induce relatively large dynamical changes of their dielectric function, which stems from their adaptive property to fill microchannels of desired shape,^{30} compliance under external stimuli, and ability to sustain changes of their physicochemical properties. Prominent examples include light-induced collective orientation of liquid crystal molecules,^{31} pressure-induced control of lasing frequency in a 1-D array of droplets,^{32} magnetic induced ferrofluid-based hyperbolic metamaterial,^{33} and chemical composition-induced changes enabling lasing frequency tuning in 1-D systems^{34} and in 2-D photonic crystals.^{35} Furthermore, recent studies introduced a novel interaction between SPPs and a thin liquid dielectric (TLD) film due to geometrical changes of the gas–fluid (or fluid–fluid) interface facilitated by the thermocapillary (TC) effect: theoretical study of self-induced focusing and defocusing effects of propagating SPPs due to nonlocal interaction^{36} (where refractive index changes extend beyond the regions of maximal optical intensity) as well as formation of an optical liquid lattice of a fixed square symmetry, and experimental demonstration of TC-assisted optical tuning of surface plasmon resonance coupling angle.^{37} While these two studies demonstrated a significant coupling between the topography of the gas–fluid interface and propagation of SPP modes, these studies did not leverage surface tension nor surface optical modes to study bandstructure and lasing modes tuning due to optically induced changes of the liquid lattice symmetry.

In this work, we theoretically demonstrate that SPPs or slab waveguide (WG) modes propagating in a TLD film, which is thinner than the penetration depth of the corresponding surface optical mode in the direction normal to the film’s surface, lead to a nonlocal and nonlinear response of the corresponding dielectric function due to optically driven surface tension effects. In particular, we take advantage of the surface tension dependence on local physicochemical conditions^{38} to show that the optically induced TC^{39} or solutocapillary (SC) effects,^{40} which stem from temperature or chemical concentration gradients, respectively, lead to self-induced changes of the dielectric function. In particular, since both surface tension effects are accompanied by flows and deformation of the TLD film, the latter lead to self-induced changes of the dielectric function that are coupled back to the propagation conditions of the surface optical mode. Importantly, the fact that both SPP and WG modes can propagate within thin dielectric films, which are thinner than the optical penetration depth of the corresponding surface optical mode into the domain outside the TLD film, facilitates a significant coupling between these surface optical modes and changes of liquid’s topography. In particular, the SPPs and slab WG modes are both guaranteed to propagate in arbitrarily thin liquid films; the former due to the fundamental capability of metal–dielectric interfaces to support SPPs,^{41} whereas the latter can be described by the transverse resonance condition of the fundamental mode.^{42} Employing this interaction for the case of interfering surface optical waves leads to a self-induced optical liquid lattice of tunable symmetry and bandstructure, which can be tuned by changing the relative propagation directions and the amplitudes of the interfering SPP or WG modes. Furthermore, applying bandstructure tuning in the case, the TLD film or the dielectric substrate admits gain properties leads to configurable DFB mechanism with gain and/or index modulation that can control the threshold condition and the corresponding lasing frequency. In particular, dielectric substrate with gain and TLD film without gain leads to index modulation, whereas the complementary case of TLD film with gain supports both index and gain modulation, as shortly discussed below.

Figure 1 presents a periodic deformation of a TLD film driven by SPPs propagating on a metal–fluid interface and forming a plasmonic liquid lattice, whereas Fig. 1(b) shows a periodic deformation of a pair of gas–fluid interfaces present in a symmetric liquid slab WG due to propagating WG modes forming a suspended photonic liquid lattice without metals. In practice, the suspended optically thin liquid film can be realized either by bracketing it with an immiscible liquid, which admits a distinct refractive index^{43} or by stabilizing liquid film with surfactants leading to an optically thin and stable liquid sheet with a pair of parallel gas–fluid interfaces (see also Ref. 44 for a recent experiments of optical guiding in soap films). Figure 1(c) shows a thin solid dielectric substrate, which hosts a TLD film on both of its sides, leading to a symmetric slab WG forming a supported photonic liquid lattice. For simplicity through the work, we consider the case where the refractive index of the solid dielectric components and of the TLD film is equal to minimize the impedance mismatch, which could be realized by employing a silicone oil and silica of refractive indices 1.4 and 1.45 (under 785-nm illumination), respectively. Under this assumption, the system shown in Fig. 1(c) allows in principle realization of the optically similar setup shown in Fig. 1(b), without employing surfactants, provided the film deformation is sufficiently small. To the best of our knowledge, the prospect of employing the prominent interaction of surface optical modes and TLD film, in order to induce and control the symmetry of 2-D optical liquid lattices, and furthermore to leverage it for lasing modes and of topological properties tuning, has not been explored to date.

This paper is structured as follows. First, we describe the nonlinear self-induced interaction between surface optical modes and the TLD film driven by either the TC or SC effects and derive the underlying complex nonlocal Ginzburg–Landau equation that governs the dynamics of the corresponding envelope function. We then solve the TLD film equation and describe tuning of the optical liquid lattices and of the corresponding symmetries due to changes of the propagation directions and amplitudes of the surface optical modes. As an illustrative example, we first consider the simpler 1-D case and analyze the threshold change and lasing frequency tuning, due to self-induced and nonself-induced amplitude modulation of the formed liquid lattice. Afterward, we present numerical simulation results of bandstructure tuning due to the breaking of 2-D hexagonal and square symmetries and demonstrate formation of Dirac points in liquid lattices with a broken hexagonal symmetry. Finally, we demonstrate that symmetry changes of 2-D optical liquid lattices leads to tuning of the corresponding lasing frequencies and the corresponding emission directions.

## 1.

## Results

## 1.1.

### Light-Induced Interaction between Surface Optical Modes and a TLD Film

The set of coupled governing equations that describes light–fluid interaction due to TLD film thickness changes includes Maxwell equations, Navier–Stokes equations for an incompressible Newtonian fluid, balance condition between viscous and surface tension stresses on the gas–fluid (or fluid–fluid) interface, and heat/mass-transport equations depending on the specific light-induced mechanism, which triggers local changes of the surface tension. We employ heat-transport and mass-transport equations, relevant for the TC effect and SC effect, respectively, which constitute the coupling mechanism between light propagation and dynamics of TLD film. Specifically, light-induced changes of the surface tension lead to deformation of the TLD film, $\eta ({\overrightarrow{r}}_{\parallel},t)$, which in turn is coupled back to the propagation of light due to the associated spatial changes of the dielectric function, $\mathrm{\Delta}{\u03f5}_{D}$, which in the leading order of $\eta /{h}_{0}$ is given as^{36}

## Eq. (1)

$$\mathrm{\Delta}{\u03f5}_{D}({\overrightarrow{r}}_{\parallel},t)=b\eta ({\overrightarrow{r}}_{\parallel},t)/{h}_{0}.$$## Eq. (2)

$$\eta ({\overrightarrow{r}}_{\parallel},t)/{h}_{0}=-\mathrm{M}\int \mathrm{d}{\overrightarrow{r}}_{\parallel}^{\prime}\mathrm{d}{t}^{\prime}\frac{1}{{\tau}_{\mathrm{th}}}{G}_{l}({\overrightarrow{r}}_{\parallel}-{\overrightarrow{r}}_{\parallel}^{\prime},t-{t}^{\prime})I({\overrightarrow{r}}_{\parallel}^{\prime},{t}^{\prime})/{I}_{0},$$^{45}and optical power of $\sim 120\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mW}$ concentrated in a region of few mm. In case of TLD film deformation driven by surface optical modes, which admit significantly smaller mode volume (of several hundreds of nm at most), we expect much lower powers needed to generate TLD film deformation. Note that in TLD film, nanoparticle transport is coupled to both temperature field and to TC flows, which require additional equation in our model; also, local changes of nanoparticles concentration are expected to introduce local changes of shear viscosity,

^{46}which are beyond the scope of this work.

For the case of SC flows, we focus on cis–trans transformation that can be described by the simple phenomenological two-state model $A\underset{{k}_{\mathrm{off}}}{\overset{{k}_{\mathrm{on}}}{\rightleftharpoons}}B$, where ${c}_{A,B}$ stands for the molar concentration of the corresponding photoactive molecule in the cis/trans states $A$ and $B$, respectively. Assuming linear dependence of the surface tension $\sigma $, with respect to small changes of ${c}_{A,B}$, and furthermore assuming for simplicity that both species admit identical molecular diffusion coefficients $D$ and similar affinity to the gas–fluid interface, leads to the following nonlocal relation between the deformation $\eta $ and the optical power $I$, which is similar to the relation Eq. (2) relevant to the TC case discussed above (see the Supplemental Material for derivation):

## Eq. (3)

$$\eta ({\overrightarrow{r}}_{\parallel},t)/{h}_{0}=-{\mathrm{M}}_{c}\int \mathrm{d}{\overrightarrow{r}}_{\parallel}^{\prime}\mathrm{d}{t}^{\prime}\frac{1}{{\tau}_{\mathrm{on}}}{G}_{l}({\overrightarrow{r}}_{\parallel}-{\overrightarrow{r}}_{\parallel}^{\prime},t-{t}^{\prime})I({\overrightarrow{r}}_{\parallel}^{\prime},{t}^{\prime})/{I}_{0},$$Employing a perturbative expansion for the SPP and slab WG modes (see the Supplemental Material), which incorporates the effects of diffraction and gain, yields the following complex nonlocal Ginzburg–Landau (GL) equation (see the Supplemental Material for derivation):

## Eq. (4)

$$2i{\beta}_{0}\frac{\partial A}{\partial z}+\frac{{\partial}^{2}A}{\partial {y}^{2}}-V[A({\overrightarrow{r}}_{\parallel})]A+i\mathrm{\Gamma}A=0,\phantom{\rule{0ex}{0ex}}V[A({\overrightarrow{r}}_{\parallel})]\equiv -\mathrm{\Delta}{\epsilon}_{d}[A({\overrightarrow{r}}_{\parallel})]=-{\chi}_{\mathrm{TC}/\mathrm{SC}}\int \mathrm{d}{\overrightarrow{r}}_{\parallel}^{\prime}{G}_{l}({\overrightarrow{r}}_{\parallel},{\overrightarrow{r}}_{\parallel}^{\prime}){|A({\overrightarrow{r}}_{\parallel}^{\prime})|}^{2}.$$^{48}as well as nanofocusing of SPPs in tapered plasmonic waveguides.

^{49}The nonlocal integral term, $V[A({\overrightarrow{r}}_{\parallel})]$, in Eq. (4) represents the self-induced potential due to refractive index changes that stem from liquid deformation extending beyond the region of maximal optical intensity. For the case of a TLD film interacting with an SPP mode, the constant ${\chi}_{\mathrm{TC}/\mathrm{SC}}=fb\mathrm{M}$ is a product of three dimensionless numbers incorporating the effects of SPP enhancement, depth averaging, and TC or SC effects, whereas for the case of a TLD film interacting with slab WG modes, the constant is ${\chi}_{\mathrm{TC}/\mathrm{SC}}=b\mathrm{M}$. Importantly, the case where a TLD film introduces nonlocal loss or gain is incorporated through a complex constant $b$, where its positive/negative imaginary part corresponds to effects of loss/gain (i.e., damped/amplified modes), respectively. The term, $i\mathrm{\Gamma}$, on the other hand represents local loss or gain, which are not related to refractive index changes due to TLD film deformation. For instance, local change of the dielectric function due to gain changes is given by $i\mathrm{\Gamma}={\epsilon}_{m}^{2}\delta {\epsilon}_{d}/{({\epsilon}_{d0}+{\epsilon}_{m})}^{2}$, where $\delta {\epsilon}_{d}$ is the optically induced dielectric function change in the gain media. In dye solutions, the latter may stem either from intensity-dependent modulations, which affect the population of the electronic states of the dye molecules or from temperature gradients formed by the heat generated by relaxation of the excited molecules (see Ref. 50 and references within).

## 1.2.

### One-Dimensional and Two-Dimensional Optical Liquid Lattices with Tunable Symmetry

Next we show that by controlling the interference pattern of surface optical modes, it is possible to form optical liquid lattices directly from a TLD film and tune their symmetries. Consider an interference of $N$-plane waves of real-valued amplitudes ${a}_{n}$ that propagate in the $y-z$ plane along the directions formed by angles ${\theta}_{n}$, which are measured relative to the positive direction of the $y$ axis. The corresponding wave vectors are given by ${\overrightarrow{k}}_{n}={k}_{0}[\mathrm{cos}({\theta}_{n}),\mathrm{sin}({\theta}_{n})]$, where ${k}_{0}=2\pi /\lambda $ and $\lambda $ are the effective wavelength that is set by the relevant frequency and the corresponding dispersion relation; i.e., for SPP or WG mode, the wavelength is given by $\lambda ={\lambda}_{0}/n$, where ${\lambda}_{0}$ is the vacuum wavelength and $n$ is the corresponding effective refractive index. The resultant optical intensity $I$ due to interference of the plane waves ${a}_{n}{e}^{i{\overrightarrow{k}}_{n}\xb7{\overrightarrow{r}}_{\parallel}}$ is then given as

## Eq. (5)

$$\begin{array}{c}I={\left|\sum _{n=0}^{N-1}{a}_{n}{e}^{i{\overrightarrow{k}}_{n}\xb7{\overrightarrow{r}}_{\parallel}}\right|}^{2}=\end{array}\phantom{\rule{0ex}{0ex}}{I}_{N}+2\sum _{n\ne m}^{N-1}{a}_{n}{a}_{m}\mathrm{cos}[({\overrightarrow{k}}_{n}-{\overrightarrow{k}}_{m})\xb7{\overrightarrow{r}}_{\parallel}];\phantom{\rule[-0.0ex]{2em}{0.0ex}}{I}_{N}\equiv \sum _{n=0}^{N-1}{|{a}_{n}|}^{2}.$$First, consider the simplest $N=2$ case that leads to a 1-D optical liquid lattice. The intensity distribution of two optical plane waves of equal intensity ${I}_{0}$ propagating along the $y-z$ plane with wave vectors ${\overrightarrow{k}}_{\pm}=k[\pm \mathrm{cos}(\theta ),\mathrm{sin}(\theta )]$ [see Fig. 2(a)] is given by $2{I}_{0}\{1+\mathrm{cos}[K(\theta )]x\}$, which upon inserting into Eq. (2) yields the corresponding 1-D deformation of the TLD film:

## Eq. (6)

$$\frac{\eta}{{h}_{0}}=-\frac{2\mathrm{M}}{{\tau}_{th}{\lambda}_{\overline{n}}}\mathrm{cos}[K(\theta )x];K(\theta )=2k\mathrm{cos}(\theta );\mathrm{\Lambda}(\theta )\equiv \frac{2\pi}{K(\theta )}=\frac{\pi}{k\mathrm{cos}(\theta )},$$Next, we consider the case $N=3$, which leads to a 2-D symmetric liquid lattice, where the relative amplitudes of the three beams satisfy ${a}_{1}={a}_{2}={a}_{3}/q\equiv a$, where $a$ and $q$ are real constants and assume the corresponding propagation directions are given by ${\widehat{n}}_{1}=[\mathrm{cos}(\theta ),-\mathrm{sin}(\theta )]$, ${\widehat{n}}_{2}=-[\mathrm{cos}(\theta ),\mathrm{sin}(\theta )]$, ${\widehat{n}}_{3}=(\mathrm{0,1})$ (which can be also written as ${\theta}_{n}=(n+1)\pi +(-1{)}^{n}\theta $ for $n=\mathrm{0,1}$ and ${\theta}_{2}=\pi /2$). Inserting the corresponding intensity distribution (see the Supplemental Material) into Eq. (2) yields the following TLD film deformation:

## Eq. (7)

$$\begin{array}{c}\frac{\eta ({\overrightarrow{r}}_{\parallel},t)}{{h}_{0}}=-\frac{{\eta}_{\mathrm{max}}}{3}\left[\mathrm{cos}\right(\frac{4\pi y}{\mathrm{\Delta}y})+2\mathrm{cos}(\frac{2\pi y}{\mathrm{\Delta}y}\left)\mathrm{cos}\right(\frac{2\pi z}{\mathrm{\Delta}z}\left)\right].\end{array}$$## Eq. (8)

$$\mathrm{\Delta}y=\frac{\lambda}{\mathrm{cos}(\theta )};\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{\Delta}z=\frac{\lambda}{1+\mathrm{sin}(\theta )},$$## Eq. (9)

$${\eta}_{\mathrm{max}}\equiv (3\overline{q}\mathrm{M})/({\lambda}_{\overline{n},0}{\tau}_{\mathrm{th}});\phantom{\rule[-0.0ex]{1em}{0.0ex}}\overline{q}=2{\lambda}_{2\overline{n},0}/{\lambda}_{\overline{m},\overline{n}},$$Expanding the deformation given by Eq. (7) up to a second-order in Taylor series around one of the maximum points [e.g., $(\mathrm{\Delta}y/\mathrm{3,0})$] yields a harmonic oscillator potential that admits a degree of anisotropy given by $3(\mathrm{\Delta}z/\mathrm{\Delta}y{)}^{2}$, which is a monotonic decreasing function of $\theta $. In particular, the case $\theta =\pi /6$, i.e., when the relative angles between the interfering plane waves all equal to $2\pi /3$, the degree of anistropy equals to unity and the fluid deformation admits a hexagonal symmetry. Values of $\theta $, which differ from $\pi /6$, result in a fluid deformation with a distorted unit cell and of broken hexagonal symmetry. Interestingly, the unit cell area $\mathrm{\Delta}y\xb7\mathrm{\Delta}z/2$ is not preserved under changes of the interfering plane-waves’ relative angle, and the minimal unit cell area corresponds to $\theta =\pi /6$. Figures 2(b)–2(d) show a fluid deformation described by Eq. (7) for the values $\theta =\pi /6$, $\theta =\pi /4$, and $\theta =\pi /3$, respectively. Other values of the parameter $q$, which are close to $q=\stackrel{\u203e}{q}$, can be used to control fluid distribution within each unit cell. Values of $q$, which differ significantly from $q=\stackrel{\u203e}{q}$ can be employed to trigger phase transition of the fluid lattice. In particular, increasingly smaller values of $q$ (i.e., larger ${a}_{\mathrm{1,2}}/{a}_{3}$) and fixed propagation directions, which do not modify $\mathrm{\Delta}y$ and $\mathrm{\Delta}z$ lead to change of the fluid deformation symmetry from hexagonal, which is composed of two equivalent interpenetrating triangular Bravais lattices that admit two primitive vectors ${\overrightarrow{V}}_{\mathrm{1,2}}$ and and two cites per unit cell,^{51} to a face centered cubic symmetry with one primitive vector per unit cell. Figures 2(d), 2(h), and 2(i) show TLD film deformation for the case $\theta =\pi /6$ and successively smaller values of $q$ leading to merging of fluid peaks (blue) within each cell and a phase transition from hexagonal lattice to a lattice with centered cubic symmetry.

A liquid lattice with a rectangular tunable symmetry can be formed by interfering four beams where the relative amplitudes satisfy ${a}_{\mathrm{0,2}}=qa$, ${a}_{\mathrm{1,3}}=a$ and the propagation angles are given by ${\theta}_{n}=(1+2n)\pi /4+(-1{)}^{n}\alpha $ ($n=\mathrm{0,1},\mathrm{2,3}$), where $\alpha $ is a real parameter. Inserting the corresponding intensity (see the Supplemental Material) into Eq. (2) yields the following TC-driven deformation of the TLD film:

## Eq. (10)

$$\begin{array}{c}\frac{\eta ({\overrightarrow{r}}_{\parallel},t)}{{h}_{0}}=-\frac{{\eta}_{\mathrm{max}}}{3}\left[\mathrm{cos}\right(\frac{2\pi y}{\mathrm{\Delta}y})+\mathrm{cos}(\frac{2\pi z}{\mathrm{\Delta}z})+\mathrm{cos}(\frac{2\pi y}{\mathrm{\Delta}y}\left)\mathrm{cos}\right(\frac{2\pi z}{\mathrm{\Delta}z}\left)\right].\end{array}$$## Eq. (11)

$$\mathrm{\Delta}y=\frac{\lambda}{2[\mathrm{cos}(\alpha )-\mathrm{sin}(\alpha )]};\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{\Delta}z=\frac{\lambda}{2[\mathrm{cos}(\alpha )+\mathrm{sin}(\alpha )]};$$## 1.3.

### Self-Induced Distributed Feedback Lasing in 1-D Optical Liquid Lattices

The self-induced 1-D deformation of the TLD film given by Eq. (6) [see also Fig. 2(a)] leads to a periodic modulation of the depth-averaged refractive index for the corresponding surface optical modes, which triggers a self-induced selection mechanism. The latter singles out a set of potential lasing modes that satisfy the minimum phase or the $q$’th-order Bragg resonant condition, given as

where $\mathrm{\Delta}$ is the mismatch parameter^{52}

^{,}

^{53}between the wave numbers ${k}_{\pm}^{L}$ of the lasing modes and $K$ is defined by Eq. (6). The case where the modes admit equal wave vectors and satisfy ${k}_{+}^{L}=-{k}_{-}^{L}\equiv {k}_{L}$ leads to a lasing mode wave vector of magnitude ${k}_{L}=qk\text{\hspace{0.17em}}\mathrm{cos}(\theta )$ and wavelength ${\lambda}_{L}=\lambda /[q\text{\hspace{0.17em}}\mathrm{cos}(\theta )]$, where the integer $q$ is determined by the spectral region of maximal spontaneous emission of the gain media and $\theta $ is the tunable angle of the interfering surface opitcal modes schematically described in Fig. 2(a). For instance, in a case of two counterpropagating ($\theta =0$) interfering TE WG modes of wavelength $1.55\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$, the resonant condition given by Eq. (12) leads to ${\lambda}_{L}=775\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$ and $q=2$. For the case $\theta =\pi /3$, identical resonant wavelength ${\lambda}_{L}=775\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$ is realized for $q=4$. Note that the effective mode index, ${n}_{\mathrm{eff}}$, defined via the relation $k={n}_{\mathrm{eff}}{k}_{0}$ cancels from both sides of the resonant condition, because in the self-induced case, the pumping beam coincides with the beam that forms the grating. In the limit case of nonself-induced lasing regime, i.e., in a case the pair of interfering surface optical plane waves with wave vectors of equal magnitude $k$ is used to form and tune the 1-D liquid lattice of periodicity $\mathrm{\Lambda}=\pi /[k\text{\hspace{0.17em}}\mathrm{cos}(\theta )]$, whereas the pumping is performed by an additional unfocused incident source of free-space wavelength ${\lambda}_{p}$ (e.g., similarly to the method employed in Ref. 29), the lasing condition takes the form $2{n}_{\mathrm{eff}}\mathrm{\Lambda}-q{\lambda}_{p}=0$. For instance, a 300-nm thick suspended liquid film of index 1.409, the corresponding ${n}_{\mathrm{eff}}$ of the TE mode is given by ${n}_{\mathrm{eff}}=1.319$, and increasing $\theta $ from 0 to $\pi /6$, modifies the lattice period $\mathrm{\Lambda}$ from 1 to $2\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ and consequently shifts the resonance condition from $q=5$ and ${\lambda}_{p}=527.6\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$ to $q=7$ and ${\lambda}_{p}=753.7\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{nm}$.

Note that, similar to other self-induced lasing mechanisms,^{54} the field intensity plays a double role; it determines both the coupling constant^{55}^{,}^{56} of the backward Bragg scattering (which provides the feedback mechanism) as well as the pump intensity necessary for lasing. Nevertheless, in our case, the amplitude of the optically induced periodic liquid lattice and consequently the coupling between the right- and the left-propagating modes is a more pronounced function of the optical intensity due to the higher difference between the typical refractive indices of liquids and gases. Notably, the relatively low optical power needed to introduce index changes induced by the TC/SC gas–liquid interface deformation [see Ref. 36, below Eq. (13) for the TC case] is expected to lead to substantially lower lasing threshold as compared to a thick liquid film or a liquid without gas–liquid interface (e.g., liquid fully occupying a microchannel), where the liquid–light interaction due to interface deformation is expected to be absent. In particular, in the absence of gas–liquid interface, index changes in liquids mostly stem from changes of polarizability, population of the electronic states (see Ref. 50 and references within), and material density; for instance, the index change due to thermo-optical effect, $\mathrm{\Delta}{n}_{\mathrm{TO}}$, is given by $\mathrm{\Delta}{n}_{\mathrm{TO}}={\alpha}_{\mathrm{TO}}\mathrm{\Delta}T$, where $\mathrm{\Delta}T$ is a temperature change and ${\alpha}_{\mathrm{TO}}$ is the thermo-optical coefficient, which is typically of the following order of magnitude ${\alpha}_{\mathrm{TO}}={10}^{-4}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{K}}^{-1}$. In the case of TC effect, the corresponding index change driven by identical temperature difference is given by $\mathrm{\Delta}{n}_{\mathrm{TC}}={\alpha}_{\mathrm{TC}}\mathrm{\Delta}T$, where ${\alpha}_{\mathrm{TC}}=\frac{3{\sigma}_{T}{d}_{z}^{2}}{2{\sigma}_{0}{h}_{0}^{2}{\pi}^{2}}\frac{{(n-1/2)}^{2}}{{n}^{4}}$ and $n$ is the integer indicating the numbers of optical periods that fit into the slot (see the Supplemental Material for derivation). Inserting the typical numbers used below Eq. (6) and taking $n=20$, we learn that ${\alpha}_{\mathrm{s}}=187.6$. The latter indicates that under similar temperature change, the amplitude of the periodically induced index and consequently also the coupling coefficient between the counter propagating modes (which directly determines the lasing threshold) is significantly higher for the TC case.

In the limit of small deformation and small intensity, the corresponding coupling constant can be written as $\kappa ={\kappa}_{0}{I}_{0}$, where ${\kappa}_{0}$ is the constant of proportionality given as

## Eq. (13)

$${\kappa}_{0}=k\overline{\mathrm{M}}{h}_{0};\phantom{\rule[-0.0ex]{2em}{0.0ex}}\overline{\mathrm{M}}=\frac{2\mathrm{Ma}}{{\tau}_{th}{\lambda}_{N}}\xb7\frac{{\alpha}_{\mathrm{th}}{d}^{2}}{{k}_{\mathrm{th}}\mathrm{\Delta}T}.$$^{54}

## Eq. (14)

$$\kappa L\text{\hspace{0.17em}}\mathrm{sinh}(i\gamma L)=\pm \gamma L;\phantom{\rule[-0.0ex]{1em}{0.0ex}}\gamma =-i{[{(\frac{g}{2}+i\mathrm{\Delta})}^{2}+{\kappa}^{2}]}^{1/2},$$## 1.4.

### Bandstructure of WG and SPP Modes in Hexagonal and Rectangular Liquid Lattices

Consider the case of interferring plane waves of the type presented in Eq. (5), which interact via nonlocal and nonlinear integral term in Eq. (4). Representing the intensity distribution given by Eq. (5) as a linear superposition of elements in the ${\phi}_{m,n}$ basis and employing orthogonality of these basis functions leads to the following dispersion relation:

## Eq. (15)

$$[2i{\beta}_{0}{k}_{z}^{(i)}+{({k}_{y}^{(i)})}^{2}]{a}_{i}=4{\chi}_{\mathrm{TC}/\mathrm{SC}}\sum _{m,n=0}^{N-1}{\mathrm{\Delta}}_{imn}\frac{{a}_{i}{a}_{m}{a}_{n}}{{\lambda}_{\overline{m},\overline{n}}};\phantom{\rule[-0.0ex]{2em}{0.0ex}}{\mathrm{\Delta}}_{imn}\equiv \frac{1}{{d}_{y}{d}_{z}}\underset{0}{\overset{{d}_{y}}{\int}}\underset{0}{\overset{{d}_{z}}{\int}}{e}^{-i{\overrightarrow{k}}_{\parallel}^{(i)}\xb7{\overrightarrow{r}}_{\parallel}}{\phi}_{m,n}({\overrightarrow{r}}_{\parallel})\mathrm{d}{\overrightarrow{r}}_{\parallel},$$## Eq. (16)

$$2i{\beta}_{0}\frac{\partial A}{\partial z}+\frac{{\partial}^{2}A}{\partial {y}^{2}}-{V}^{(0)}({\overrightarrow{r}}_{\parallel})A=0;\phantom{\rule{0ex}{0ex}}{V}^{(0)}({\overrightarrow{r}}_{\parallel})\equiv -\mathrm{\Delta}{\epsilon}_{d}=-{\chi}_{\mathrm{TC}/\mathrm{SC}}\int \mathrm{d}{\overrightarrow{r}}_{\parallel}^{\prime}{G}_{l}({\overrightarrow{r}}_{\parallel},{\overrightarrow{r}}_{\parallel}^{\prime}){|{A}^{(0)}({\overrightarrow{r}}_{\parallel}^{\prime})|}^{2}.$$## Eq. (17)

$${V}^{(0)}({\overrightarrow{r}}_{\parallel})={V}_{0}[\mathrm{cos}(2\overline{n}\beta x)+2\text{\hspace{0.17em}}\mathrm{cos}(\overline{n}\beta x)\mathrm{cos}(\overline{m}\beta y)];\phantom{\rule[-0.0ex]{0em}{0.0ex}}{V}_{0}\equiv \overline{q}b\mathrm{M}/({\lambda}_{\overline{n},0}{\tau}_{\mathrm{th}})=b{\eta}_{\mathrm{max}}/3,$$To get an insight into light propagation in a hexagonal crystal, one usually adopts the tight binding approximation,^{57} which describes the dynamics of light as hopping between nearest lattice sites. This approximation is applicable to cases when the potential well at each site is sufficiently deep, leading to localized intensity of the Bloch modes around the sites. For instance, one could in principle expand the potential Eq. (17) around the peak up to second-order,^{57} and then approximate it by the exponential function, to obtain closed-form expressions for the orbitals as a function of the constant ${V}_{0}$ (see the Supplemental Material). Since the 2-D hexagonal fluidic lattice admits two sites per one unit cell, the corresponding Wannier function is a linear combination of the two associated orbitals. By taking the continuous limit of the discrete system,^{57} the Schrödinger equation around high symmetry points in the k-space can be represented as a pair of Dirac equations. ^{58} The latter implies the presence of a Dirac point around the corresponding high symmetry points (which usually appears around $K$ points), though it may also emerge in the inner regions of the Brillouin zone.

To determine the dispersion relation for SPP or slab WG modes in the regime that is not satisfied by the tight binding approximation and captures the effect of broken hexagonal and cubic symmetries, we turn to numerical simulations. Figure 3 presents simulation results of the bandstructure for three interfering beams, with angles $\theta =\pi /6$, $\pi /4$, and $\pi /3$, which led to a formation of the suspended TLD film configuration [shown in Fig. 1(b)] with fluid topography shown in Figs. 2(a)–2(c), respectively. Interestingly, optical liquid lattice with increasingly higher broken hexagonal symmetry, due to modified interference pattern, shows emergence of Dirac points for the case $\theta =\pi /4$ and $\theta =\pi /3$, i.e., regions with linear dispersion relation, which do not appear in the case $\theta =\pi /6$. The latter shows formation of two Dirac points [marked with black arrows in Figs. 3(c) and 3(f)] and admits the values 352.3 and 354.3 THz for TE mode and 350.9 and 347.4 THz for the TM mode. Notably, while Dirac cones are known to emerge in all-dielectric, solid-made, and nonreconfigurable setups, in our setup they emerge in reconfigurable optical liquid lattices. In particular, taking advantage of the compliant property of the liquid film allows position tuning of the Dirac points under optically induced surface tension stresses, which is markedly different than other recently suggested mechanisms to control topological properties of light in liquids, such as electrical tunability of liquid crystals.^{59}

## 1.5.

### Self-Induced Distributed Feedback Lasing and Tuning in Two-Dimensional Liquid Lattices

Figure 4 shows simulation results of the lasing frequencies and the corresponding wave vectors in the k-space, within plasmonic and suspended photonic liquid lattices of tunable symmetry formed by interfering surface optical waves described by angles $\theta $ and $\alpha $, described in Fig. 2. Figure 4(a) presents lasing frequencies of slab WG TE and TM modes (blue and red circles, respectively) in suspended photonic liquid lattices for several values of the interference angle $\theta $, as well as lasing frequencies of SPP modes propagating within plasmonic liquid lattices formed by four interfering SPPs with a relative interference angle set by $\alpha $ (red squares). Interestingly, unlike the 1-D case, the lasing frequency is not a monotonic function of the relative angle between the interfering waves. For simplicity, in our simulations, the laser material is set as a linear gain Lorentz model,^{60} described by a dielectric function given by $\u03f5(f)={\epsilon}_{l}+{\epsilon}_{L}{\omega}_{0}^{2}/[{\omega}_{0}^{2}-4\pi i{\delta}_{0}f-{(2\pi f)}^{2}]$ with resonance frequency centered at ${\omega}_{0}/(2\pi )=374.97\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{THz}$ (800 nm), liquid permittivity ${\epsilon}_{l}={n}_{l}^{2}={(1.409)}^{2}=1.9852$, Lorentz linewidth $\delta /(2\pi )=52.521\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{THz}$, and Lorentz permittivity ${\epsilon}_{L}=-0.0075$; the metal is taken as gold.

## 2.

## Discussion

In this work, we theoretically investigated optical properties of configurable liquid-made lattices, formed by interfering SPPs or slab WG modes. We leveraged the underlying complex nonlocal and nonlinear interaction described by the GL and Schrödinger equations, which capture the effect of the self-induced action of the optical mode on itself due to geometrical changes of the gas–liquid interface, to predict formation of optical liquid lattices and bandstructure tuning due to symmetry breaking of hexagonal symmetric and square symmetric lattices as well as phase transition effects between hexagonal symmetric to face centered symmetric lattices. We then applied the bandstucture tuning to demonstrate control over various properties of the lasing systems, such as gain threshold, lasing frequency, and emission direction of the corresponding lasing mode. Notably, the self-induced lasing threshold of the TLD film interacting with the surface optical mode is expected to admit much lower values as compared to a similar liquid system without gas–fluid interface and therefore has the prospect to serve as a future Lab-on-a-Film bio-sensing platform, which integrates liquid delivery with self-induced DFB lasing mechanism. Interestingly, the formation of a graphene-like liquid lattice with tunable Dirac points in lattices with hexagonal broken symmetry is substantially different from other configurable platforms introduced to date for the formation of Dirac points in optical lattices, such as schemes that incorporate cold atoms^{61} and Fermi gases.^{62} In addition, since metamaterials can be utilized as an optical computational platform, as was recently demonstrated in a solid-made and non-reconfigurable setup for solution of linear equations,^{63} the adaptive property of optical liquid lattices has the potential to allow reconfigurable computation of computationally challenging problems of systems of linear and nonlinear algebraic equations [e.g., Eq. (15)] and also has the potential to serve as an emulator of many-body quantum mechanical problems, such as electron propagation in an atomic lattice, including topological edge-state effects. We hope that our work will stimulate future experimental and theoretical studies to realize optical liquid lattices and to explore underlying nonlinear light–liquid interaction mechanisms that include also birefringence and magneto-optical effects.

## Acknowledgments

SR cordially thanks Dr. Anna Rubin-Brusilovski for valuable and inspiring discussions. This work was supported by the Defense Advanced Research Projects Agency (DARPA) DSO’s NAC (HR00112090009) and NLM, the U.S. Office of Naval Research (ONR) Multidisciplinary University Research Initiative (MURI), the U.S. National Science Foundation (NSF) Grant Nos. CCF-1640227, the Semiconductor Research Corporation (SRC), and the Cymer Corporation. The authors have no relevant financial interests in the manuscript and no other potential conflicts of interest to disclose.

## References

## Biography

**Shimon Rubin** received his BSc and MSc degrees from Technion, Israel in 2001 and 2004, respectively and his PhD on theoretical physics from Ben-Gurion University of the Negev in 2012. Between 2013 and 2016 he served as postdoctoral fellow at Technion at the Department of Mechanical Engineering, and since 2017 he has served as a postdoctoral researcher and assistant project scientist at the Department of Electrical and Computer Engineering (ECE) at University of California San Diego (UCSD).

**Yeshaiahu Fainman** received his MSc and PhD degrees from Technion, Israel in 1979 and 1983, respectively. He joined the ECE department at UCSD in July 1990, following a faculty appointment at the University of Michigan, and now serves as a Cymer chair and distinguished professor, directing the Ultrafast and Nanoscale Optics Group in research of nanophotonics. He has contributed over 280 manuscripts and 450 conference proceedings and is a fellow of OSA, IEEE, and SPIE, and a recipient of numerous prizes and awards.