Nonlinear, Tunable and Active Optical Metasurface with Liquid Film

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 2D metasurfaces and the underlying interaction regimes between surface optical modes and liquids. In this work, 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 lead to formation of 2D optical liquid-made lattices/metasurfaces with tunable symmetry and which can be leveraged for tuning of lasing modes. Furthermore, we show that the symmetry breaking of the 2D optical liquid lattice leads to phase transition and tuning of its topological properties which allows to form, destruct and move 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.

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 selfinduced 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.

Light-Induced Interaction between Surface Optical
Modes and a TLD Film The set of coupled governing equations that describes lightfluid 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 masstransport equations, relevant for the TC effect and SC effect, respectively, which constitute the coupling mechanism between light propagation and dynamics of TLD film. Specifically, lightinduced changes of the surface tension lead to deformation of the TLD film, ηðr ∥ ; tÞ, which in turn is coupled back to the propagation of light due to the associated spatial changes of the dielectric function, Δϵ D , which in the leading order of η∕h 0 is given as 36 Here, h 0 is the initial thickness of a flat and undeformed TLD film, which we assume to be thinner than the penetration depth of the corresponding optical surface mode into the bulk, and b is the mode-dependent coefficient we discuss below. First, we describe the case of a TC-driven interaction. Applying quasistatic temperature field distribution for the thermal transport and the thin film deformation linear equations yields the following nonlocal relation between the deformation η and the optical intensity I ≡ jEj 2 : where M ≡ Ma · χ∕2; and G l is the corresponding Green's function of the thin film equation (see Supplemental Material of Ref. 36). Here, Ma ¼ σ T ΔTh 0 ∕ðμD th Þ is the dimensionless Marangoni number that represents the ratio between the surface tension stresses due to the TC effect and dissipative forces due to fluid viscosity and thermal diffusivity; χ ≡ α th d 2 I 0 ∕ðk th ΔTÞ is the dimensionless intensity of the heat source; I 0 is the characteristic optical intensity; D th ¼ k th ∕ðρc p Þ is the heat diffusion coefficient; ρ, c p , k th , and α th are the mass density, specific heat, heat conductance, and optical absorption coefficient, respectively; τ th ¼ d 2 ∕D th is the typical time scale; d is the typical length scale along the in-plane direction. Importantly, Eqs. (1) and (2) are valid for both SPP and slab WG modes. In particular, the optical absorption coefficient α th , which drives thermal interaction between SPPs and TLD film or WG and TLD film could stem from ohmic losses in metals or from optical absorption in dielectrics, respectively. In fact, optical absorption in dielectrics could be enhanced by mixing the liquid with strong nanoabsorbers, as was demonstrated by doping liquid with concentration of 2-to 4-nm-diameter CdSe or CdTe nanoparticles of concentration 10 22 m −3 (i.e., ∼0.1% of the total liquid volume), triggering TC effect in capillaries under 514-nm illumination, 45 and optical power of ∼120 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 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 σ, 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 η 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): Here, M c ≡ Ma c · c 0 ∕Δc whereas Ma c ≡ h 0 σ c Δc∕ðμDÞ is the dimensionless SC Marangoni number, which represents the ratio between the surface tension stresses due to the SC effect and dissipative forces due to fluid viscosity and molecular diffusivity; τ c ¼ d 2 ∕D is the typical diffusion time; and τ on ¼ 1∕k on is the typical transformation rate for a given optical intensity (see the Supplemental Material). 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): A complex GL equation with a local potential was shown to capture the effects of nonlinear Kerr media (see Ref. 47 and references within), propagating SPPs adjacent to solid dielectric with gain, 48 as well as nanofocusing of SPPs in tapered plasmonic waveguides. 49 The nonlocal integral term, V½Aðr ∥ Þ, 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 χ TC∕SC ¼ fbM 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 χ TC∕SC ¼ bM. 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Γ, 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 where δε 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).

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 θ n , which are measured relative to the positive direction of the y axis. The corresponding wave vectors are given bỹ k n ¼ k 0 ½cosðθ n Þ; sinðθ n Þ, where k 0 ¼ 2π∕λ and λ 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 λ ¼ λ 0 ∕n, where λ 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 ik n ·r ∥ is then given as In this work, we focus on the cases N ¼ 2,3; 4, which, as shown shortly below, correspond to optical liquid lattices with 1-D translation, 2-D hexagonal, and 2-D rectangular symmetries, respectively. Furthermore, we allow nonequal relative propagation angles between adjacent wave vectors and nonequal amplitudes, leading to broken symmetries and to changes of the symmetry group of the formed optical liquid lattices. 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 vectorsk AE ¼ k½AE cosðθÞ; sinðθÞ [see Fig. 2(a)] is given by 2I 0 f1 þ cos½KðθÞxg, which upon inserting into Eq. (2) yields the corresponding 1-D deformation of the TLD film: where λn is a constant (see the Supplemental Material). Importantly, the amplitude of the deformation is proportional to the Marangoni constant and the periodicity ΛðθÞ can be tuned by the angle θ. Employing the definitions for M; χ, and τ th near Eq. (2), the modulation amplitude of the cosine deformation described in Eq. (6) can be written as 2M∕ðτ th λnÞ ¼ ½3σ T α th d 4 x ∕ðσ 0 k th h 2 0 π 4 ÞI 0 . Inserting the following typical parameters that can be satisfied by silicone oil on gold substrate: and n ¼ 20; we learn that 2M∕ðτ th λnÞ ¼ 1.5 · 10 −7 I 0 . Consequently, using Eq. (6), the optical intensity required to generate periodic deformation η of amplitude 450 nm, which is similar to TLD film thickness (i.e., η∕h 0 ≈ 1), is of the order of magnitude I 0 ¼ 6.5 · 10 6 Wm −2 .
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 ≡ a, where a and q are real constants and assume the corresponding propagation directions are given byn 1 ¼ ½cosðθÞ; − sinðθÞ,n 2 ¼ −½cosðθÞ; sinðθÞ, n 3 ¼ ð0,1Þ (which can be also written as θ n ¼ ðn þ 1Þπþ ð−1Þ n θ for n ¼ 0,1 and θ 2 ¼ π∕2). Inserting the corresponding intensity distribution (see the Supplemental Material) into Eq. (2) yields the following TLD film deformation: Here, Δy and Δz are the corresponding periodicities along the y and z directions, respectively which are set by the angle θ and the wavelength of the optical surface wave λ but do not depend on the amplitude a nor on q, and η max is a constant given as η max ≡ ð3qMÞ∕ðλn ;0 τ th Þ;q ¼ 2λ 2n;0 ∕λm ;n ; where m; n are integers that satisfy m ¼ 2d y ∕Δy and n ¼ 2d z ∕Δz (see the Supplemental Material for more details). Expanding the deformation given by Eq. (7) up to a secondorder in Taylor series around one of the maximum points [e.g., ðΔy∕3,0Þ] yields a harmonic oscillator potential that admits a degree of anisotropy given by 3ðΔz∕ΔyÞ 2 , which is a monotonic decreasing function of θ. In particular, the case θ ¼ π∕6, i.e., when the relative angles between the interfering plane waves all equal to 2π∕3, the degree of anistropy equals to unity and the fluid deformation admits a hexagonal symmetry. Values of θ, which differ from π∕6, result in a fluid deformation with a distorted unit cell and of broken hexagonal symmetry. Interestingly, the unit cell area Δy · Δz∕2 is not preserved under changes of the interfering plane-waves' relative angle, and the minimal unit cell area corresponds to θ ¼ π∕6. Figures 2(b)-2(d) show a fluid deformation described by Eq. (7) for the values θ ¼ π∕6, θ ¼ π∕4, and θ ¼ π∕3, respectively. Other values of the parameter q, which are close to q ¼ q, can be used to control fluid distribution within each unit cell. Values of q, which differ significantly from q ¼ q can be employed to trigger phase transition of the fluid lattice. In particular, increasingly smaller values of q (i.e., larger a 1,2 ∕a 3 ) and fixed propagation directions, which do not modify Δy and Δ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Ṽ 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 θ ¼ π∕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 0,2 ¼ qa, a 1,3 ¼ a and the propagation angles are given by θ n ¼ ð1 þ 2nÞπ∕4 þ ð−1Þ n α (n ¼ 0,1; 2,3), where α is a real parameter. Inserting the corresponding intensity (see the Supplemental Material) into Eq. (2) yields the following TCdriven deformation of the TLD film: Here, Δy and Δz are the corresponding periodicities along the y; z directions given as Δy ¼ λ 2½cosðαÞ − sinðαÞ ; Δz ¼ λ 2½cosðαÞ þ sinðαÞ ; (11) and η max is a constant given by η max ≡ 3Mðq − 3∕4Þ∕ð2λn ;0 τ th Þ, where m; n are integers that satisfy m ¼ 2d y ∕Δy and n¼2d z ∕Δz, and the constant q satisfies the quadratic equation ð1∕2−q 2 ∕4Þ∕ ð2λm ;n Þ¼ðq−3∕4Þ∕λn ;0 . Since λ m;n and λ n;0 are both positive, the latter quadratic equation is guaranteed to admit two different and positive solutions for q. Figures 2(d)-2(f) show a TLD film deformation described by Eq. (10) for increasingly higher values of α (α ¼ 0 deg, 10 deg, 20 deg), leading to an increase of Δy and to a decrease of Δz as follows from Eq. (11). Similarly, the case of SC-driven deformation of TLD film is obtained by substituting the corresponding intensity into Eq. (3), leading to identical expressions for the deformations given by Eqs. (7) and (10), where the ratio M∕τ th in η max is replaced by M c c b ∕τ c Δc.

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 Δ is the mismatch parameter 52,53 between the wave numbers k L AE 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 − ≡ k L leads to a lasing mode wave vector of magnitude k L ¼ qk cosðθÞ and wavelength λ L ¼ λ∕½q cosðθÞ, where the integer q is determined by the spectral region of maximal spontaneous emission of the gain media and θ is the tunable angle of the interfering surface opitcal modes schematically described in Fig. 2(a). For instance, in a case of two counterpropagating (θ ¼ 0) interfering TE WG modes of wavelength 1.55 μm, the resonant condition given by Eq. (12) leads to λ L ¼ 775 nm and q ¼ 2. For the case θ ¼ π∕3, identical resonant wavelength λ L ¼ 775 nm is realized for q ¼ 4. Note that the effective mode index, n eff , defined via the relation k ¼ n 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 Λ ¼ π∕½k cosðθÞ, whereas the pumping is performed by an additional unfocused incident source of free-space wavelength λ p (e.g., similarly to the method employed in Ref. 29), the lasing condition takes the form 2n eff Λ − qλ p ¼ 0. For instance, a 300nm thick suspended liquid film of index 1.409, the corresponding n eff of the TE mode is given by n eff ¼ 1.319, and increasing θ from 0 to π∕6, modifies the lattice period Λ from 1 to 2 μm and consequently shifts the resonance condition from q ¼ 5 and λ p ¼ 527.6 nm to q ¼ 7 and λ p ¼ 753.7 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 gasliquid 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 gasliquid 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 gasliquid 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, Δn TO , is given by Δn TO ¼ α TO ΔT, where ΔT is a temperature change and α TO is the thermo-optical coefficient, which is typically of the following order of magnitude α TO ¼ 10 −4 K −1 . In the case of TC effect, the corresponding index change driven by identical temperature difference is given by 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 α 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 κ ¼ κ 0 I 0 , where κ 0 is the constant of proportionality given as Here, h 0 is the mean thickness of the TLD film and k is a modedependent coefficient; e.g., for WG TE-TE coupling, the expression for k is given in Ref. 56 (see the Supplemental Material). The coupling coefficient affects the oscillation condition and for DFB laser cavities it can be written as the following complex equation: 54 where the corresponding solution for g and Δ must satisfy simultaneously vanishing of the real and the imaginary parts of Eq. (14). While index modulation typically supports lasing at the edges of Brillouin zone, the case of a gain coupling, which is described by Eq. (14) with κ → iκ, results in lasing at the center of the Brillouin zone. In our study, both can be realized depending on the dielectric properties of the TLD film. Figure  S3 in the Supplemental Material shows the solution of Eq. (14) in the index and gain modulation regimes, as a function of increased intensity leading to decreasing of the lasing threshold.

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 φ m;n basis and employing orthogonality of these basis functions leads to the following dispersion relation: Δ imn a i a m a n λm ;n ; where the nonlocal term is reduced to a cubic nonlinear term and n; m satisfy k Here, V ð0Þ ðr ∥ Þ ¼ −χ TC∕SC a ð0Þ m a ð0Þ n φ m;n ðr ∥ Þ∕λ m;n is the self-induced potential function and a ð0Þ are the lowest order nonperturbed envelope functions. For the specific case of three interfering beams with propagation directions specified following Eq. (5) and amplitudes that satisfy the condition specified following Eq. (7), the TLD film deformation is given by Eq. (7) and the corresponding induced optical potential for TC case is given as V ð0Þ ðr ∥ Þ ¼ V 0 ½cosð2nβxÞ þ 2 cosðnβxÞ cosðmβyÞ; where the constant V 0 depends on the optical intensity and the Marangoni constant. In particular, for a positive Marangoni constant σ T > 0 (M > 0) and b > 0 the TLD film deformation given by Eq. (7) admits peak values of TLD film deformation in a triangular lattice and minimum values in a hexagonal lattice, and vice versa for the induced potential V 0 ðr ∥ Þ. The opposite case of M < 0, which could be realized either by σ T < 0 and b > 0, or by σ T > 0 and b < 0, yields an opposite behavior, i.e., liquid accumulation and potential minima around triangular lattice points, and liquid wells and potential maxima around hexagonal lattice points.
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 , and (f) π∕3. Two Dirac points emerge at the angle θ ¼ π∕3 at frequencies around 353 and 349 THz, respectively, marked at (c) and (f) by black arrows. In the simulation, the refractive index is 1.409, the mean TLD film thickness is 450 nm, and the peak-to-peak undulation amplitude is 300 nm. 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 θ ¼ π∕6, π∕4, and π∕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 θ ¼ π∕4 and θ ¼ π∕3, i.e., regions with linear dispersion relation, which do not appear in the case θ ¼ π∕6. 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

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 θ and α, 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 θ, 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 α (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 2 with resonance frequency centered at ω 0 ∕ð2πÞ ¼ 374.97 THz (800 nm), liquid permittivity ε l ¼ n 2 l ¼ ð1.409Þ 2 ¼ 1.9852, Lorentz linewidth δ∕ð2πÞ ¼ 52.521 THz, and Lorentz permittivity ε L ¼ −0.0075; the metal is taken as gold.

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 gasliquid 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 selfinduced 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 nonreconfigurable 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 edgestate 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.