Temporal dynamics of light-written waveguides in unbiased liquid crystals

TL;DR

This paper investigates the temporal formation dynamics of light-written waveguides in nematic liquid crystals, revealing that higher powers cause fluctuations and longer formation times due to self-steering effects.

Contribution

It provides a theoretical analysis of the temporal dynamics of light-written waveguides in liquid crystals, highlighting the effects of spatial walk-off and nonlocality on formation time.

Findings

Higher power induces beam fluctuations before stabilization.

Longer formation times occur at increased power levels.

Spatial walk-off influences waveguide formation dynamics.

Abstract

The control of light by light is one of the main aims in modern photonics. In this context, a fundamental cornerstone is the realization of light-written waveguides in real time, resulting in all-optical reconfigurability of communication networks. Light-written waveguides are often associated with spatial solitons, that is, non-diffracting waves due to a nonlinear self-focusing effect in the harmonic regime. From an applicative point of view, it is important to establish the temporal dynamics for the formation of such light-written guides. Here we investigate theoretically the temporal dynamics in nematic liquid crystals, a material where spatial solitons can be induced using continuous wave (CW) lasers with few milliWatts power. We fully address the role of the spatial walk-off and the longitudinal nonlocality in the waveguide formation. We show that, for powers large enough to induce light self-steering, the beam undergoes several fluctuations before reaching the stationary regime, in turn leading to a much longer formation time for the light-written waveguide.

Equations

$$2i{k}_0{n}_e({\theta }_b)\left(\frac{\partial A}{\partial z}+\tan \delta ({\theta }_b)\frac{\partial A}{\partial y}\right)+\frac{{\partial }^{2}A}{\partial x^2}$$

$$+D_y\frac{{\partial }^{2}A}{\partial y^2}+k_0^2\Delta n_e^2(\theta )A=0.$$

$$\nu \frac{\partial \theta }{\partial t}=K{\nabla }^2\theta +\frac{{ϵ}_0{ϵ}_a}{4}\sin \left[2\left(\theta -{\delta }_b\right)\right]\left(|{\mathbf{E}}_t{|}^2-|{\mathbf{E}}_s{|}^2\right)$$

$$+\frac{{ϵ}_0{ϵ}_a}{2}\cos \left[2\left(\theta -{\delta }_b\right)\right]\text{Re}\left(E_t{E}_s^{\ast }\right),$$

$$\theta (t)={\theta }_0+\frac{{ϵ}_0{ϵ}_a}{4\nu }\sin \left[2({\theta }_0-{\delta }_0)\right]|{\mathbf{E}}_t{|}^2(t-t_0),$$

$$\Delta t=\frac{\nu }{K}\frac{|{\mathbf{E}}_t{|}^2}{\left|{\nabla }^2|{\mathbf{E}}_t{|}^2\right|}.$$

$$\min (\Delta t)\approx \frac{\nu }{K}\frac{w^2}{80}.$$

$$\frac{\partial f(\mathbf{r},t)}{\partial t}=D{\nabla }^{2}f(\mathbf{r},t)+F(\mathbf{r},t),$$

$$\tilde{f}(\mathbf{k},t)=\tilde{f}(\mathbf{k},0)e^{-{\mathbf{k}}^{2}Dt}+{\int }_0^t{e}^{-{\mathbf{k}}^{2}D(t-t^{\prime })}\tilde{F}(\mathbf{k},t^{\prime })d{t}^{\prime }.$$

$$\tilde{f}(\mathbf{k},t)=\tilde{f}(\mathbf{k},0)e^{-{\mathbf{k}}^{2}Dt}+\frac{\tilde{F}(\mathbf{k})}{{\mathbf{k}}^{2}D}\left(1-e^{-{\mathbf{k}}^{2}Dt}\right).$$

$${\tilde{f}}_{stationary}(\mathbf{k})=\frac{\tilde{F}(\mathbf{k})}{{\mathbf{k}}^{2}D},$$

$$\tilde{f}(\mathbf{k},t)=\left(1-{\mathbf{k}}^{2}Dt\right)\tilde{f}(\mathbf{k},0)+\tilde{F}(\mathbf{k})t,$$

$$f(\mathbf{r},t)=f(\mathbf{r},0)+\left(D{\nabla }^{2}f(\mathbf{r},0)+F(\mathbf{r})\right)t.$$

$$2i{k}_0{n}_e({\theta }_b)\left(\frac{\partial A}{\partial z}+\tan \delta ({\theta }_b)\frac{\partial A}{\partial y}\right)+D_y\frac{{\partial }^{2}A}{\partial y^2}$$

$$+k_0^2\Delta n_e^2(\theta )A=0.$$

$$\nu \frac{\partial \psi }{\partial t}$$

$$+\frac{{ϵ}_0{ϵ}_a}{4}\sin \left[2\left({\theta }_0+\psi -{\delta }_b\right)\right]\left(|{\mathbf{E}}_t{|}^2-|{\mathbf{E}}_s{|}^2\right)$$

$$+\frac{{ϵ}_0{ϵ}_a}{2}\cos \left[2\left({\theta }_0+\psi -{\delta }_b\right)\right]\text{Re}\left({\mathbf{E}}_t{\mathbf{E}}_s^{\ast }\right),$$

Full text

I Introduction

Nowadays, in optical communication systems light beams carry the information via amplitude or phase modulations originating from an electric signal, either a voltage or a current. The signal processing is usually made by electronic circuits: the consequent transduction from optical to electrical signal, back and forth, represents today the most important limitation on the maximum achievable speed in optical communications. Accordingly, a big effort has been dedicated to realize all-optical systems, where direct optical signal processing and light-controlling-light schemes can be realized Willner et al. (2014). All-optical signal processing typically (see Refs. Zhang et al. (2012); Zhao et al. (2016) for counterexamples) implies a nonlinear response in the material, that is, changes in the optical properties of the medium due to light Willner et al. (2014).

Although in some cases freely propagating beams are used to transmit information Chan (2006), typically the optical beams are localized in a specific spatial region using waveguides, the most famous being fiber optics Yariv (1997). Optical waveguides allow to minimize the footprint of photonic devices, one of the most important advantages of integrated optics with respect to free space optics. Starting from the concepts of all-optical signal processing and integrated optics, the next breakthrough is to realize guiding structures defined by light in real time, paving the way to all-optical control of the geometry and topology of optical networks. The natural candidates for this are spatial solitons (SSs), waves subject to spatial localization owing to the action of nonlinear self-focusing Bjorkholm and Ashkin (1974); Stegeman and Segev (1999); Kivshar and Agrawal (2003). In fact, spatial solitons are usually associated with a light-written waveguide inside the material Kelley (1965); Snyder et al. (1991); Duree et al. (1993); Rotschild et al. (2005); Man et al. (2013); Kelly et al. (2016); Bezryadina et al. (2017); Salazar-Romero et al. (2017), capable to guide as well as probe beams at other wavelengths and low power.

One of the most used material for the generation of stable two-dimensional SSs is nematic liquid crystals (NLCs) Peccianti and Assanto (2012); Assanto (2012). In NLCs several different mechanisms can lead to nonlinear optical effects, according to the NLC composition, the optical excitation(s) and the environmental conditions Simoni (1997). For CW inputs, the strongest contribution usually comes from reorientational nonlinearity, consisting in a collective rotation of the molecules when subject to an optical torque Tabiryan and Zeldovich (1980). Reorientational nonlinearity in NLC is very large, allowing the observation of nonlinear effects with input powers of few milliWatts. On the other side, the time response is quite slow, ranging from milliseconds to seconds, according to the geometry and size of the NLC cell Simoni (1997). The effect of self-focusing induced by reorientational nonlinearity on long propagation distances (i.e., beyond the plane wave approximation) has been first studied by Braun in 1993 Braun et al. (1993). In the latter experiment, the presence of the Fréedericksz transition (FT) inhibited the observation of a stable SS. Stable SSs in NLCs have been found in 2000 applying a bias to the NLC cell in order to change the initial direction of the NLC molecules and avoid the appearance of the FT Peccianti et al. (2000). Later on, such approach has been generalized to unbiased cells by changing the anchoring conditions of the NLC molecules at the edge of the sample Peccianti et al. (2004a). Unlike Townes soliton (i.e., SS in a local Kerr media Chiao et al. (1964)), SS in NLCs, often called nematicons Peccianti and Assanto (2012), are stable due to the spatially nonlocal nonlinear response Suter and Blasberg (1993); Bang et al. (2002); Conti et al. (2003) stemming from the strong elastic interaction between adjacent molecules Peccianti and Assanto (2012).

Here we are interested in the investigation of the temporal response in the formation of nematicons in unbiased NLC cells Peccianti et al. (2004a). From an applicative point of view, the knowledge of the dynamical behavior of nematicons is fundamental to understand the maximum speed achievable in the modification of the guiding structure. The dynamics of nematicon formation has been already studied in Ref. Beeckman et al. (2005a); Strinić et al. (2005) in the case of a biased cell. The absence of bias makes the self-written waveguides qualitatively different. As a matter of fact, the static (or low-frequency Peccianti et al. (2000)) electric field creates a transverse gradient in the refractive index distribution, the maximum being at the center of the cell. The bias-induced inhomogeneity works both as a strongly multimodal waveguide Beeckman et al. (2010) and a potential landscape where the soliton starts to oscillate Peccianti et al. (2004b); Beeckman et al. (2005b). Such dynamics is strongly affected by spatial walk-off Peccianti et al. (2004a) as well, a fundamental ingredient in the propagation of nematicons Braun et al. (1993); Alberucci et al. (2010) neglected in Ref. Beeckman et al. (2005a).

II Geometry and model

We consider a planar cell of thickness $L_{x}=100~{}\mu$m, length $L_{z}$ and infinitely extended along $y$. Two glass slides, parallel to the plane $yz$, provide the NLC confinement. The two slides are treated in order to provide a uniform alignment of the NLC molecules in the absence of external excitations, the direction in the bulk being at $\theta_{0}=45^{\circ}$ with respect to $\hat{z}$. We assume strong anchoring conditions, that is, the molecular alignment in proximity of the NLC/glass interface is not affected by the impinging optical field Simoni (1997); Peccianti and Assanto (2012). A fundamental Gaussian beam at $\lambda=$1064 nm is injected into the sample with wavevector $\bm{k}=k_{0}\hat{z}=\frac{2\pi}{\lambda}\hat{z}$ parallel to $\hat{z}$. Two other very thin glass slides, inducing a molecular alignment parallel to $\hat{y}$, are placed at the input ($z=0$) and the output ($z=L_{z}$) interfaces to avoid the meniscus that forms at NLC/air interface. The interfaces induce a perturbation on the molecular distribution on distances $\approx L_{x}$, which we will neglect in this work. Noteworthy, in the paraxial limit the perturbation close to the input interface does affect neither the input polarization nor the input wavevector, the main effect being a slight transverse shift owed to the $z-$dependent walk-off angle.

The mean molecular direction is provided by a vector field, called the director $\hat{n}$, and represented by the long axis of the green ellipses in Fig. 1. Here we are interested in reorientational nonlinearity, that is, rotation of the director induced by the optical field Simoni (1997); Assanto (2012). Thus, we will focus on $y$-polarized inputs ($x$-polarization induces reorientation only above a power threshold Durbin et al. (1981)). Due to the NLC anisotropy, the beam inside the cell propagates (i.e., the Poynting vector $\bm{S}$) forming a walk-off angle $\delta$ with $\hat{z}$. Hereafter we will also neglect effects associated with changes in temperature (in this paper we focus our attention on undoped NLC away from the isotropic-nematic transition; we also assume that the saturation of the reorientational nonlinearity is not achieved) Assanto (2012); Alberucci et al. (2017), considering the latter constant in time and space. Due to the chosen input polarization, the director is constrained to lie on the plane $yz$, thus the angle $\theta$ formed with respect to $\hat{z}$ describes univocally $\hat{n}$. Optically, the system behaves like an inhomogeneous positive uniaxial material, with refractive indices equal to $n_{\bot}$ and $n_{\|}$ for electric fields oscillating normal and parallel to $\hat{n}$, respectively. The extraordinary waves experience a refractive index $n_{e}(\theta)=\left(\cos^{2}\theta/n^{2}_{\bot}+\sin^{2}\theta/n^{2}_{\|}\right)^{-1/2}$. We also define the optical anisotropy $\epsilon_{a}=n_{\|}^{2}-n_{\bot}^{2}$.

After setting $H_{x}=Ae^{ik_{0}n_{e}(\theta_{b})z}$, the nonlinear light propagation in the paraxial approximation ($\partial^{2}_{z}A=0$) and in the harmonic regime (electromagnetic field $\propto e^{-i\omega t}$) can be described by Alberucci et al. (2010)

[TABLE]

In Eq. (1) we introduced the nonlinear index well $\Delta n_{e}^{2}(\theta)=n_{e}^{2}(\theta)-n_{e}^{2}(\theta_{b})$ and the diffraction coefficient $D_{y}=n_{e}^{2}(\theta_{b})/\left(\epsilon_{\bot}+\epsilon_{a}\cos^{2}\theta_{b}\right)$, where $\theta_{b}(z)$ is the average angle perceived by the beam on each section $z=constant$. Equation (1) accounts for nonperturbative nonlinear effects, including self-steering and nonlinearity saturation Alberucci and Assanto (2010).

In writing Eq. (1) we assumed that the medium is stationary, that is, the reorientation angle $\theta$ does not depend on time $t$. In our case this hypothesis clearly does not hold valid. In fact, the reorientation of the NLC molecules, under the single elastic constant approximation and neglecting inertial effects in the NLC dynamics Au et al. (1991), is governed by

[TABLE]

where we defined $\delta_{b}=\delta(\theta_{b})$, $K$ is the Frank elastic constant, and $\nu$ is the rotational viscosity. The subscripts $t$ and $s$ indicate the component of the electric field parallel to $\hat{t}(z)=\hat{y}\cos\delta_{b}(z)-\hat{z}\sin\delta_{b}(z)$ (direction of the electric field for the plane wave with $\bm{k}\|\hat{z}$) and $\hat{s}(z)=\hat{y}\sin\delta_{b}(z)+\hat{z}\cos\delta_{b}(z)$ (direction of the Poynting vector $\bm{S}$), respectively. Equation (2) thus provides a time-dependent $\theta$, seemingly in disagreement with the assumptions behind Eq. (1). Nonetheless, the time response of the optical field to material changes is much faster than typical NLC response times, thus we can safely suppose a quasi-stationary regime for the electromagnetic field, since the latter adapts instantaneously to the refractive index landscape.

Equation (2) must be solved jointly with the boundary conditions $\theta=\theta_{0}$ at all the edges (i.e., $x=0$, $x=L_{x}$, $|y|\rightarrow\infty$, $z=0$ and $z=L_{z}$) of the NLC layer. For the sake of simplicity, we introduce the nonlinear perturbation on the director angle $\psi=\theta-\theta_{0}$, corresponding to a vanishing condition at the boundaries.

Let us first discuss qualitatively solutions of Eq. (2) for a fixed optical field with a spatial extension smaller than $L_{x}$ (illumination starting at $t=t_{0}$) with $\bm{E}_{s}=0$. In the stationary regime (i.e., $t\rightarrow\infty$), it is well known that Eq. (2) provides a nonlinear perturbation extending on an area of transverse size comparable with $L_{x}$, regardless of the input beam width $w_{in}$ Khoo et al. (1987); Alberucci and Assanto (2007). If $w_{in}\ll L_{x}$, nonlinear light propagation is in the highly nonlocal regime. However, during the transient, the behavior of $\psi$ is different. Specifically, for short time intervals after the illumination, in Eq. (2) the elastic forces can be neglected, providing (see also Eq. (10) in the Appendix)

[TABLE]

where $\delta_{0}=\delta(\theta_{0})$ and where we supposed small perturbations, that is, $\psi\ll\theta_{0}$. According to Eq. (3), in this regime the NLC behaves spatially like a local material Henninot et al. (2002); Warenghem et al. (2006), whereas it shows a highly nonlocal response in the time domain Conti et al. (2010). Physically, the optically-induced perturbation has not enough time to spread away from the optical beam due to the slow response time of the intermolecular elastic forces Beeckman et al. (2005a); Henninot et al. (2007). Later on, the perturbation recovers its nonlocal nature, thus providing a stable optical soliton.

To evaluate the interval of validity $\Delta t$ for Eq. (3), we can substitute Eq. (3) into the RHS of Eq. (2) and find when the optical and the elastic torque are equal to each other in absolute value. We find

[TABLE]

According to Eq. (4), $\Delta t$ depends on the NLC material (parameter $\nu/K$, here the diffusion coefficient) and on the shape of the optical illumination, but not on the incident optical power. In particular, the larger the spatial derivatives of the intensity distribution the shorter is the time required for the appearance of spatially nonlocal effects.

III Solution in a longitudinally-invariant geometry

To verify the validity of Eq. (4), we consider the simplest case of a structure invariant along $z$ Au et al. (1991); Veltri et al. (2007). Hereafter, we will consider the NLC mixture E7 featuring $\nu=0.2$ Pa s and $K=12$ pN. At $\lambda=1064~{}$nm the E7 refractive indices are approximately $n_{\|}=1.7$ and $n_{\bot}=1.5$. For the optical illumination we make the ansatz $|\bm{E}_{t}|^{2}=\frac{4Z_{0}P}{\pi n_{e}w^{2}}\exp{(-2r^{2}/w^{2})}$ (we set $r=\sqrt{x^{2}+y^{2}}$), yielding $\Delta t=\nu w^{2}\left/\left(16K\left|\frac{r^{2}}{w^{2}}-\frac{1}{2}\right|\right)\right.$. Given that the transition time changes across the beam cross-section, we need to consider the minimum of $\Delta t$ on the transverse plane $xy$, the latter being achieved when $r\rightarrow\infty$, with the constraint that the field should not be negligibly small. Taking conventionally the point $r=2.35w$ (such a choice comes from the direct comparison with exact numerical simulations, see below), we find

[TABLE]

Thus, Eq. (3) holds valid on the temporal interval $[t_{0},\ t_{0}+T]$, where we set $T=\min(\Delta t)$.

Thus, according to Eq. (5), the narrower the beam the shorter is the validity range of Eq. (3). For example, for $w=2~{}\mu$m Eq. (5) provides $\min(\Delta t)\approx$ 1 ms.

Figure 2(a-b) compare the all-optical perturbation $\psi(x,y,t)$ given by Eq. (3) (dashed lines) and exact solutions (solid lines) computed via numerical simulations of Eq. (2). For quantitative comparison, we take the maximum of $\psi(x,y,t)$ at a fixed instant, let us call it $\psi_{m}(t)$. As predicted, for short times values for $\psi_{m}$ provided by the exact solution of Eq. (2) and approximate solution from Eq. (3) perfectly overlap [Fig. 2(c)], on a temporal interval $T$ in good agreement with Eq. (5) [Fig. 2(d)]. Numerical simulations plotted in Fig. 2(a-b) show how the transition from the local to the nonlocal regime takes place. Discrepancies first arise in correspondence to the peak and the tails of the field, where $|\nabla^{2}|\bm{E}_{t}|^{2}|$ is larger. The net effect is smoothing out the peak of the perturbation, with tails extending in time until reaching the closest boundary Alberucci and Assanto (2007). Numerical simulations confirm that the temporal dynamics of the maximum reorientation angle $\psi_{m}$ (i.e., $\psi_{m}(t)/\psi_{\infty}$, where $\psi_{\infty}=\lim_{t\rightarrow\infty}\psi_{m}(t)$ is the reorientation in the stationary regime) is independent from the input power as long as the small perturbation condition $\psi\ll\theta_{0}$ holds true. For example, for $\theta_{0}=\pi/4$ and $w=6~{}\mu$m, up to $P=5~{}$mW no appreciable variations occur in the curves $\psi_{m}(t)/\psi_{\infty}$ for different powers, whereas at $P=10~{}$mW a maximum discrepancy of 2.8$\%$ is observed.

Above (see results plotted in Fig. 2) we considered a fixed optical excitation. The next step is to account for the influence of the molecular rotation on light propagation Bloisi et al. (1988), accounted by Eq. (1), but conserving the $z$-invariance of our system, similarly to the approach employed in Ref. Alberucci et al. (2014). We thus consider the average value of the beam across $z$, neglecting diffraction or soliton breathing Alberucci et al. (2015). Specifically, we approximate the nonlinear index well with a cylindrically-symmetric parabolic index well, in turn providing a Gaussian profile $I=\frac{2P}{\pi w^{2}(t)}\exp{\left[-\frac{2r^{2}}{w^{2}(t)}\right]}$ for the beam, where $w(t)$ is now time-dependent due to the self-focusing. We thus make the approximation $\Delta n_{e}^{2}\approx 2n_{e}(\theta_{b})\left.\frac{dn_{e}}{d\theta}\right|_{\theta_{b}}\psi_{2}r^{2}$. The coefficient $\psi_{2}$ is computed by a best-fitting procedure for $\psi$ in the interval $x\in[-\sqrt{2}w(t),\ \sqrt{2}w(t)]$ Ouyang et al. (2006).

To calculate the temporal dynamics of the strongly coupled light-matter system, we use the following method. At the initial time $t=0$, we consider an impinging Gaussian beam of waist $w_{0}=w(t=0)$. We then compute one temporal step of Eq. (2). From the knowledge of the new rotation angle $\theta$, we can compute $\psi_{2}$, the latter being necessary to compute the new beam width. This basic procedure is iterated on the temporal interval of interest. Hence, at every step we are able to account for the effect of self-focusing on the intensity distribution. Finally, in the numerical algorithm we set the constraint $w(t)\leq w_{0}$ so that the beam width cannot be larger than the initial value. Figure 3 shows the results for $P=5~{}$mW and $w_{0}=6~{}\mu$m. In agreement with Eq. (4), the optical perturbation $\psi$ is local at the beginning (few tens of milliseconds), whereas diffusion comes into play at later times. After the stationary regime is achieved, the profile of the reorientation angle is dictated by the closest boundary (see Section V.2) Alberucci et al. (2010); Minovich et al. (2007). On the other hand, the beam profile in the first milliseconds is not changing because the predicted soliton would be wider than the input beam.

Figure 4 shows in more detail the dynamics for different input powers. When the self-focusing kicks in, the beam starts to shrink towards the stationary soliton solution (for $P=5~{}$mW it is $w_{s}\approx 2.8~{}\mu$m). Noteworthy, the beam reaches its final size after about 1 s [Fig. 4(b)], whereas the maximum and the width of the reorientation angle $\psi$ keep increasing on a longer temporal interval [Fig. 4(a), Fig. 5(a)].

For $P=1~{}$mW, the soliton width is larger than $w_{0}$, the latter meaning no changes in the optical beam according to our simplified model. Thus, the reorientation follows the case for fixed beam previously investigated (see Fig. 2). For $P=5~{}$mW, self-focusing takes place given that the soliton width in the stationary regime is $2.8~{}\mu$m. Beam narrowing starts after 42 ms, leading to a fast decrease in the beam width. After about 1 s, the nonlinear confinement parameter $\psi_{2}$ (thus, the beam width $w$) reaches the stationary value. Despite that, the maximum rotation angle keeps increasing (achieving $95\%$ of the stationary value at $t\approx 20~{}$s), but without appreciable effects on the self-trapped beam. Reorientation in time is slightly steeper for $P=5~{}$mW than for $P=1~{}$mW due to the narrower beam, leading to a faster reorientation [see Eq. (8)]. When $P=20~{}$mW, both reorientation and beam focusing dynamics are steeper than for $P=5~{}$mW (now the $95\%$ of the stationary value is achieved at $t\approx 12~{}$s). Self-focusing now starts after 8 ms. Due to the fact that $\psi$ is now comparable with $\theta_{0}$ (i.e., we are in the nonperturbative regime), the shape of the normalized reorientation curve $\psi_{m}/\psi_{\infty}$ now strongly differs from the two previous cases. The nonperturbative behavior of the nonlinearity is responsible for a non-monotonic behavior in time for $\psi_{2}$ and $w$, as well. Summarizing, the self-focusing accelerates the reorientation dynamics owing to a stronger external torque associated with a narrower optical beam, in agreement with Eq. (4).

In Fig. 5 we quantitatively analyze the transition from a local to a nonlocal response with time. Figure 5(a) shows the width of the optically-induced perturbation $\psi$, defined through the second central moment as $w_{\psi}=2\sqrt{\int\psi\eta^{2}d\eta/\int\psi d\eta}\ (\eta=x,y)$ with respect to time. The perturbation width is monotonically increasing both for $P=1$ and 5 mW, whereas for 20 mW an initial shrinking is observed due to the stronger amount of self-focusing. In the stationary regime $w_{\psi}$ is the same for input powers larger than few mWs (the shape of the solution of Eq. (2) is almost independent on the soliton width $w$), whereas at $1~{}$mW the size of the beam leads to a wider reorientation in space. Although the asymmetric boundary conditions, beam width $w_{\psi}$ is almost the same along $x$ and $y$, i.e., the nonlinear index well is cylindrically symmetric within a good accuracy. Figure 5(b) shows the temporal evolution of nonlocality $\sigma$ defined as the ratio between the widths of the optical perturbation and of the intensity profile, $\sigma=w_{\psi}/w$. At the initial time $t=0$, the medium is local and $\sigma=1$, in agreement with Eq. (3). Nonlocality then increases monotonically in time due to the diffusion effect in agreement with Eq. (8), the larger the power the steeper the increase is.

IV Formation of the waveguides: the longitudinal dynamics

In Section III we described the temporal evolution of light self-trapping neglecting the dynamics of the system along $z$. Such simplified model, despite providing a qualitative picture of the dynamics, does not account for relevant physical effects, such as self-steering Alberucci et al. (2010), finite boundary conditions along $z$ McLaughlin et al. (1995), oscillation in the beam width Conti et al. (2004). The exact 3D numerical solution of Eqs. (1-2) is computationally highly demanding, thus we will account for the longitudinal evolution of the system by using a simplified (1+1)D model (see Section V.3). In the simplified model we define an effective power $P_{2D}$, which is about 3$\sim$10 times smaller than the real power $P$ Alberucci et al. (2010). In this section all the simulations are carried out considering an input Gaussian beam with planar phase front in $z=0$ and waist $w_{\mathrm{in}}=2~{}\mu$m. Finally, the cell length along the propagation distance $z$ is set to $1$ mm.

Snapshots of typical light distribution on the plane $yz$ are shown in Fig. 6. Videos of the temporal dynamics for $P_{2D}=0.5,~{}1$ and 3 mW on short (time window of 1 s) and long scales (tens of seconds) are reported in Visualization 1-3 and 4-6, respectively. For any input power, the amount of light self-localization grows with time. The effect of nonlinearity on the beam width $w(z,t)$ is plotted in Fig. 7. At low power ($P_{2D}=0.1~{}$mW), the net effect of nonlinearity is a decrease in the spreading angle. At $P_{2D}=0.5~{}$mW, the beam undergoes a strong focusing effect in the first 2 s, whereas for longer times slight adjustments towards the stationary condition occur. The final condition is a breather, that is, a self-localized mode subject to quasi-periodic oscillations during propagation Conti et al. (2004); Kivshar and Agrawal (2003). At $P_{2D}=1~{}$mW, spatial self-localization occurs faster, but the overall behavior is very close to the case $P_{2D}=0.5~{}$mW. In all of these cases, the beam width undergoes a continuous transition from a linear spreading (diffraction) to a localized breather. Light propagation for $P_{2D}=3~{}$mW is more complicated. In fact, in this case the optical perturbation $\psi$ becomes comparable with $\theta_{0}$, thus light is able to change its own path by changing the walk-off angle Piccardi et al. (2010). The effect is already visible at $t=1~{}$s, where the beam is propagating along a curved trajectory (see Fig. 6). A longer temporal interval is necessary before reaching the stationary regime: after $20~{}$s the light trajectory is not a straight line yet, as it should be in the stationary regime Alberucci et al. (2010). As a matter of fact, at the beginning the beam oscillates in a seemingly random fashion in time. Then the peaks and valleys on the beam width computed in the stationary regime emerge one by one from the beam wiggling, starting from the input interface up to the output section (last panel in Fig. 7).

From an engineering point of view, we are particularly interested in the minimum time required to form the waveguide. To address this point, the short-time dynamics of the beam width is plotted in Fig. 8 (see also Visualization 1-3). Regardless of the input power, self-confinement takes place first close to the input interface, eventually moving to larger $z$ with time. The time required for the formation of waveguides strongly depends on the input power. We will consider the waveguide formed when the beam width at the output ($z=1$ mm) is lower than 10$~{}\mu$m. When $P_{2D}=0.5$ mW [Fig. 8(a)], the waveguide is not yet formed in the first two seconds, whereas for $P_{2D}=1$ mW waveguiding is achieved after 1 s [Fig. 8(b)]. For $P_{2D}>3~{}$mW, localization takes place at $t=0.1~{}$s [Fig. 8(c-d)]. Important to stress, the guides formed in this transient regime differ from the stationary value, the difference increasing with the input power. In fact, in Fig. 8(c-d) we notice that for $t>0.5~{}s$ the guide at large $z$ is destroyed due to trajectory oscillations associated with self-steering Braun et al. (1993).

To summarize, the waveguides can be formed in a temporal interval of hundreds of milliseconds, however, using powers large enough to induce oscillations in the beam path on longer times. The drawback is a longer relaxation time before the stationary regime is reached. Let us then discuss the nature of the oscillations in the beam trajectory (see Visualization 4-6). Figure 9 shows the behavior of the light path, defined as the center-of-mass of the beam. When nonlinear changes in the walk-off are negligible [Fig. 9(a-b)], the trajectories slowly relax towards the stationary condition. When self-steering is appreciable [Fig. 9(c-d)], at very short times the beam propagates along straight lines determined by the linear value of the walk-off, $\delta(\theta_{0})$. As time goes on, the reorientation angle becomes large enough to misplace the beam with respect to the induced index well, the latter being inhibited to follow instantaneously the changes in the light distribution owing to the slow response of the material. The net effect is that the light beam starts to oscillate around the index well created at the previous time by the light itself. At later times and starting from the input interface, the beam slowly relaxes towards the final trajectory, corresponding to a straight line with a lower slope than the initial one (walk-off $\delta$ encompasses a maximum in proximity of $\theta=45^{\circ}$).

Next step is to analyze the features of the waveguide and their evolution in time. Figure 10 reports the maximum director angle $\theta$ computed on each plane $z=\text{const}$ (dubbed $\theta_{m}(z,t)$) and the width of the waveguide $w_{\psi}(z,t)$. The magnitude of $\theta_{m}$ versus time closely resembles  trends shown in Fig. 4, with larger powers corresponding to faster reorientation. For $P_{2D}=0.5~{}$mW, $\theta_{m}$ decays with $z$ owing to the beam spreading in propagation, see the second panel in Fig. 7. For the same reason, the width of the guide increases with $z$. For larger input powers the maximum reorientation angle is almost constant along $z$ due to the filtering action of the spatial nonlocality along $z$. When self-steering occurs ($P_{2D}=3$ mW) and in accordance with the oscillation in the beam path, $\theta$ varies rapidly close to the output interfaces, stabilizing after about 20 s. Regardless of the input power and in agreement with the simplified model discussed in Section III, nonlocality increases monotonically with time. Thus, the waveguides observed at short time intervals correspond actually to a different and narrower spatial profile with respect to the stationary case (see Section V.2).

V Discussion of the results and future perspectives

In this article we investigated the temporal dynamics associated with the formation of light-written waveguides in liquid crystals. We showed how the degree of spatial nonlocality evolves in time from a local (the waveguide shape is identical to the intensity profile) to a highly nonlocal regime (the shape and extension of the waveguide is dictated by the sample geometry), such result holding valid in thermal materials as well Rothschild et al. (2006); Kaminer et al. (2007); Navarrete et al. (2017). An interesting consequence of our results is the possibility to control the waveguide spatial extension varying the repetition rate of a train of pulses.

With respect to the reorientational nonlinearity in NLC, two main regimes have been found. When the input power is not large enough to induce self-steering (the perturbative regime), the nonlocal waveguide is formed in few seconds, and the formation time is shorter for increasing powers. Despite that, the full stationary regime is achieved in tens of seconds. The multi stage dynamics arises from the different time constants associated with every spatial mode of the nonlinear perturbation (the wider the mode the longer the response time is). Importantly, these results hold valid in thermal materials such as lead glasses Dabby and Whinnery (1968); Rotschild et al. (2005); Alfassi et al. (2007) or, more in general, in the presence of a nonlinear mechanism based upon diffusion (we stress that in NLC more complicated thermal effects can appear owing to convective heat exchange or owing to the presence of the nematic-isotropic phase transition Simoni (1997); Shih et al. (2010)). In the nonperturbative regime, the power is large enough to induce a power-dependent beam trajectory. The nonperturbative regime is typical of the reorientational nonlinearity in NLC, and it can be described only accounting for the power-dependent spatial walk-off and the longitudinal evolution of the self-written waveguide. In the nonperturbative regime, the beam starts to fluctuate during the transient regime owing to the misalignment between the beam and the waveguide. As a matter of fact, the stationary regime is reached after a longer time interval with respect to the perturbative regime. Analogously to the perturbative case, a local waveguide can be written after fractions of seconds.

Our results demonstrate that temporal modulation of the impinging beam represents an additional degree of freedom in the control of light-written waveguides in nonlocal nonlinear media Henninot et al. (2007); Vocke et al. (2017). In a broader perspective, our results represent the first step towards a full investigation of the interplay between temporal and spatial effects in highly nonlocal materials. With respect to the field of light self-localization in nematic liquid crystals, future developments include the assessment of the role played by the strong scattering associated with the collective molecular fluctuations, the latter being associated with a temporal instability for self-trapped beam at large enough powers Braun et al. (1993); Valkov et al. (1985); Bolis et al. (2017). Finally, in future works it will be interesting to pursue the generalizations of the phenomena discussed here to the case of smectic phase, where much faster time responses with respect to the nematic phase can be achieved, but at the expense of stronger elastic interactions Clark and Lagerwall (1980).

Funding. Deutsche Forschungsgemeinschaft (DFG) (GRK 2101); Academy of Finland, Finland Distinguished Professor grant No. 282858; ERC Adv. grant No. 694683 (PHOSPhOR); FCT grant No. SFRH/BPD/77524/2011

Appendix

V.1 Time-dependent solution of a forced diffusion equation

Let us consider the normalized equation

[TABLE]

where $F(\bm{r},t)$ is the forcing term, in the most general case dependent on time as well. Using the 3D spatial Fourier transform, $\tilde{f}(\bm{k},t)=(2\pi)^{-\frac{3}{2}}\iiint{f(\bm{r},t)e^{i\bm{k}\cdot\bm{r}}d^{3}\bm{r}}$, we find the solution in the transformed domain

[TABLE]

If $F$ is independent on time, we get

[TABLE]

In the stationary regime, $t\rightarrow\infty$

[TABLE]

corresponding to the solution of the Poisson equation. Boundary conditions, i.e., the shape of the sample, affect Eq. (9) by selecting the allowed spatial frequencies $\bm{k}$. The stationary solution provided by Eq. (9) does not depend on the initial condition $f(\bm{r},0)$. For short time intervals, Eq. (8) yields

[TABLE]

corresponding in the real domain to

[TABLE]

V.2 Shape-preserving solitons in the stationary regime

The fundamental shape-preserving solitons of Eqs. (1-2) in the stationary regime (i.e., $\partial_{t}\theta=0$) are found by setting $A=A_{s}(x,y-\tan\delta_{s}z)e^{i\beta_{s}z}$ and $\theta=\theta_{s}(x,y-\tan\delta_{s}z)$. The ansatz corresponds to a soliton with planar phase-front and propagating along a walk-off angle $\delta_{s}$, the latter being also dependent on the input power in the non-perturbative regime. Then Eqs. (1-2) turn into a nonlinear eigenvalue problem parametrized with respect to the input power $P$. The eigenfunction $A_{s}$ and the eigenvalue $\beta_{s}$ of the electromagnetic problem are found as the corresponding eigenvector and eigenvalue of the discrete operator corresponding to Eq. (1). The optical profile found is then substituted into Eq. (2) and the new director rotation is obtained. The new material parameters are then substituted back into Eq. (1), and the procedure is iterated until convergence is achieved.

Figure 11(a) and (b) show the beam width and the walk-off angle $\delta_{s}$ versus power $P$. The beam width is not monotonically shrinking owing to the saturation effect related to large reorientation angles, in turn weakening the effective nonlinearity. Figure 11(c) and (d) show respectively the functions $A_{s}$ and $\theta_{s}$ versus the transverse coordinates for $P=5.1~{}$mW. Shapes of the soliton changes slightly with power, conserving a quasi-Gaussian profile. The reorientation also maintains its transverse profile, the latter being determined mainly by the boundary conditions along $x$, i.e., the shortest size of the sample Hutsebaut et al. (2005); Minovich et al. (2007). Despite the asymmetric boundary conditions, the reorientation angle in proximity to the cell center is cylindrically symmetric.

V.3 Effective (1+1)D model

Numerical simulations of Eqs. (1-2) are highly demanding from the computational time point of view. An effective bidimensional problem can be solved in behalf of Eqs. (1-2). The model is Alberucci et al. (2010)

[TABLE]

with an effective power $P_{2D}$ lower than the real one. The screening term in Eq. (13) is inserted in order to conserve the nonlocality degree of the original system Kaminer et al. (2007). According to the results shown in Section V.1, the response time depends critically on the shape of the Green function for the material response. Thus, the simplified model conserves the response time of the original 3D model.