Numerical modeling of optical modes in topological soft matter

Vector and vortex laser beams are desired in many applications and are usually created by manipulating the laser output or by inserting optical components in the laser cavity. Distinctly, inserting liquid crystals into the laser cavity allows for extensive control over the emitted light due to their high susceptibility to external fields and birefringent nature. In this work we demonstrate diverse optical modes for lasing as enabled and stablised by topological birefringent soft matter structures using numerical modelling. We show diverse structuring of light -- with different 3D intensity and polarization profiles -- as realised by topological soft matter structures in radial nematic droplet, in 2D nematic cavities of different geometry and including topological defects with different charges and winding numbers, in arbitrary varying birefringence fields with topological defects and in pixelated birefringent profiles. We use custom written FDFD code to calculate emergent electromagnetic eigenmodes. Control over lasing is of a particular interest aiming towards the creation of general intensity, polarization and topologically shaped laser beams.

Typically, laser cavities are designed to generate simple Gaussian beams and then the output beam is manipulated and transformed [13,14] to create vortex and vector beams. However, vector beams can be created already within the laser cavity by inserting various optical components, including birefringent components [15] and liquid crystals [16]. Vortex beams can be created directly within the laser cavity by breaking the orbital angular momentum (OAM) degeneracy of a standard laser cavity, which requires a very specific setup with a series of optical elements [17]. An ideal laser cavity would be geometrically simple and easy to manufacture, but optically capable of generating arbitrary-vector or vortex-beams, which moreover, would be tunable and switchable. Creating such cavities is an open challenge.
Due to their unique properties, such as birefringence, self-assembly characteristics and ability to easily tune their optical properties by applying external fields [18][19][20], liquid crystals are employed in a variety of photonic devices. Configuration of liquid crystals is characterized by a nematic director, a headless vector with n and −n fully equivalent, which points in the direction of average molecular orientation. The alignment of molecules is reflected in coupling of the material with electromagnetic fields as the direction of optical axis corresponds to the director. Liquid crystals can be used to generate vortex beams by passing light through a q-plate [21], a radial nematic liquid crystal droplet [22] or through spontaneously formed topological defects [23] and vector beams by passing a light beam along the disclination line [24].
Light-matter interaction can be bidirectional, with complex beams also influencing the liquid crystal material by optomechanical effects [25], generation of liquid crystal structures [26] and induction of liquid-crystal vortices which consequently create optical vortices [27]. Liquid crystal based lasers have been demonstrated in different liquid crystal phases [28] and different geometries [29], including Fabry-Pérot cavities [30,31]. Most studies on liquid crystal lasers however focus on the tunability of the lasing wavelength and not the shaping of the laser beam. Distinctly, in our recent work in [31], the focus was on how to achieve lasing of selected passive modes and tune them in selected liquid crystal structures.
When the director field is forced in a configuration that cannot be continuously transformed into a uniform state, topological defects appear. Rather recently, different structures of highly nonuniform nematic director were realised that include colloidal crystal [32], surface imposed defect structure [33], nematic drops with handles [34] and even knotted and linked nematic fields [35,36]. Moreover, if chiral nematic (cholesteric) liquid crystal is used in addition to point and linear defect structures in droplets [37,38],other non-trivial and highly twisted configurations including topological solitons [39,40] like torons [41], hopfions [42] and skyrmions [43] can be achieved in planar cells. Structures can be further stabilized or tailored by use of patterned surface anchoring, metasurfaces and two-photon polymerization [44][45][46]. In highly twisted nematic fields and especially in the vicinity of topological defects nematic director field is highly distorted, leading to rapid changes in optical axis direction. Moreover, since electromagnetic fields are polar and nematic director is a non-polar vector, their topological constraints are in principle different, therefore interaction between light and material becomes highly non-trivial. By inserting such topological structures into a lasing cavity, generation of complex modes is expected.
In this paper we numerically calculate diverse complex optical vector modes from generalised lasing cavity based on self-assembled liquid crystal topological structures in a surrounding Fabry-Pérot cavity. We propose a few cavity designs for lasing applications based on liquid crystal structures inside the cavity. We show that well known cylindrical vector modes in a birefringent medium decouple into modes with polarization governed by nematic director field and that different topologies can be patterned into electric field. Theoretical description of liquid crystal polarization alignment to polarization projection is provided. Finally, we give few ideas for creating complex vector modes. More generally, the proposed design based on soft laser cavities has no principal limitation for realisation of arbitrary lasing beam of any intentsity or polarization profile.

Methods
All the numerical results presented in this article were calculated by a custom written Finite Difference Frequency Domain (FDFD) code, capable of solving a system with arbitrary profile of materials' dielectric permittivity tensor in the form of an eigenproblem. The Master equation for magnetic field was formulated from Maxwell curl equations: where the eigenvalue represents the frequency of the optical mode. Electric field E in every point of the cavity was calculated from the (nodal) eigenvector H.
Simulations ran on Intel Xeon nodes with 190GB RAM in Matlab R2019a environment. A fixed number of eigenmodes around the desired wavelength was calculated in each run, where in the material was kept large enough to obey the standard rule of thumb determining the voxel size of /10 or smaller. Perfectly matched layers (PML) with the thickness larger than /2 were used to truncate the domain and simulate infinite boundary conditions in transverse directions. Perfect electric conductor (PEC) boundary conditions were used on the top and bottom boundaries to reproduce the effects of perfect mirrors. Applying these symmetries does not affect any results and conclusions, and is used for the only reason of reducing the computational complexity of the problem (especially the needed computer memory). We also performed selected calculations in full geometries (i.e. without applied symmetries) and identical solutions were obtained. Quality factors were calculated directly from the real and imaginary part of and served as a tool to extract localized modes and classify them. Note, that the presented methodological approach does not include any formulation of gain materials for lasing, thus returning passive resonant modes. Nevertheless, using same methodological approach, good agreement between -in principle-passive resonant modes and experimentally observed lasing modes is observed in such topological soft matter structures [31].
Dielectric profile of the cavity was calculated from the nematic order parameter tensor : where is the degree of order and are components of the nematic director. e (1) ⊥ n is the secondary director and e (2) = n × e (1) . The second term in Eq. (2) accounts for biaxiality , which quantifies fluctuations around the secondary director e (1) . The anisotropic dielectric tensor was calculated from as [47]: where¯ is the average dielectric permittivity and mol a = ( − ⊥ )/ is the molecular dielectric anisotropy for a degree of order . Nematic director field, which directly affects the dielectric tensor, was determined analytically for topological defects with different topological charges in 3D and winding numbers in 2D, following Refs. [24,48].
Nematic soft matter structures have the distinct capability to provide topologically non-trivial birefringent patterns, that affect and possibly control the polarization, intensity and even angular momentum of the light, either under transmission, in resonators or even in laser cavities [16,22,31,49]. In this work, the non-trivial topology of the nematic fields is imposed by the presence of topological defects, i.e. the +1 radial point defects in the radial nematic droplet, and differentwinding number disclination lines that span across the nematic layer. These nematic defects strongly affect the topology of the beam polarization profiles-the 2D topological winding number of the polarization in the real-space-but in given structures not the angular momentum of the beams (which is zero), as the cavity does not distinguish between the OAM modes with the opposite handedness due to their identical spatial intensity distributions, radii of curvature on the wavefront and Gouy phase shifts [17]. Note, however, that more generally, especially chiral nematic structures could potentially give rise also to vortex beams with non-zero angular momentum momentum [50,51].

Modes in isotropic and radial nematic droplet
First, we analyzed eigenmodes of a Fabry-Pérot (FP) resonator containing a radial nematic droplet and compared them with the modes of a FP resonator containing a droplet of an optically isotropic material. The droplet was suspended in an isotropic medium and enclosed by two reflective surfaces, simulated by PEC boundary, as shown in Figure 1(a,b). To achieve sufficient resolution, simulation box was split into eight octants and only one octant with the size of 55 × 55 × 50 pixels,which roughly corresponds to 5.5 × 5.5 × 5 , was simulated at a time (scheme shown in Fig. 1(a)). Three different types of octants with PEC boundary conditions in vertical direction and various combinations of even and odd boundary conditions in lateral directions, shown in Fig. 1(c), were used to cover all possible symmetries of the cross sections of the resulting modes. In each case, two of the lateral boundary surfaces included PML boundary conditions, representing the outer-infinite-boundary of the lasing cavity. Other two lateral boundary conditions were either both set to PMC/PEC or one to PMC (perfect magnetic conductor) and the other to PEC to enforce even or odd boundary condition for electric field, respectively. Electric field outside the simulation box was obtained by using symmetry transformations.
Nematic director field with the radial configuration and a +1 defect in the middle of the droplet was determined by analytical formula. Radius of the droplet was set to 40 pixels (roughly 4 ), meaning that the droplet filled 80 % of the FP layer thickness. Ordinary ( o ) and extraordinary ( e ) refractive indices of liquid crystal and the surrounding isotropic material ( out ) were set to o = 1.54, e = 1.71, and out = 1.47, respectively, generally corresponding to the indices of widely used 5CB liquid crystal [52,53] and glycerol. In the case of the isotropic droplet, its refractive index was set to approximate index of the molten 5CB liquid crystal (i.e. 5CB in the isotropic phase), namely i = 2 3 o + 1 3 e ∼ 1.6. 3D electric field intensity distributions of selected modes are shown in Figure 1(d) and Figure  2(a). Optical eigenmodes of a FP resonator containing optically isotropic spherical droplet or particle are cylindrical vector beams [54] and are shown in Figure 2(b). Following the notation of vector modes in optical fibers, modes in panels (i) and (ii) are known as TE 01 and TM 01 modes, while modes in panels (iii) and (iv) are known as hybrid modes-HE 11 and HE 21 -and contain both, radial and azimuthal electric field contributions [55,56]. Mode in panel (v) is one of the higher order hybrid modes, where polarization axis performs multiple turns around the center of the beam.
Mode profiles in Figure 2(b) are shown in the plane close to the top reflecting surface of the cavity, where the refractive index is homogeneous in the lateral direction. Therefore the profile can not be identified as one of the modes that occur in the core region of the step-index optical fibers. In fact, it has been pointed out by Hall in his original article on vector beam solutions of Maxwell's equations [57] that the fiber mode would excite a vector Bessel-Gauss free space mode when reaching the end of the fiber-here approximated by a droplet in the center of the cavity-and propagating into the free space (i.e. homogeneous media). The profiles of the observed simulated vector beams, which are characterized by spatially varying polarization and uniform phase profile, are quantitatively compared to vector Bessel-Gauss (vBG) modes, also referred to as free space solutions in the following, by calculating the overlap integrals between two beam cross sections E sim and E an as [58]: where E sim is the cross section of numerically simulated mode in the intensity maximum closest to the top reflective surface and E an is the cross section of the analytical free space vector beam solution at the distance max from the waist of the beam, where max equals the distance between the center of the droplet and the plane of the E sim cross section.
Overlap integrals were calculated for selected vBG mode numbers and a range of parameters 0 ∈ (0.1 , 10 ) and 0 ∈ (0.1, 10), by varying both in steps of 0.1. Resulting overlap integrals for vBG mode number , for which the best agreement is achieved, are shown in Figure 2(c). For all five modes, the agreement is above 95 % (99 % for (i) and (iii), 98.5 % for (ii), 98 % for (iv) and 95 % for (v)), implicating that vBG beam profiles can in fact be used to describe the modes of the isotropic droplet.
If the birefringent droplet with radial profile of the optical axis/director is inserted in to FP cavity, two families of vector modes with space varying polarizaton profiles are observed: one with polarization purely parallel (radial polarization) and one with polarization purely perpendicular to the director field (azimuthal polarization) of the nematic droplet. Cross sections of mode profiles with the highest Q-factors from both families in the plane of the intensity maximum closest to the upper reflective surface are shown shown in Figure 2(d). In the presence of birefringence the optical modes of the isotropic droplet separate into purely radial and azimuthal contributions. TE 01 and TM 01 modes are still found as elementary radial and azimuthal solutions. All the hybrid modes decouple into solutions with multiple intensity maxima and polarization in either radial or azimuthal direction. We label the observed vector beams as VB radial and VB azimuthal modes, where is the mode number, which is reflected in the 2 intensity maxima in the azimuthal direction and is associated with the mode number of the beam emerging from the isotropic droplet (as determined by overlap integrals in Fig. 2(c)). Phase profiles of the modes are uniform over the cross sections, meaning that the modes do not possess phase singularities (i.e. optical vortices) and do not carry orbital angular momentum (OAM). Importantly, localization of both polarizations is ensured by refractive index of the surrounding material, which is lower than both, ordinary and extraordinary, indices of liquid crystal.
This discussion shows that the observed vector beams from a laser cavity with an embedded radial nematic droplet show similarities with vBG beams, but up to a varying degree of similarity. Similarly, vector Laguerre-Gaussian [59] solutions (i.e. not the scalar Laguerre-Gauss beams [4] with the helical phase profile) could be used, however neither of those sets can perfectly describe the emerging modes by only using a single vBG or vLG solution. Therefore, the modes in this work are labeled simply as vector beam modes with the mode numbers (VB ), meaning that they are in principle unique, as determined by the actual optically anisotropic cavity setup.

2D cavities
In order to further study the effects of cavity shape, director field configuration and different topology, we assumed that the Fabry-Perot cavity was filled entirely with liquid crystal. Inside a certain cylindrical or rectangular region the optical axis was directed parallel to the top and bottom mirrors (i.e. in-plane). Outside the region the director was gradually turned in the direction normal to the mirrors (see Fig. 3(a)-(c), outline of the cylindrical region is market with orange dashed outline). For simplicity, the director field was invariant in the direction normal to the mirrors, meaning the dielectric profile only varies in two dimensions, hence we name such profiles 2D cavities (although the third dimension is still needed to confine the light). The setup automatically ensured change in refractive index value and localized the modes with extraordinary polarization (i.e. polarization parallel to the director field), which was previously done by use of surrounding isotropic material, meaning that the cavity shape for the extraordinary polarization was determined by its corresponding refractive index. No symmetry boundary conditions were used. Computational domain had size of 80 × 80 × 30 pixels (≈ 8 × 8 × 3 ). Refractive indices were kept the same as in the case of the droplet.
Results for cylindrical and rectangular cavity with the homogeneous in-plane director field are shown in Figure 3(d) and (e). The shape of the intensity profile is clearly directed by the shape of the in-plane director field region, while the direction of the optical axis is directly imprinted in the polarization profile. Insets show that modes with the rotated intensity profile still keep the same direction of polarization.
In Figure 3(f) the in-plane director formed a defect line with a winding number of +1 spanning between top and bottom mirrors. Resulting modes are highly similar to the modes, that emerge from a radial nematic droplet, confirming that the polarization is mainly determined by the local direction of in-plane projection of the nematic director. It is noteworthy that the director manipulation as described here, will only localize modes polarized in the direction of the director field, as the perpendicular polarization sees no refractive index contrast, therefore no azimuthally polarized modes are found.
In the next step a defect line with winding number of +1/2 was studied, a structure which is allowed in a non-polar vector field of nematic director, but not in real vector electric field. Interestingly, resulting modes once again possess a polarization resembling the underlying director field (Fig. 4), but with odd number of intensity maxima in the azimuthal direction. Since the direction of the polarization of the emitted modes matches the optical axis direction of the liquid crystal, the liquid crystal acts effectively as a polarizer. Therefore we presuppose that the electric field polarization of the modes emerging around -invariant topological nematic defects E pol can be effectively described as a product of free-space solutions of vector Helmholtz equation E free (for example Hermite-Gaussian, vector Bessel-Gaussian or vector Laguerre-Gaussian beams) and the polarizer with the optical axis in the direction of the director field around the defect P as: We use analytical beam profiles at = 0, since the waist of the emerging beam is in this case expected to be centered at the end of the cylindrical region which acts as a waveguide. Calculated results for +1/2 defect are shown in Fig. 4, where the +1/2 polarizer field is written as P 1/2 ( , ) = cos(− 1 2 )e + sin(− 1 2 )e . With the appropriate selection of elementary solutions' parameters for vector Bessel-Gauss beams, good matching with simulations in terms of polarization direction can be achieved. Matching is once again quantified by overlap integrals. For mode numbers > 2 matching of up to 90 % is achieved. For smaller mode numbers, matching of about 30 %-70 % is found. The origin of lower matching likely lies in emergence of additional intensity maxima in the lateral direction (1 in the case of = 0 and = 2, 2 in the case of = 1), which cannot be described well by the vector beams, which are solutions of the Helmholtz equation in the free space, while we deal with the spatially varying refractive index. Maxima in the lateral direction emerge due to the fact that vector beams with low mode numbers have smaller radii (as seen from Fig. 2(b)) and do not fit into the cylindrical region for a wavelength, determined by the distance between mirrors. Modes with low mode numbers and single maximum in the lateral direction are expected to be found at a different ratio between the radius of the cylindrical region and the distance between mirrors or different wavelength. To generalize, local direction of spatially varying optical axis in -invariant birefringent cavity determines the polarization of the emerging mode and can be well described by using free space solutions and Equation (5). On the other hand the shape of its intensity profile is determined by the shape of the cavity in terms of refractive index and matching with free space solutions is lower.

Constructing complex modes
Above presented principles provide an interesting platform for custom design of modes with more complex intensity and polarization profiles. Figure 5(a)-(f) shows that bringing two cylindrical regions with ±1 defect lines in vicinity of each other makes their modes merge into a mode with non-trivial polarization winding. The distance between the defects at which the transition occurs is here mainly governed by the size of cylindrical region, which determines the size of the mode. When the ratio / < 2, regions begin to penetrate into each other and modes start to merge. However until / reaches value of 1 features of isolated modes are still visible (Fig. 5(b)). When / ≈ 1 the cylindrical regions are effectively merged into one, which is reflected in the shape of corresponding modes. Intensity profiles of isolated modes become less prominent. When the / < 1 the shape of the region is becoming more and more cylindrical which is reflected in the shape of the mode, which becomes less elongated and resembles the mode of the isolated defect. However the complex polarization patterns are governed by nematic director field with two topological defects (Fig. 5(c)). We expect the most interesting effects to occur when the distance between the defects is close to the wavelength of the light (here R/4, Fig. 5(d)). At this point it is hard to recognize modes of isolated defects, shape of the mode is again almost cylindrical, but with a polarization profile which is not observed in the modes of isolated defects. When the distance gets even smaller the wavelength of the light would be too large to be effected by the distortions of the director field. The defects would annihilate when the distance between their cores would get smaller than the correlation length in the LC. To showcase this effects we plot the intensity profile in the cross section at = /2 (marked by white dashed lines in Fig. 5(b)-(d)), where is the size of the system in direction (Fig. 5(e)) and integral of the intensity in the direction (Fig. 5(f)).
A distinct advantage of using liquid crystals as a medium inside cavity is their tunability and susceptibility to external fields. Figure 5(g)-(h) show that localized modes can be found in a continuous nematic director field. A sinusoidal modulation of angle between the mirror plane normal and director field with different periods in orthogonal lateral directions was used. As shown in the first panel, such configuration generates areas of in-plane director field with different profiles. Once again, modes in those regions occur and possess polarization which is parallel or perpendicular to the director field. Director configuration and consequently polarization of emerging modes could be further tailored by use of external fields, patterned surfaces or even by insertion of colloidal particles. Moreover, areas of in-and out-of-plane director field can be also achieved in a pixelated system where the direction is decided separately in each individual pixel with its size comparable to . As shown in Figure 5(i) intensity profile of the beam can be, in principle, shaped into custom patterns, corresponding to pixels with in-plane director. The wavelengths of all three modes shown differ by less than 3 %, size of a pixel is approximately 0.7 × 0.7 .

Discussion
In this work, the topology of the resonator-embedded liquid crystal structures combined with the actual geometry of the structures (droplet, layer,...) within the resonator is transferred distinctly into the topology of the light polarization. We extend the work of Ref. [31], which recently reported on experimental realization of lasing from liquid crystal structures inserted into Fabry-Pérot (FP) cavity, and here, explain the basic concepts how to design birefringent cavities and their belonging passive modes, distinctly focusing on the role of birefringence in the origin of petal-like modes, their relation to modes of isotropic spheres in FP cavities reported by Ref. [54] and free space vector beams. Additionally, we study 2D cavities as a simplified system, which allows us to use a wider variety of elementary defects in the nematic director with different winding numbers. Lastly we show--as a perspective and possible experimental challenge-that the arbitrary mode profiles can be created in different simplified (pixelated) design or as combination of previously presented cavity profiles.
The calculated passive modes can be used to anticipate the type and profiles of the active (lasing) modes which will emerge from a lasing cavity in the presence of gain media. However, the actual lasing mode profile -as above the lasing threshold-will importantly depend also on gain and pump distribution. Therefore, combinations of these parameters together with the geometry effectively determined within the passive modes will select the actual lasing modes, i.e. which will have the lowest lasing thresholds and what will be the mode profile above the lasing threshold.
Importantly, the spatial profiles of the light polarization presented here are shown to carry a nontrivial (i.e. nonzero) topological charge-near-field winding of the polarization in the real space at the top of the cavity, which would be an emitting facet of an active laser-(but not the orbital angular momentum-OAM) which is distinctly determined and coupled to the topological charge of the nematic topological defects embedded in the FP cavity (e.g. Fig. 4). Especially in the droplet resonators, such winding of the polarization is expected to be observed also in the far field as the calculated modes show similarities with vector Bessel Gauss beams, quantified by calculated overlap integrals, for which the diffraction and propagation solutions in free space are known to preserve the polarisation topology [60]. More generally, this shows a two way coupling between the topology of the LC structures and the topology of the emitted light (polarization). However in order to fully understand the propagation and far field of the general modes created though general liquid crystal topological structures, the actual lasing modes would need to be calculated first and then propagated to the far field regime.
The general approach of soft matter photonics has selected experimental chalenges which one could experience when realizing liquid crystal FP lasers, like photobleaching of dies, alignment of emission and absorption spectra of the materials with the modal spectrum of the cavity, width of the emission spectra of die, compared to spectral distance between the modes (mode selection), dependence of Q-factors and consequently stability of the belonging modes on the used materials' parameters and their birefringence etc. A range of these challenges were mitigated and/or solved in different setups, with examples of topological soft matter applications including for example cholesteric band edge lasers [28], blue phase lasers [61], high Q-factor WGM resonators [19] and others [62].
In broader terms the soft matter photonics has gained significant interest in recent years as it can offer alternatives to address some limitations of conventional solid state photonics [63]. Compared to their solid state counterparts, the soft photonic elements are well reconfigurable as they are able to adapt and respond to different changes in the surrounding environment [63]. In liquid crystal systems, the nematic birefringent profile which determines the optical properties can be controlled by applying external electric, magnetic or light fields, by confining the material into cavities and enforcing different surface (anchoring) conditions. Soft photonic elements also can be biocompatible [64,65] and can self assemble into 3D [66] and hierarchical structures [67], even with very fine feature sizes and high surface smoothness [19], which makes them useful for plasmonic elements [68], metamaterials [69] and laser cavities [70]. Some limitations of the soft matter approach to structuring light include photobleaching [71], possible material instability due to flows and evaporation, less control on smaller scales, and complex self-assembly mechanisms [72].

Conclusion
In conclusion, we have demonstrated that birefringent nematic structures inserted into Fabry-Pérot cavity highly impact the optical modes compared to the modes that would emerge from similar cavity based on isotropic components, especially, by decoupling the modes into families with predominant perpendicular or parallel polarization. Also, we show that the nematic director field configuration gets directly patterned in the modal electric field polarization, including through different topologies and cavity shapes, as determined by the (possibly polarization dependent) effective refractive index profiles. We have also shown that complex laser vector beam profiles can be constructed by use of space variant nematic fields, pixelated systems and topological defects, enabling the possibility to tune the light beam shape and polarization by external stimuli.
More generally, the work provides an insight into the control of the intensity and polarization of the passive light modes by using topological soft matter structures based on birefringent profiles of liquid crystals inserted into the Fabry-Pérot cavities. In addition to studies of active modes and lasing thresholds, an interesting extended approach would also be to use more complex-possibly chiral-liquid crystals structures, such as torons and cholesteric droplets to perform either as resonators or laser cavities, possibly also leading to structured light with topological helical light modality with net nonzero angular momentum. Soft photonic elements could also perform as building blocks of photonic chips [73,74] or light-by-light controllable logical devices [75,76].

Funding
This research was supported by the Slovenian Office of Science (ARRS) through grants J1-1697, P1-0099, N1-0195 and EU Horizon 2020 ERC advanced grant LOGOS.