Physics-based multistep beam propagation in inhomogeneous birefringent media

Guilhem Poy; Slobodan Žumer; Slobodan Žumer

doi:10.1364/OE.400984

1. Introduction

The beam propagation method (BPM) is one of the most wide-spread simulation toolbox for light propagation together with the Finite-Difference-Time-Domain (FDTD) method. These two methods allow the simulation of the propagation of optical fields in complex three-dimensional geometries and are thus highly relevant for photonics applications. While FDTD fully solves Maxwell equations and does not rely on any approximation other than numerical discretization, it generally requires a high number of mesh points over the distance of a wavelength in order to be accurate, which put constraints on the size of the computational mesh. On the other hand, BPM allows accurate simulation of light propagation on much coarser meshes provided that light don’t deviate too much from a given reference direction of propagation. BPM is very well established in isotropic media, with a number of sub-class of methods based on Finite-Difference [1], Finite-Element [2] or Fast-Fourier-Transform [3]. Among these methods, paraxial BPM makes use of the Slowly-Varying-Amplitude-Approximation to yield methods which are very efficient but limited in terms of accuracy, while wide-angle BPM is a class of methods which are more accurate (but a bit more expansive) [4]. In isotropic media, one of the major advantage of BPM in comparison to FDTD is that it is generally based on easy-to-invert sparse matrices representing differential operators on 2D meshes, thus leading to very efficient numerical schemes even in complex 3D systems.

In birefringent media such as liquid crystals, the anisotropic nature of the permittivity tensor complicates the derivation of paraxial and wide-angle BPM. Some approaches assume that some components of the permittivity tensor are zero [5] or the existence of translational invariance axis [6] to simplify the derivation, while other approaches do not make any assumptions on the permittivity tensor but use the Slowly-Varying-Amplitude-Approximation to get paraxial schemes [7]. In our opinion, the most general method that was proposed for light propagation in birefringent media is the wide-angle BPM of Vanbrabant et al. [8], which only assumes slow variation of the permittivity tensor along the reference propagation direction. The latter scheme is unfortunately based on very complicated nonsparse matrices, which limit the efficiency of the method and make its numerical implementation a delicate task.

In this paper, we present a unified theoretical framework for paraxial and wide-angle coherent beam propagation in inhomogeneous birefringent media. Our main contribution is to show that with a minimal set of physical assumptions, one can retrieve efficient BPM schemes based on sparse matrices representing 2D operators with an elegant physical interpretation, without having to make strong assumptions on the symmetry properties or smoothness of the permittivity tensor.

The plan of the paper is as follows. In Sec. 2, we derive the paraxial and wide-angle propagation equation for optical fields in anisotropic media and also present a method to take care of the continuity equation for electric and magnetic fields at interfaces of discontinuity of the permittivity. In Sec. 3, we validate our scheme by calculating the computational error on three simple birefringent systems, and then present a more complex example of application by demonstrating that topological solitons in frustrated cholesteric liquid crystal can be used as light waveguides for integrated liquid-crystal-based photonics circuit. Last, we give our concluding remarks and suggest possible extensions of our work with application to microscopy, photonic bandgap and nonlinear optics.

2. Beam propagation methods and continuity equations

2.1 General propagation equations and approximations

We consider the time-harmonic Maxwell electric field $\boldsymbol {E}(x,y,z)\exp \left (-ik_0 c t\right )$ (with $k_0$ the wavevector in empty space) propagating dominantly along the $z$ direction, and explicitly split the longitudinal degree of freedom ($E_z$) from the transverse degrees of freedom ($E_x$ and $E_y$). Using a compact notation based on dimensionless space coordinates $\{X=k_0x,Y=k_0y,Z=k_0z\}$ and Einstein’s convention on repeated Greek indices $\alpha ,\beta \in \{X,Y\}$ (i.e. Greek indices are only used for the transverse degrees of freedom), we find that the well-known wave equation for $\boldsymbol {E}$ [8] can be expressed under the following form:

(1)$$ \left[\left(\partial^{2}_Z+\Delta_\perp\right)\delta_{\alpha\beta}- \partial_\alpha\partial_\beta+\epsilon_{\alpha\beta}\right]E_\beta + \left[\epsilon_{\alpha z}-\partial_\alpha\partial_Z\right]E_z = 0, $$

(2)$$ \left[\epsilon_{z \beta}-\partial_Z\partial_\beta\right]E_\beta + \left[\Delta_\perp+\epsilon_{zz}\right]E_z = 0 , $$

where $\boldsymbol \epsilon$ is the permittivity tensor, $\delta _{\alpha \beta }$ is the Kronecker delta and $\Delta _\perp \equiv \partial _X^{2}+\partial _Y^{2}$.

Similarly to the general procedure in isotropic media [9], we want to eliminate $E_z$ from Eqs. (1) and (2) to obtain a simple propagation equation involving only the transverse field $E_\beta$. Before doing this, let us split the birefringent media in a series of slabs normal to the $z$-direction and assume that the permittivity tensor is $z$-independent in each slab – thereby greatly simplifying the derivation of the propagation equation in each slab. This approach is valid only if a sufficiently high number of slabs is chosen (an estimate of the computational error as a function of the slab thickness $\Delta _z$ will be given later on) and if the continuity equations for the electric field at the interface between each slabs are properly taken into account. In this section, we will only consider the propagation equations inside a slab, while the jump interface conditions are described later in a dedicated subsection. The geometry of the problem is schematically represented on Fig. 1.

Fig. 1. Schematic representation of the geometry of our computational scheme: the arbitrary birefringent system is split into a set of slabs normal to $z$ and we assume that in each slab the permittivity is $z$-independent; inside a slab, forward and backward propagating modes have their evolution controlled by two independent Schrödinger equations; continuity equations allow the calculation of the transfer of these modes through interfaces between different slabs.

Download Full Size | PDF

Eliminating $E_z$ from Eqs. (1) and (2), we therefore obtain after elementary algebra:

(3)$$ \boldsymbol{0} = \left[\mathbf{I}\partial_Z^{2}+\mathbf{S}\partial_Z+\mathbf{R}\right]\boldsymbol{E}_\perp, $$

(4)$$ E_z = \mathcal{H}^{-1}\left[\partial_Z\partial_\beta-\epsilon_{z\beta}\right]E_\beta, $$

where $\mathbf {I}$ is the identity matrix and $\boldsymbol {E}_\perp$ is the vector $\{E_x,E_y\}$. Equations (3) and (4) introduce the differential operators $\mathcal {H}$, $\mathbf {S}$ and $\mathbf {R}$, which only depends on the transverse coordinates $X$ and $Y$ and acts on arbitrary functions symbolized by a bullet $\bullet$ in the following. These operators are defined as:

(5)$$\mathcal{H}\,\bullet = \left[\epsilon_{zz}+\Delta_\perp\right]\bullet, $$

(6)$$S_{\alpha\beta}\,\bullet = \left[\delta_{\alpha\gamma}+ \partial_\alpha\frac{1}{\epsilon_{zz}}\partial_\gamma\right]\left[ \partial_\gamma\mathcal{H}^{-1}\epsilon_{z\beta} + \epsilon_{\gamma z}\mathcal{H}^{-1}\partial_\beta\right]\bullet, $$

(7)$$R_{\alpha\beta}\,\bullet = \left[\delta_{\alpha\gamma}+ \partial_\alpha\frac{1}{\epsilon_{zz}}\partial_\gamma\right]\left[ \delta_{\gamma\beta}\Delta_\perp-\partial_\gamma\partial_\beta+ \epsilon_{\gamma\beta}-\epsilon_{\gamma z}\mathcal{H}^{-1}\epsilon_{z\beta}\right] \bullet. $$

In all equations involving differential operators, we assume right-associativity ($D_1D_2\,\bullet \equiv D_1\left (D_2\left [\bullet \right ]\right)$ with $D_{1,2}$ two differential operators) and non-commutativity of operators involving the permittivity tensor (which can vary in space).

The reader can easily guess that the numerical discretization of the differential operators $\mathbf {S}$ and $\mathbf {R}$ cannot be done efficiently with simple sparse matrices, since their definitions involve multiple operator products and inverses. The authors of Ref. [8], who obtained a propagation equation similar to our Eq. (3) in the framework of a finite element BPM, made the same observation and deduced that the complex non-sparsity of the involved differential operators is the price to pay to get an accurate propagation method which takes into account the full anisotropy of the permittivity tensor $\boldsymbol \epsilon$. Here, we will show that an efficient and sparse BPM can be retrieved under a minimal set of approximations while retaining all effects related to the anisotropy of $\boldsymbol \epsilon$.

More specifically, we will assume in the remainder of this article the following physical limitation:

(H1) The permittivity contrast $\eta$ is small: $\exists \;\epsilon _\textrm {ref}$ s.t. $\eta ^{2}\ll 1$ with $\eta \equiv \left \lVert {\boldsymbol \epsilon -\epsilon _\textrm {ref}\mathbf {I}}\right \rVert _\infty /\epsilon _\textrm {ref}$.

where $\left \lVert {\bullet }\right \rVert _\infty$ is the uniform norm (maximum over all components and all space points of a tensor field) and $\epsilon _\textrm {ref}$ represents an average permittivity. Hypothesis (H1) is not very restrictive if light propagation is considered in pure birefringent media such as liquid crystals: even in liquid crystals of high birefringence $\eta$ is still small (e.g. we calculate $\eta <0.3$ in high birefringence isothiocyanato tolane liquid crystals [10], while in common liquid crystals used in display applications $\eta$ is typically of the order of $0.1$ [11]). In LC-based photonics devices confined with isotropic media, (H1) is still valid as long as the refractive indices of the confining materials have values comprised between (or similar to) the extraordinary and ordinary indices of the LC (which is true for typical thermotropic liquid crystals and silica glasses [12]).

Using (H1), the definitions of the differential operators $\mathbf {R}$ and $\mathbf {S}$ can be simplified to:

(8)$$\mathbf{S} = 2\left(\mathbf{W}+\mathbf{W}'\right)+\mathcal{O}\left({\eta^{2}}\right), \quad\quad\quad\quad \mathbf{R} = \tilde{\boldsymbol\epsilon}+\mathbf{L}+\mathcal{O}\left({\eta^{2}}\right),$$

where $\mathbf {W}'$, defined in Eq. (S1) of Supplement 1, is an unimportant differential operator which will be neglected later on with an additional physical assumption, and the new differential operators $\mathbf {W}$, $\mathbf {L}$ and $\tilde {\boldsymbol \epsilon }$ are defined as:

(9)$$ \mathbf{W} = \left[\mathbf{I}+\mathbf{W}^{(1)}\right]^{-1}\mathbf{W}^{(0)}, $$

(10)$$ W^{(0)}_{\alpha\beta} = \frac{1}{2}\left[ \partial_\alpha\frac{\epsilon_{z\beta}}{\epsilon_{zz}} + \frac{\epsilon_{\alpha z}}{\epsilon_{zz}}\partial_\beta\right], \quad\quad\quad\quad\quad\quad\quad\;\, W^{(1)}_{\alpha\beta} = \frac{\delta_{\alpha\beta}\Delta_\perp- \partial_\alpha\partial_\beta}{\epsilon_\textrm{ref}}, $$

(11)$$ L_{\alpha\beta} = \delta_{\alpha\beta}\Delta_\perp- \partial_\alpha\partial_\beta+\partial_\alpha\frac{1}{\epsilon_{zz}} \partial_\gamma\tilde\epsilon_{\gamma\beta}, \quad\quad\quad\quad \tilde\epsilon_{\alpha\beta} = \epsilon_{\alpha\beta}- \frac{\epsilon_{\alpha z}\epsilon_{z\beta}}{\epsilon_{zz}}. $$

In Eq. (8) and in the rest of the article, we abusively use the Bachmann-Landau notation $\mathcal {O}\left ({\eta ^{2}}\right )$ — usually used to characterize the asymptotic behavior of scalar functions — for differential operators using the following rule: for any differential operator $\mathbf {M}$, we write $\mathbf {M}=\mathcal {O}\left ({\mu }\right )$ with $\mu$ a small parameter if and only if $\mathbf {M}\boldsymbol {E}_\perp =\mathcal {O}\left ({\mu }\right )\boldsymbol {E}_\perp$.

The hypothesis (H1) can also be used to transform the second-order wave Eq. (3) into a simpler first-order Schrödinger-like equation $\left (i\partial _Z+\mathbf {T}\right )\boldsymbol {E}_\perp =\boldsymbol {0}$. To find the expression of the differential operator $\mathbf {T}$, we inject the formal solution of the previous Schrödinger equation into the original wave Eq. (3), and find that $\mathbf {T}$ must be a solution of the following quadratic equation:

(12)$$-\mathbf{T}^{2}+i\mathbf{S}\mathbf{T}+\mathbf{R}=\mathbf{0}.$$

Since $\mathbf {R}=\mathcal {O}\left ({1}\right )$ and $\mathbf {S}=\mathcal {O}\left ({\eta }\right )$ we can use a perturbative approach by first neglecting $\mathbf {S}$ and then finding a first order correction in $\mathbf {S}$ – which involves a continuous Lyapunov equation whose formal solution can be find in Ref. [13]. We find that two solutions are allowed:

(13)$$\mathbf{T}^{(\pm)}\equiv\pm\sqrt{\tilde{\boldsymbol\epsilon}+\mathbf{L}}+ i\mathbf{W}+i\mathbf{W}^{\prime\prime}+\mathcal{O}\left({\eta^{2}}\right),$$

where $\sqrt {}$ corresponds here to the functional square root. The functional square root of a differential operator $\mathbf {A}$ is defined as any operator $\mathbf {B}$ verifying $\mathbf {B}^{2}=\mathbf {A}$. Here, we impose the additional constraint that the eigenvalues of $\mathbf {B}$ must have positive real part. The differential operator $\mathbf {W}''$ is defined in Eq. (S2) of Supplement 1, and its complicated form should not frighten the reader, as we will show in a moment that similarly to $\mathbf {W}'$ it can be absorbed in a $\mathcal {O}\left ({\eta ^{2}}\right )$ term using an additional physical assumption. For now, let us take a bit of distance from the mathematics and try to summarize our results by giving physical meaning to the involved differential operators.

Since we found two solutions for the differential operators $\mathbf {T}$, there exists two class of propagating solution in the birefringent slabs introduced earlier: forward propagating modes $\boldsymbol {E}_\perp ^{(+)}$ and backward propagating modes $\boldsymbol {E}_\perp ^{(-)}$ which obey the following Schrödinger-like equation:

(14)$$\partial_Z\boldsymbol{E}_\perp^{(\pm)} = i\mathbf{T}^{(\pm)}\boldsymbol{E}_\perp^{(\pm)}.$$

These two types of mode are represented schematically in Fig. 1.

The definitions of the operators $\mathbf {T}^{(\pm )}$ in Eq. (13) involves three operators with very clear physical meaning:

• $\tilde {\boldsymbol \epsilon }$ is the dominant contribution to the evolution of the optical phase in Eq. (14), and for this reason will be called phase operator. It does not contain any derivatives and is therefore a pointwise matrix operator: $[\tilde {\boldsymbol \epsilon }\boldsymbol {E}_\perp ](\boldsymbol {r})$ only depends on the value of $\boldsymbol {E}_\perp$ at point $\boldsymbol {r}$. Interestingly, this operator also appears in the context of the well-known Jones method (see Eq. (2) in Ref. [14] with $k_x=k_y=0$).
• $\mathbf {L}$ is spectrally equivalent to a Laplacian (for homogeneous isotropic media, it can even be verified that $\mathbf {L}$ is the transverse Laplacian), and therefore will be called diffraction operator by analogy to the formalism of beam propagation in isotropic media [9]. It is responsible for the redistribution of energy in the transverse plane as light is propagated along $z$ (e.g. the natural spread of a Gaussian beam as it is transmitted through a homogeneous medium). It fully includes the effect of the anisotropy of the permittivity tensor up to order 1 in $\eta$.
• Finally, $\mathbf {W}$ is associated with terms of the type $\epsilon _{\alpha z}\partial _\beta$ in Eq. (14), which impose that energy and phase do not flow in the same direction when $\epsilon _{\alpha z}\neq 0$ – a typical property of birefringent media. For this reason, $\mathbf {W}$ will be called walk-off operator in the following. Similarly to the diffraction operator, $\mathbf {W}$ is accurate up to order 1 in $\eta$.

Figure 2 shows the typical physical effect associated with these three operators. As for the operator $\boldsymbol {W}''$ (which is defined as a function of $\mathbf {W}'$), it can be viewed as a correction to the walk-off operator $\boldsymbol {W}$ which starts becoming relevant when the optical fields have non-negligible high-frequency spatial components and when $\epsilon _{xz}$ and $\epsilon _{yz}$ have sharp variations or discontinuities (in which case the commutator terms in the definition of $\mathbf {W}'$ and $\mathbf {W}''$ becomes non-negligible). This can be shown by expanding the definition of $\boldsymbol {W}''$ in Supplement 1 as a series of derivatives of increasing order, and observing that the lowest order for the derivatives is 3 – contrary to $\mathbf {W}$ which is associated with a first-order derivative at the lowest order.

Fig. 2. (a) Color-graded representation of the real part of $E_y$ when a plane wave propagates through a system of birefringent slabs with homogeneous permittivity tensor. In this particular system, our formalism in the main text is fully equivalent to the simple Jones method and solely relies on the phase operator $\tilde {\boldsymbol \epsilon }$. (b) Same as (a) with a Gaussian beam propagating through a uniaxial media with tilted optical axis $\boldsymbol {n}$; the Poynting vector $\boldsymbol {S}_p$ is not aligned with the wavevector $\boldsymbol {k}$ due to the action of the walk-off operator $\mathbf {W}$ in the propagation Eq. (14). (c) Natural spread of the intensity of a Gaussian beam in a homogeneous birefringent media, which is physically modeled by the diffraction operator $\mathbf {L}$ in the propagation Eq. (14). In (a,b,c), the white bar represents $0.5$ µm, and in (a,c) the optical axis is normal to the plane of the figure.

Download Full Size | PDF

Since the numerical implementation of $\boldsymbol {W}''$ is not possible using sparse matrices – contrary to $\tilde \epsilon$, $\mathbf {L}$, and the operators $\mathbf {W}^{(0,1)}$ defining $\mathbf {W}$ – we will chose to neglect this operator as mentioned previously either by assuming that the optical fields do not have high-frequency components (paraxial regime of propagation, see Sec. 2.2) or by assuming arbitrary optical fields and smooth variations of $\epsilon _{xz}$ and $\epsilon _{yz}$ (wide-angle regime of propagation, see Sec. 2.3).

2.2 Paraxial scheme

In this subsection, we assume that the optical fields propagate paraxially along $z$. Based on the angular spectrum approach [15], this assumption is fully equivalent to the following physical limitation:

(H2.a) Optical fields have negligible high-frequency spatial components in the $x$ and $y$ directions (whose typical order of magnitude can be evaluated using the uniform norm of the transverse gradient $\nabla _\perp =\{\partial _x,\partial _y\}$ of the transverse field): $\rho \equiv \left \lVert {\nabla _\perp \boldsymbol {E}_\perp }\right \rVert _\infty /\left [k_0\left \lVert {\boldsymbol {E}_\perp }\right \rVert _\infty \right ]\ll 1$.

Using hypothesis (H2.a) and the main result of Ref. [16] to expand $\sqrt {\tilde {\boldsymbol \epsilon }+\mathbf {L}}$ around $\sqrt {\tilde {\boldsymbol \epsilon }}$, we find that the expressions of the forward and backward propagation operators $\mathbf {T}^{(\pm )}$ are simplified to:

(15)$$\mathbf{T}^{(\pm)}=\pm\left(\sqrt{\tilde{\boldsymbol\epsilon}}+\mathbf{L}'\right)+i\mathbf{W}^{(0)}+ \mathcal{O}\left({\eta^{2},\rho^{3}}\right),$$

where we introduced the rescaled diffraction operator $\mathbf {L}'=\left [ \tilde {\boldsymbol \epsilon }^{-1/2}\mathbf {L}+\mathbf {L}\tilde {\boldsymbol \epsilon }^{-1/2}\right ]/4$.

Here, we should emphasize that the obtained paraxial beam propagation operator in Eq. (15) does not depend on the choice of a reference index – contrary to most paraxial BPM in isotropic media [1,9] and anisotropic media [5,7]. In all of these cited documents, the choice of a reference index $n_\textrm {ref}$ for the phase evolution is done in order to skip the use of a square root in the beam propagation equation, but can lead to inaccurate results when multiple modes with different effective indices propagate in the considered structure [9]. In our formalism, the $\sqrt {\tilde {\boldsymbol \epsilon }}$ terms in the definition of $\mathbf {L}'$ can be viewed as a structure-dependent local reference index allowing to get an accurate diffraction of the fields. In addition and as already mentioned in Fig. 2, the phase evolution based on $\sqrt {\tilde {\boldsymbol \epsilon }}$ in Eq. (15) is fully equivalent to the Jones method, which gives exact solution of Maxwell equations for plane waves propagating in homogeneous birefringent systems. This is a major advantage in comparison to schemes using a reference index, since such schemes cannot provide an exact evolution of the optical phase even in homogeneous birefringent systems.

Numerically, optical fields are propagated in each slab using a Strang splitting of the evolution operators associated with Eq. (14), accurate at order 3 [17] in the renormalized slab thickness $\Delta _Z=k_0\Delta _z$ ($\Delta _z$ being the slab thickness). For the sake of brevity, we only give the evolution equation for forward-propagating modes (similar expressions for the backward propagating modes can be obtained by using the definition of $\mathbf {T}^{(-)}$ in Eq. (15)):

(16)$$\boldsymbol{E}^{(+)}_\perp\big|_{Z+\Delta_Z} = \mathbf{P}^{(\epsilon)}\mathbf{P}^{(\rm diff)}\mathbf{P}^{(\epsilon)} \boldsymbol{E}^{(+)}_\perp\big|_Z + \mathcal{O}\left({\eta^{2},\rho^{3},\Delta_Z^{3}}\right),$$

where the evolution operators $\mathbf {P}^{(\epsilon )}$ and $\mathbf {P}^{(\rm diff)}$ are defined as:

(17)$$\mathbf{P}^{(\epsilon)} = \exp\left[\frac{i\Delta_Z}{2}\sqrt{\tilde{\boldsymbol\epsilon}}\right], \quad\quad\quad\quad \mathbf{P}^{(\rm diff)} = \exp\left[i\Delta_Z\left(\mathbf{L}'+i\mathbf{W}^{(0)}\right)\right].$$

The phase evolution operator $\mathbf {P}^{(\epsilon )}$ can be calculated analytically using the Cayley-Hamilton theorem since $\tilde {\boldsymbol \epsilon }$ is a pointwise 2x2 matrix operator. The numerical implementation of the evolution operator $\mathbf {P}^{(\rm diff)}$ is more tricky but can nevertheless be done efficiently using a combination of finite-difference discretization, (1,1) Padé approximant for the exponential, and alternative direction implicit (ADI) splitting identical to Ref. [9]. The advantage of the ADI method, accurate at order $3$ in $\Delta _Z$ in this context, is that it is associated with a very efficient linear complexity $\mathcal {O}(N)$ (with $N$ the number of degrees of freedom in the transverse place) based on the Thomas matrix inversion algorithm. Note that our numerical implementation includes transparent boundary conditions as described in Ref. [18] in order to avoid reflections on the sides of the mesh.

2.3 Wide-angle scheme

We now examine the wide-angle regime of propagation by removing any constraints on the variations of the optical fields (which are therefore allowed to propagate at any angle with respect to $z$) and replacing (H2.a) with the following physical assumption:

(H2.b) There are no discontinuities or sharp variations of $\epsilon _{xz}$ and $\epsilon _{yz}$ in the $x$ and $y$ directions:
$\forall \;\alpha \in \{x,y\}$, $\Delta _\perp \epsilon _{\alpha z}-\epsilon _{\alpha z}\Delta _\perp =\mathcal {O}\left ({\eta ^{2}}\right )$.

Using (H2.b) allows us to absorb the complicated operator $\mathbf {W}''$ (correction to the walk-off operator) in the $\mathcal {O}\left ({\eta ^{2}}\right )$ term in Eq. (13):

(18)$$\mathbf{T}^{(\pm)}\equiv\pm\sqrt{\tilde{\boldsymbol\epsilon}+\mathbf{L}}+ i\mathbf{W}+\mathcal{O}\left({\eta^{2}}\right).$$

The hypothesis (H2.b) put a rather stringent constraint on the optical axis orientation at interface of discontinuity of the permittivity. Nevertheless, for LC-based photonics devices based on simple planar or homeotropic anchoring at the confining parallel plates, (H2.b) can be always made true by rotating the sample in the simulation such that $\epsilon _{xz}=\epsilon _{yz}=0$ at the sample plate interfaces. For more complicated systems in which such a simplification is not possible (e.g. patterned anchoring such as in Ref. [19], or localized defects with nonzero $\epsilon _{xz}$ and/or $\epsilon _{yz}$ inside the core), hypothesis (H2.b) is broken and wide-angle corrections to the walk-off operator will be inaccurately modeled if our scheme is still used, but the other operators for the diffraction and phase evolution are still accurately modeled. More precisely, an error of order $\mathcal {O}\left ({\eta '\rho ^{3}}\right )$ will be introduced in Eq. (18) if hypothesis (H2.b) is not fulfilled, with $\eta '=\max \left [\left \lVert {\epsilon _{xz}}\right \rVert _\infty , \left \lVert {\epsilon _{xz}}\right \rVert _\infty \right ]/\epsilon _\textrm {ref}$ and $\rho \equiv \left \lVert {\nabla _\perp \boldsymbol {E}_\perp }\right \rVert _\infty / \left [k_0\left \lVert {\boldsymbol {E}_\perp }\right \rVert _\infty \right ]$. For localized structures such as defects, we expect that the impact of this error term will be negligible if the beam size is much bigger than the width of the defect.

Similarly to the paraxial scheme, the evolution equation for the optical fields relies on a Strang splitting, which we only give for forward-propagating modes for brevity:

(19)$$\boldsymbol{E}^{(+)}_\perp\big|_{Z+\Delta_Z} = \mathbf{P}^{(w}\mathbf{P}^{(r)}\mathbf{P}^{(w)} \boldsymbol{E}^{(+)}_\perp\big|_Z + \mathcal{O}\left({\eta^{2},\Delta_Z^{3}}\right),$$

where the evolution operators $\mathbf {P}^{(w)}$ and $\mathbf {P}^{(r)}$ are defined as:

(20)$$\mathbf{P}^{(w)} = \exp\left[-\frac{\Delta_Z}{2}\mathbf{W}\right], \quad\quad\quad\quad \mathbf{P}^{(r)} = \exp\left[i\Delta_Z\sqrt{\tilde{\boldsymbol\epsilon}+\mathbf{L}}\right].$$

The numerical implementation of $\mathbf {P}^{(w)}$ is based on finite-difference discretization and a (1,1) Padé approximant of the exponential accurate at order 3 in $\Delta _Z$:

(21)$$\mathbf{P}^{(w)} = \left[\mathbf{I}+\mathbf{W}^{(1)}+\frac{\Delta_Z}{4}\mathbf{W}^{(0)}\right]^{-1} \left[\mathbf{I}+\mathbf{W}^{(1)}-\frac{\Delta_Z}{4}\mathbf{W}^{(0)}\right] + \mathcal{O}\left({\Delta_Z^{3}}\right).$$

The numerical implementation of $\mathbf {P}^{(r)}$ is also based on finite-difference discretization and a (1,1) Padé approximant for the exponential, supplemented by the rotated branch cut rational approximation of the square root [20] and the factorization of the Padé approximant into simpler terms linear in $\mathbf {R}\equiv \tilde {\boldsymbol \epsilon }+\mathbf {L}$ [21]. The action of matrix inverses of the form $\left (\mathbf {I}+\mathbf {A}\right )^{-1}$ (where $\mathbf {A}$ can be a linear combination of $\mathbf {W}^{(0)}$ and $\mathbf {W}^{(1)}$ or a term proportional to $\mathbf {R}$) are iteratively calculated using a BiCGStab solver with a naive preconditioner of the form $\mathbf {I}-\mathbf {A}$. We recall that preconditionners are mathematical tools allowing to boost the convergence of iterative matrix inverse solver by providing an easy-to-calculate estimation of the matrix inverse; here, our naive preconditioner simply corresponds to a first-order Taylor expansion of matrix inverses of the form $\left (\mathbf {I}+\mathbf {A}\right )^{-1}$. In all numerical experiments performed for this paper, the order of the rotated branch cut rational approximation of the square root functional was 4, which ensured an accurate evaluation of the evolution operator $\mathbf {P}^{(r)}$. Similar to the paraxial scheme, we include transparent boundary conditions in the numerical implementation of the evolution operators.

2.4 Interface conditions for the optical fields

In the last two sub-sections, we presented the general formalism allowing to propagate optical fields inside each birefringent slab of Fig. 1 in two separate regimes – paraxial and wide-angle. Since discontinuity jumps of the permittivity are allowed at each interface between the birefringent slabs, all that remains to do is formulate the continuity equations for the optical fields at these interfaces. Such continuity equations allows the calculation of the forward/backward modes of the $(n+1)$-th slab as a function of the forward/backward modes of the $n$-th slab. From the continuity of the transverse electric and magnetic fields at a slab interface, Maxwell-Faraday equation ($ik_0c\boldsymbol {B}=\nabla \times \boldsymbol {E}$) and Eqs. (4) and (14), we obtain the following set of equations:

(22)$$ \boldsymbol{E}^{(+)}_{n}+\boldsymbol{E}^{(-)}_{n} = \boldsymbol{E}^{(+)}_{n+1}+\boldsymbol{E}^{(-)}_{n+1}, $$

(23)$$ \mathbf{N}^{(+)}_{n}\boldsymbol{E}^{(+)}_{n}-\mathbf{N}^{(-)}_{n}\boldsymbol{E}^{(-)}_{n} = \mathbf{N}^{(+)}_{n+1}\boldsymbol{E}^{(+)}_{n+1}-\mathbf{N}^{(-)}_{n+1}\boldsymbol{E}^{(-)}_{n+1}, $$

where a subscript $n$ (resp. $n+1$) denote a quantity associated with the $n$-th slab (resp. $(n+1)$-th slab). The operators $\mathbf {N}^{(\pm )}$ are defined in Eqs. (S3) and (S6) of Supplement 1, here we will simply mention that they simplify to $\sqrt {\tilde {\boldsymbol \epsilon }}$ (i.e. the phase operator) when the optical and permittivity fields are invariant by translation in the $x$ and $y$ directions.

Physically speaking, the operators $\mathbf {N}^{(\pm )}$ in Eqs. (22) and (23) play the role of effective index in the direction normal to the interface for the forwards and backward propagating modes, and allow to generalize the well-known Fresnel equations for arbitrary optical fields (not simply plane waves) at interfaces between arbitrary anisotropic media (to the best of our knowledge, such a generalisation was only done for isotropic media, see for example Ref. [9]). For example, assuming an interface between two anisotropic media 1 (input media) and 2 (output media) with an incident ($\boldsymbol {E}_i$), reflected ($\boldsymbol {E}_r$) and transmitted field ($\boldsymbol {E}_t$), we calculate from Eqs. (22) and (23):

\begin{aligned} \boldsymbol{E}_t &= \left[\mathbf{N}^{(+)}_2+\mathbf{N}^{(-)}_1\right]^{-1} \left[\mathbf{N}^{(+)}_1+\mathbf{N}^{(-)}_1\right]\boldsymbol{E}_i, \\ \boldsymbol{E}_r &= \left[\mathbf{N}^{(+)}_2+\mathbf{N}^{(-)}_1\right]^{-1} \left[\mathbf{N}^{(+)}_1-\mathbf{N}^{(+)}_2\right]\boldsymbol{E}_i, \end{aligned}

which directly gives the well-known Fresnel amplitude coefficients $t=2n_1/(n_1+n_2)$ (transmission) and $r=(n_1-n_2)/(n_1+n_2)$ (reflection) in the case of a normally incident plane wave and isotropic media of indices $n_{1}$ and $n_{2}$.

Note that the current numerical implementation of our BPM schemes takes into account Eqs. (22) and (23) only at the entrance and exit of the birefringent sample (interfaces between an isotropic external medium and the birefringent medium), where we make the approximation $\boldsymbol {N}^{(\pm )}\approx \sqrt {\tilde {\boldsymbol \epsilon }}$ mentioned above (no significant variation of the fields in the $x$ and $y$ directions) and where we assume zero backward-propagating fields on the output side of the interface. In the bulk of the birefringent medium we assume $\boldsymbol {E}_n^{(+)}\approx \boldsymbol {E}_{n+1}^{(+)}$, i.e. that the birefringent medium do not significantly contribute to reflections. This approximation is valid only if the permittivity tensor is smoothly varying along the $z$ direction, which we will assume in the remainder of this article. In the future, it would be interesting to fully implement Eqs. (22) and (23) in the bulk of the birefringent sample by using the bi-directional beam propagation scheme described in Ref. [9], thereby allowing to model light propagation in sample with photonics bandgap or significant reflective properties (e.g. Fabry-Perrot or Distributed Bragg mirrors).

3. Numerical results and discussion

Now that our two beam propagation schemes — paraxial and wide-angle — for birefringent media are introduced, we discuss associated numerical results for liquid-crystal-based photonics systems. We first validate our methods by calculating the computational error inside homogeneous birefringent samples and one-dimensional diffraction gratings, and then show a real-life example of application of our numerical code by demonstrating that topological solitons in frustrated cholesterics can be used as light waveguides.

3.1 Validation of the scheme

To validate the two schemes introduced in Sec. 2, we calculate the computational error between numerically calculated optical fields and reference solutions on three systems:

• System A: homogeneous birefringent layer with optical axis $\boldsymbol {n}=\boldsymbol {e}_x$ and input beam polarisation $\boldsymbol {e}_x+\boldsymbol {e}_y$;
• System B: 1D diffraction grating with in-sample-plane rotation of the optical axis $\boldsymbol {n}=\cos (2\pi y/p)\boldsymbol {e}_x+\sin (2\pi y/p)\boldsymbol {e}_y$ and input beam polarisation $\boldsymbol {e}_x+i\boldsymbol {e}_y$;
• System C: 1D diffraction grating with cholesteric-like rotation of the optical axis $\boldsymbol {n}=\cos (2\pi y/p)\boldsymbol {e}_z+\sin (2\pi y/p)\boldsymbol {e}_x$ and input beam polarisation $\boldsymbol {e}_x$.

For each system, the wavelength in empty space is $\lambda =500~$nm and the input profile for the transverse optical fields is a diffraction-limited highly-focused profile for system A and a Gaussian profile with a waist of $1$ µm for system B and C. We assume that $x$ is an axis of invariance for the optical and permittivity fields and propagate the optical fields along the $z$-axis on a computational box of sizes $L_y=6$ µm and $L_z=3$ µm. The computational error $\mu$ is calculated with the following formula:

(24)$$\mu\equiv||\boldsymbol{E}^\textrm{(num)}_\perp-\boldsymbol{E}^\textrm{(ref)}_\perp||_2/ ||\boldsymbol{E}^\textrm{(ref)}_\perp||_2,$$

where an index “(num)” (resp. “(ref)”) indicates a numerically calculated (resp., reference) solution. For system A, the reference solution is calculated analytically from Maxwell equations in Fourier space (which is possible because system A is homogeneous). For system B and C, the reference solution is calculated from a Finite-Difference-Time-Domain (FDTD) simulation on a very fine mesh ($\Delta _y=\Delta _z=6~$nm), since no known analytical solution are readily available on these systems. In every simulations, the permittivity tensor is assembled from the optical axis field assuming uniaxial birefringence: $\boldsymbol \epsilon = n_o^{2}\mathbf {I}+(n_e^{2}-n_o^{2})\boldsymbol {n}\otimes \boldsymbol {n}$, with $n_e=1.75$ ($n_o=1.5$) the extraordinary (ordinary) index. All FDTD simulations were performed using the open-source software MEEP developed by the MIT [22] with MPI parallelisation. Paraxial and wide-angle beam propagation simulations were run on a custom-developed C++ code based on the schemes of Sec. 2 and parallelized with OpenMP. Transparent (resp., perfectly matched layer) boundary conditions ensured that outgoing fields were absorbed on the side of the computational box in BPM (resp., FDTD) simulations. Finally, analytical solutions of Maxwell equations in Fourier space were calculated directly using the Python library numpy.

The error $\mu$ can include two contributions: numerical error due to the finite-difference discretization of BPM or Yee-lattice discretization of FDTD, and scheme error due to the approximations (H1) and (H2) made in Sec. 2. For system A and B, the latter error is exactly zero for our wide-angle scheme because the optical axis is $z$-independent and always orthogonal to the propagation axis $z$; this can be verified either by repeating the calculations of Sec. 2 with $\epsilon _{\alpha z}=0$ and observing that no error terms are needed to reach the final propagation equation, or simply comparing our final propagation equation with the exact wave equation in planar liquid-crystal structures [23]. Note that this means that whenever the walk-off operator $\mathbf {W}$ is exactly zero and the permittivity tensor is $z$-independent, our wide-angle scheme is expected to be identical to the exact wave equation in birefringent media without having to assume hypotheses (H1) and (H2.b). We conclude that in system A and B, the computational error in wide-angle simulations is only due to the numerical discretization, while for system C — for which the walk-off operator is nonzero — the error $\mu$ also includes a $\mathcal {O}\left ({\eta ^{2}}\right )$ in addition to the discretization error, as shown in Eq. (19). Paraxial simulations always include a scheme error in addition to the discretization error since for this class of BPM, hypothesis (H2.a) is always needed whatever the system in consideration.

In Fig. 3, we show color-graded images of the optical field solutions and plots of the computational error as a function of the mesh spacing $\Delta _y$ for systems A, B and C. Note that we impose $\Delta _z=\Delta _y/3$ for beam propagation simulations and $\Delta _z=\Delta _y$ for FDTD simulations. We used smaller axial steps for BPM because we empirically noticed that the naive preconditioners of the wide-angle scheme were more efficient in this case. Figure 3(d) shows that the numerical error scales quadratically with the mesh spacing for the FDTD and wide-angle beam propagation simulation, as expected for centered finite-difference and Yee-lattice discretization [22]. In Fig. 3(e), this quadratic scaling is only seen transiently for systems B and C because the error $\mu$ includes constant contributions from the reference FDTD simulation (which is only approximate even with very fine computational meshes) and from the scheme error as discussed above for system C. As for paraxial simulations, they are always less accurate than the wide-angle simulations and saturates for small mesh spacings because they rely on an approximation (see hypothesis (H2.a)). Nevertheless, our two schemes allows to reach computational errors comparable to $1$–$5~\%$, even in systems where the propagation equations derived in Sec. 2 are not exact but rely on approximations. This validates the accuracy of our two schemes.

Fig. 3. (a,b,c) Color-graded plot of the real part of $E_x$ in systems A, B and C respectively. The optical axis variations are represented schematically with cylinders below the plots, and the axes orientations on the right apply to all plots. The white bars represent $0.5$ µm. (d) Computational error $\mu$ as a function of the mesh spacing $\Delta _y$ for system A. The reference solution for the calculation of the error is analytical. (e) Same as (d) for systems B and C. The reference solution is calculated from a FDTD simulation on a very fine mesh. In (d,e), dashed (dotted) curves correspond to wide-angle (paraxial) beam propagation simulations and the solid curve correspond to a FDTD simulation.

Download Full Size | PDF

We did not make a full theoretical analysis of the stability of our schemes, but can at least point out that our paraxial scheme was always found stable in our numerical experiments (including a wide range of geometries such as topological solitons, non-local optical solitons, cholesteric droplets, photo-patterned liquid crystal samples$\ldots$). This is expected since this scheme is based on an exact calculation of the phase evolution and a mature ADI scheme already widely used for isotropic media [9]. The situation is a bit more complex for our wide-angle BPM, which necessitates some tweaking of the scheme parameters to get stable propagation when the mesh is very fine (typically $\Delta _y<0.1$ µm). We believe the stability of this scheme could be improved with better preconditioning for the iterative solvers and more robust absorbing boundary conditions such as perfectly matched layers [24].

We conclude this section by giving typical computational running times of our numerical schemes. With a threshold accuracy of $1~\%$, the typical running time for our wide-angle scheme is $\sim 10~$s for the systems studied here, while the running time of the FDTD code MEEP was $\sim 1300~$s. These times were obtained on a state-of-the-art computational station with 12 CPU cores (i7-7800X), and are typically multiplied by a factor $10$ for three-dimensional simulations. Of course, these results highly depend on the hardware on which the codes are running, but as a general rule, our wide-angle beam propagation scheme runs $\sim 100\times$ faster than the FDTD code MEEP, and our paraxial scheme is $\sim 5\times$ faster than our wide-angle scheme — but is limited in terms of accuracy as shown in Figs. 3(d) and (e). We deduce that our two schemes are particularly suited for efficient BPM in anisotropic media, as long as the optical fields propagates dominantly along $z$. Typically, we found that our paraxial scheme can reliably propagates fields up to an angle of deviation of $5^{\circ}$ with respect to $z$, while our wide-angle scheme is accurate up to angles of $20-30^{\circ}$, depending on the mesh spacing. Thus, users of our beam propagation schemes should always choose the reference propagation axis $z$ in order to minimize angular deviation and keep these angular limits in mind when evaluating the accuracy of their simulations.

We remark that these angular limits for both our schemes were obtained from the simulations on system A by imposing a threshold of accuracy of $5~\%$ and calculating the errors of individual transverse Fourier modes, which correspond to tilted plane waves in a homogeneous birefringent slab. In principle, higher angles of deviation could be covered with higher-order Padé approximant in our wide-angle scheme; however, for extremely fine mesh and highly non-paraxial fields, the simple preconditioners used in our wide-angle scheme are not robust enough to efficiently invert the matrices of propagation. This technical limitation could be lifted in principle with better preconditioning.

3.2 Waveguides based on topological solitons

Now that the validity of our code was established on simple 2D simulations, we turn our focus to a more complex 3D system. The goal of this section is to show that topological solitons in frustrated cholesterics can be used as light waveguides. In usual unconfined cholesterics — a chiral liquid crystal phase with orientational order — the director field (which here is identical to the optical axis $\boldsymbol {n}$) rotates around a unique spatial direction, thus forming the so-called cholesteric helix texture. When confined between two plates treated for an homeotropic anchoring ($\boldsymbol {n}$ constrained to be normal to the plate), the equilibrium cholesteric helix cannot be formed and more complicated states emerge. Among these frustrated states, we will focus here on cholesteric fingers of the second kind (CF2), which correspond to line-like structures embedded in a uniform homeotropic background [25]. These structures can be addressed as topological solitons because they cannot be continuously deformed into the homogeneous homeotropic state, and typically appears in samples with a ratio $h/P\sim 1$, with $h$ the sample thickness and $P$ the cholesteric pitch — the periodicity of the equilibrium cholesteric helix texture in unconfined domains. We will first consider straight CF2 associated with an axis of invariance, and will give some preliminary results with curved CF2 at the end of this section.

Since CF2s are associated with a localized modulation of the permittivity tensor $\boldsymbol \epsilon$, they make good candidates for the guiding of light. In Fig. 4(a), we show the director field profile of a straight CF2 in a plane orthogonal to the axis of invariance $z$, as calculated from a minimization of the free energy of the cholesteric phase [26] using the material constants of the liquid crystal E7 at ambient temperature [27,28]: $K_1=11~$pN, $K_2=7~$pN, $K_3=18~$pN, $n_e=1.746$, $n_o=1.522$. We also assumed that the confining plates sandwiching the cholesteric layer had a refractive index $n_p=1.51$, as measured for typical crown glasses [29]. In all simulations of this section, we always impose $h/P=1$. Note that since CF2s do not contain any topological defects, they are scale-invariant and only depends on the ratio $h/P$. As a consequence, the profile of Fig. 4(a) can be used to model arbitrarily-thick samples with the appropriate scaling operation as long as the sample thickness and cholesteric pitch fulfill $h/P=1$.

Fig. 4. (a) Representation of the transverse profile of the director field of a $z$-invariant CF2 with cylinders. The thickness of the liquid crystal layer is $h$ and homeotropic boundary conditions are imposed on the top and bottom confining surface along $x$. Light eigenmodes are calculated inside the dashed region. Color-graded representation of $\textrm {Re}~E_x$ (b), $\textrm {Im}~E_x$ (c), $\textrm {Re}~E_y$ (d) and $\textrm {Im}~E_y$ (e) are shown below for the fundamental eigenmode in a $30$ µm-thick sample.

Download Full Size | PDF

We calculated eigenmodes of the wide-angle propagation operator $\mathbf {T}^{(+)}$ defined in Eq. (18) using the ARPACK++ library. We restricted the calculation of eigenmodes in the dashed region of Fig. 4(a) by imposing zero field values outside this region, for a reason that will become clear later on. A typical result of this calculation for a sample of thickness $h=30$ µm is shown in Figs. 4(b)–(e), where we plotted color-graded representations of the transverse optical field components for the fundamental eigenmode, i.e. the eigenmode with the highest eigenvalue. As visible, light is confined inside a banana-shaped domain which correspond to the dark-blue cylinders in Fig. 4(a). The polarisation of this eigenmode is mostly-aligned with the director field along $\boldsymbol {e}_x$, as expected since the birefringence is positive here ($n_e>n_o$) and since index-gradient type waveguides rely on permittivity profile with a negative curvature — in other words, the effective refractive index must be at its highest value at the center of the guided mode.

However, we emphasize that the eigenmodes calculated inside the dashed domain of Fig. 4(a) are not eigenmodes of the full liquid crystal sample because of the artificial zero-field boundary conditions that we imposed. When calculating eigenmodes on the complete computational mesh rather than the localized dashed domain of Fig. 4(a), we noticed that the optical fields always contained a small fraction of intensity scattering out of the CF2, which suggest that the calculated eigenmodes are not strictly speaking waveguide modes but rather leaky modes, i.e. confined modes of light which slowly loose energy when propagating along the axis of invariance $z$. This suggestion is reinforced by the fact that the $\epsilon _{xx}$-profile along the $y$-axis in the mid-plane of the sample has a W-shape, a specific shape which is known to lead to leaky waveguides and leaky modes in isotropic media [30]. To definitely demonstrate that the fundamental eigenmode shown in Figs. 4(b)–(e) is indeed a leaky mode, we used it as input for a typical wide-angle beam propagation simulation along the $z$-axis. By denoting $\boldsymbol {E}_\perp ^{(0)}$ the transverse optical field at the input, we then calculated the fraction $\phi (z)$ of lost intensity for the fundamental mode after a distance of propagation $z$:

(25)$$\phi(z)\equiv 1-\left|\frac{\int\boldsymbol{E}_\perp(x,y,z)^{*}\cdot\boldsymbol{E}_\perp^{(0)}(x,y){\textrm{d}}x{\textrm{d}}y} {\int\boldsymbol{E}_\perp^{(0)}(x,y)^{*}\cdot\boldsymbol{E}_\perp^{(0)}(x,y){\textrm{d}}x{\textrm{d}}y}\right|^{2}.$$

In Fig. 5(a), we plot $\phi$ as a function of $z$ for samples of different thicknesses. After a transient regime delimited by a dashed line on Fig. 5(a), the fraction of lost intensity $\phi$ is approximately linear in $z$, which indicates that the fundamental mode of Figs. 4(b)–(e) is linearly loosing energy. The magnitude of this energy leaking depends on the sample thickness: the thinner the sample, the higher the linear loss. To better characterize this thickness-dependence effect, we introduce the linear scattering loss $\Gamma$:

(26)$$\Gamma\equiv-\frac{{\textrm{d}}}{{\textrm{d}}z}\left[10\log_{10}(1-\phi)\right].$$

Note that $\Gamma$ can be well-approximated by $(10/\log 10)({\textrm {d}}\phi /{\textrm {d}}z)$ since $\phi \ll 1$. In Fig. 5(b), we plot $\Gamma$ as a function of the sample thickness $h$. As shown, numerical data points are well-fitted with an exponential law of the type $\Gamma _0\exp (-h/l_0)$, with $\Gamma _0\approx 94.5~$db/cm and $l_0\approx 4.28$ µm. Note that this law is very similar to the well-known expression of the transmission coefficient of a quantum particle tunneling through a rectangular potential barrier. We also remark that the typical order of magnitude of the scattering loss $\Gamma$ is small with respect to losses due to thermal orientational fluctuations of the director field (typically $40~$db/cm [31]), which indicates that the leaky modes calculated here should in principle be visible in actual experiments.

Fig. 5. (a) Fraction of lost intensity of the fundamental eigenmode as a function of propagation distance $z$ inside a straight CF2. From right to left, $h=10$–$30$ µm by increments of $2$ µm. The dashed line symbolizes the transition to the linear loss regime. (b) Scattering loss $\Gamma$ as a function of sample thickness $h$. The error bars are obtained from the standard deviation of the slope of $f$ above the dashed line of (a). The solid line correspond to a fit with an exponential law.

Download Full Size | PDF

We close this section by providing an outlook on possible applications of CF2-based light waveguides. Up until now, we focused on CF2 with an axis of invariance along $z$, but more complex geometries are entirely possible since Ackerman et al. developed an experimental technique allowing to write in real-time cholesteric fingers with arbitrary shapes in the $yz$ plane using scanning lasers [32]. This technique makes CF2s a very attractive target for writing in real-time photonics circuits such as directional couplers, provided that losses due to curvature of the CF2 are not too important. We provide here a few preliminary results concerning the role of the curvature which could guide future work in this line of research. In Fig. 6, we show how light is guided inside sine-like CF2s with various in-plane curvature in a sample of thickness $h=30$ µm. In this figure, each bended section of the CF2 has a length $\pi R/9$ with $R$ the in-plane radius of curvature. The optical micrographs in Figs. 6(a)–(c) were simulated using the open-source software Nemaktis (see https://github.com/warthan07/Nemaktis) which is based on our paraxial BPM to propagate optical fields in virtual birefringent samples as in a real polarised optical microscope. As shown in Figs. 6(d)–(f), light is well confined as long as the radius of curvature of the CF2 is greater than $\sim 2h$, but this lower bound can be extended to smaller radius of curvature if the length of the bended part of the CF2 is reduced. We also expect that the quality of light confinement should be greatly increased when using high-birefringence liquid crystals.

Fig. 6. (a–c) Simulated crossed-polarizers optical micrographs of sine-like CF2s with a typical radius of curvature of $100$ µm, $70$ µm and $40$ µm respectively. (d–f) Associated $x$-averaged beam intensity when the fundamental eigenmode of Figs. 4(c)–(e) is sent inside the CF2s of (a–c). The white bars represent $30$ µm.

Download Full Size | PDF

4. Conclusion

We introduced a general theoretical framework for modeling beam propagation in arbitrary birefringent media, and derived a paraxial and wide-angle BPM from this framework assuming two reasonable mathematical hypotheses based on physical limitations (small index contrast, paraxial propagation or smooth variations of the longitudinal components of the permittivity tensor). One of the main attractive feature of our methods is that they are based on 2D differential operators with a clear physical interpretation and easy numerical implementation based on sparse matrices. The current numerical implementation of our schemes is only based on forward-propagating modes, but we also theoretically showed how the continuity equations for the electric and magnetic fields at interfaces of discontinuity of the permittivity can be taken care of.

We validated the accuracy of our schemes by showing that the computational error scales quadratically with the mesh spacing in systems where our beam propagation equations are expected to be exact, and showed that the computational error is still reasonably low in systems where our framework is only approximate. Finally, we applied our wide-angle beam propagation scheme to a promising system for tunable light waveguiding and manipulation, namely topological solitons in frustrated cholesteric liquid crystals. We showed that light is efficiently guided along straight topological soliton with low scattering loss, and examined how the curvature of these topological solitons — whose trajectory can be written on-demand with scanning lasers — affects the confinement of light. This study therefore provide the first step towards reconfigurable topological-soliton-based photonics components such as directional couplers.

We also underline once again that when the walk-off terms are inexistent ($\epsilon _{\alpha z}=0$, i.e. the optical axis is orthogonal or parallel to $z$), our wide-angle scheme is rigorously exact up to discretization errors. A very interesting consequence of this observation, as brought forward by one of the referee, is that our wide-angle scheme can also be used to model high index contrast LC devices (e.g. with silicon cladding) as long as the hypotheses (H1) and (H2.b) of Sec. 2 are fulfilled in the bulk of the liquid crystal and the optical axis is almost orthogonal or parallel to $z$ at the interface between the LC and the confining isotropic media.

We finally mention possible extensions and applications of our work. First, implementing bi-directional BPM based on the continuity equations of Sec. 2.4, as mentioned previously, would facilitate the study of photonic bandgap in birefringent media. Other possible extensions of our work include more efficient preconditionners for the evolution operators — an interesting challenge for specialists in numerical analysis — and perfectly matched layers to better absorb outgoing fields on the side of the computational box [24]. Second, we mention that our schemes can also be used to model nonlinear optics based on the elastic response of liquid crystals to optical perturbations [23,33]. Last, we emphasize that our paraxial scheme is particularly suited for microscopy simulations, since typical microscopy objectives with low numerical aperture filter-out spatial high-frequencies — which therefore do not need to be modeled accurately. We plan to provide a detailed description of this application to microscopy in a future paper.

Funding

Javna Agencija za Raziskovalno Dejavnost RS (P1-0099); Horizon 2020 Framework Programme (CA17139, MSCA 834256).

Acknowledgements

The authors acknowledge funding from the ARSS (Javna Agencija za Raziskovalno Dejavnost RS) through grant P1-0099 and from the European Union’s Horizon 2020 program through the Marie Skłodowska-Curie grant agreement No. 834256 and COST action EUTOPIA CA17139

Disclosures

The authors declare no conflicts of interest.

See Supplement 1 for supporting content.

References

1. J. Yamauchi, J. Shibayama, and H. Nakano, “Modified finite-difference beam propagation method based on the generalized Douglas scheme for variable coefficients,” IEEE Photonics Technol. Lett. 7(6), 661–663 (1995). [CrossRef]

2. K. Fan, W. Cai, and X. Ji, “A full vectorial generalized discontinuous Galerkin beam propagation method (GDG-BPM) for nonsmooth electromagnetic fields in waveguides,” J. Comput. Phys. 227(15), 7178–7191 (2008). [CrossRef]

3. J. Lopez-Dona, J. Wanguemert-Perez, and I. Molina-Fernandez, “Fast-Fourier-based three-dimensional full-vectorial beam propagation method,” IEEE Photonics Technol. Lett. 17(11), 2319–2321 (2005). [CrossRef]

4. G. R. Hadley, “Wide-angle beam propagation using Padé approximant operators,” Opt. Lett. 17(20), 1426 (1992). [CrossRef]

5. C. Xu, W. Huang, J. Chrostowski, and S. Chaudhuri, “A full-vectorial beam propagation method for anisotropic waveguides,” J. Lightwave Technol. 12(11), 1926–1931 (1994). [CrossRef]

6. E. E. Kriezis and S. J. Elston, “Wide-angle beam propagation method for liquid-crystal device calculations,” Appl. Opt. 39(31), 5707 (2000). [CrossRef]

7. G. D. Ziogos and E. E. Kriezis, “Modeling light propagation in liquid crystal devices with a 3-D full-vector finite-element beam propagation method,” Opt. Quantum Electron. 40(10), 733–748 (2008). [CrossRef]

8. P. J. M. Vanbrabant, J. Beeckman, K. Neyts, R. James, and F. A. Fernandez, “A finite element beam propagation method for simulation of liquid crystal devices,” Opt. Express 17(13), 10895 (2009). [CrossRef]

9. G. Lifante Pedrola, Beam Propagation Method for Design of Optical Waveguide Devices (John Wiley & Sons, 2015).

10. S. Gauza, H. Wang, C.-H. Wen, S.-T. Wu, A. J. Seed, and R. Dabrowski, “High Birefringence Isothiocyanato Tolane Liquid Crystals,” Jpn. J. Appl. Phys. 42(Part 1, No. 6A), 3463–3466 (2003). [CrossRef]

11. J. Li, C.-H. Wen, S. Gauza, R. Lu, and S.-T. Wu, “Refractive Indices of Liquid Crystals for Display Applications,” J. Disp. Technol. 1(1), 51–61 (2005). [CrossRef]

12. M. Rubin, “Optical properties of soda lime silica glasses,” Sol. Energy Mater. 12(4), 275–288 (1985). [CrossRef]

13. Z. Gajic and M. T. J. Qureshi, Lyapunov matrix equation in system stability and control (Courier Corporation, 2008).

14. M. Schubert, C. Cramer, J. A. Woollam, C. M. Herzinger, B. Johs, H. Schmiedel, and B. Rheinländer, “Generalized transmission ellipsometry for twisted biaxial dielectric media: application to chiral liquid crystals,” J. Opt. Soc. Am. A 13(9), 1930 (1996). [CrossRef]

15. J. W. Goodman, Introduction to Fourier optics (Roberts and Company Publishers, 2005).

16. P. Del Moral and A. Niclas, “A Taylor expansion of the square root matrix function,” J. Math. Anal. Appl. 465(1), 259–266 (2018). [CrossRef]

17. I. Farago and A. Havasi, “On the convergence and local splitting error of different splitting schemes,” Prog. Comput. Fluid Dyn. 5(8), 495 (2005). [CrossRef]

18. G. R. Hadley, “Transparent boundary condition for beam propagation,” Opt. Lett. 16(9), 624–626 (1991). [CrossRef]

19. I. Nys, V. Nersesyan, J. Beeckman, and K. Neyts, “Complex liquid crystal superstructures induced by periodic photo-alignment at top and bottom substrates,” Soft Matter 14(33), 6892–6902 (2018). [CrossRef]

20. F. A. Milinazzo, C. A. Zala, and G. H. Brooke, “Rational square-root approximations for parabolic equation algorithms,” J. Acoust. Soc. Am. 101(2), 760–766 (1997). [CrossRef]

21. G. R. Hadley, “Multistep method for wide-angle beam propagation,” Opt. Lett. 17(24), 1743 (1992). [CrossRef]

22. A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, “Meep: A flexible free-software package for electromagnetic simulations by the fdtd method,” Comput. Phys. Commun. 181(3), 687–702 (2010). [CrossRef]

23. M. A. Karpierz, “Solitary waves in liquid crystalline waveguides,” Phys. Rev. E 66(3), 036603 (2002). [CrossRef]

24. K. Saitoh and M. Koshiba, “Full-vectorial finite element beam propagation method with perfectly matched layers for anisotropic optical waveguides,” J. Lightwave Technol. 19(3), 405–413 (2001). [CrossRef]

25. P. Oswald and P. Pieranski, Nematic and cholesteric liquid crystals: concepts and physical properties illustrated by experiments (CRC, 2006).

26. G. Poy, F. Bunel, and P. Oswald, “Role of anchoring energy on the texture of cholesteric droplets: Finite-element simulations and experiments,” Phys. Rev. E 96(1), 012705 (2017). [CrossRef]

27. H. L. Ong, M. Schadt, and I. F. Chang, “Material Parameters and Intrinsic Optical Bistability in Room Temperature Nematics RO-TN-200, -201, -403, E7, m₁, and m₃,” Mol. Cryst. Liq. Cryst. 132(1-2), 45–52 (1986). [CrossRef]

28. P. D. Brimicombe, C. Kischka, S. J. Elston, and E. P. Raynes, “Measurement of the twist elastic constant of nematic liquid crystals using pi-cell devices,” J. Appl. Phys. 101(4), 043108 (2007). [CrossRef]

29. H. N. Ritland, “Relation Between Refractive Index and Density of a Glass at Constant Temperature,” J. Am. Ceram. Soc. 38(2), 86–88 (1955). [CrossRef]

30. J. Hu and C. R. Menyuk, “Understanding leaky modes: slab waveguide revisited,” Adv. Opt. Photonics 1(1), 58 (2009). [CrossRef]

31. C. Hu and J. R. Whinnery, “Losses of a nematic liquid-crystal optical waveguide*,” J. Opt. Soc. Am. 64(11), 1424 (1974). [CrossRef]

32. P. J. Ackerman, Z. Qi, Y. Lin, C. W. Twombly, M. J. Laviada, Y. Lansac, and I. I. Smalyukh, “Laser-directed hierarchical assembly of liquid crystal defects and control of optical phase singularities,” Sci. Rep. 2(1), 414 (2012). [CrossRef]

33. G. Poy, A. J. Hess, I. I. Smalyukh, and S. Žumer, “Chirality-enhanced periodic self-focusing of light in soft birefringent media,” Phys. Rev. Lett. (to be published).

Physics-based multistep beam propagation in inhomogeneous birefringent media

Abstract

1. Introduction

2. Beam propagation methods and continuity equations

2.1 General propagation equations and approximations

2.2 Paraxial scheme

2.3 Wide-angle scheme

2.4 Interface conditions for the optical fields

3. Numerical results and discussion

3.1 Validation of the scheme

3.2 Waveguides based on topological solitons

4. Conclusion

Funding

Acknowledgements

Disclosures

References

Supplementary Material (1)

Cited By

Figures (6)

Equations (27)

Optics Express