Helicoids in chiral liquid crystals under external fields

TL;DR

This paper investigates how cholesteric liquid crystals form various helicoidal structures under external magnetic fields and confinement, providing a phase diagram of transitions from nematic to helicoidal phases.

Contribution

It offers new analytical and numerical insights into the phase transitions and configurations of cholesteric liquid crystals under external fields and confinement.

Findings

Phase diagram of nematic to helicoidal transitions

Characterization of helicoidal configurations with disclinations

Analytical and numerical methods applied

Abstract

Cholesteric Liquid Crystals (CLCs), subject to externally applied magnetic fields and confined between two parallel planar surfaces with strong homeotropic anchoring conditions, are found to undergo transitions to different types of helicoidal configurations with disclinations. Analytical and numerical studies are performed in order to characterise their properties. In particular, we produce a phase diagram for the transitions from the nematic state to the helicoidal phases in terms of the molecular chirality and the strength of the applied magnetic field.

Equations

$$\mathcal{B}=\{(x,y,z)\in {\mathbb{R}}^3,\mid z\mid \le \frac{L}{2}\},$$

$$\mathbf{n}\left(\mathbf{r}\right)=\left(\sin \theta (\mathbf{r})\cos \varphi (\mathbf{r}),\sin \theta (\mathbf{r})\sin \varphi (\mathbf{r}),\cos \theta (\mathbf{r})\right),$$

$${\mathcal{E}}_{FO}=$$

$$\mathbf{n}\left(\mathbf{r}\right)=\left(0,-\sin \theta \left(x,z\right),\cos \theta \left(x,z\right)\right),$$

$$E_{FO-2d}$$

$${\partial }_x^2\Theta +{\partial }_z^2\Theta ={\Lambda }^2\sin \Theta ,\Theta =2\theta ,\Lambda =\sqrt{\frac{{\chi }_a}{K}}\text{H},$$

$$\Theta =4\arctan \left[X\left(x\right)Z\left(z\right)\right],$$

$$\frac{X^{\prime \prime }}{X{Z}^2}+X^2\left(\frac{Z^{\prime \prime }}{Z}-\frac{Z^{\prime 2}}{Z^2}\right)=X{X}^{\prime \prime }-2X^{\prime 2}+{\Lambda }^2{X}^2+\frac{Z^{\prime \prime }}{Z^3}-\frac{{\Lambda }^2}{Z^2},$$

$$2X{X}^{\prime }{\left[\frac{Z^{\prime \prime }}{Z}-\frac{Z^{\prime 2}}{Z^2}\right]}^{{}^{\prime }}-\frac{2Z^{\prime }}{Z^3}{\left[\frac{X^{\prime \prime }}{X}\right]}^{{}^{\prime }}=0.$$

$$\frac{1}{X^{\prime }X}{\left[\frac{X^{{}^{\prime \prime }}}{X}\right]}^{{}^{\prime }}=4a,$$

$$\frac{Z^3}{Z^{{}^{\prime }}}{\left[\frac{Z^{{}^{\prime \prime }}}{Z}-2{\left(\frac{Z^{{}^{\prime }}}{Z}\right)}^2\right]}^{{}^{\prime }}=4a,$$

$$X^{{}^{\prime \prime }}=2a{X}^3+bX,$$

$$\frac{Z^{{}^{\prime \prime }}}{Z}-2{\left(\frac{Z^{{}^{\prime }}}{Z}\right)}^2+\frac{2a}{Z^2}+d=0,$$

$$(X^{{}^{\prime }}{)}^2=a{X}^4+b{X}^2+c,$$

$$(Z^{{}^{\prime }}{)}^2=g{Z}^4+d{Z}^2+a.$$

$$g=c,b={\Lambda }^2-d.$$

$$(X^{{}^{\prime }}{)}^2=a{X}^4+({\Lambda }^2-d)X^2+c,$$

$$(Z^{{}^{\prime }}{)}^2=c{Z}^4+d{Z}^2+a,$$

$$X\left(x\right)=\pm \sqrt{\frac{{\Lambda }^2-d}{a}}\text{csch}\left(\sqrt{{\Lambda }^2-d}x\right),$$

$${\theta }_n$$

$$({\partial }_x-ı{\partial }_y)\left(\frac{{\Theta }_i-\Theta }{2}\right)$$

$$({\partial }_x+ı{\partial }_y)\left(\frac{{\Theta }_i+\Theta }{2}\right)$$

$$\bar{\Theta }=\Theta +4\arctan \left[\frac{{\beta }_2+{\beta }_1}{{\beta }_2-{\beta }_1}\tan \left(\frac{{\Theta }_2-{\Theta }_1}{4}\right)\right].$$

$${\theta }_0=-\text{sign}\left[\cos t\right]\pi$$

$$\Delta E_{2\pi }=E_{2\pi }-E_{\mathrm{nem}}\approx$$

$$\Delta E_{2\pi }=2KL(q_H-q_0)\left(\pi -\frac{1}{L}\sqrt{\frac{K}{{\mathcal{E}}_{max}}}\left(\pi +\log (4)\right)\right),$$

$$q_H$$

$$\theta\left(x,z=\pm\frac{L}{2}\right)=\left\{\begin{array}[]{cc}\pi&x<0\\ 0&x>0\end{array}\right..$$

$$W$$

$$\Omega \left(\lambda \right)=\frac{ı\Lambda }{2}\left(\frac{1}{\lambda }-\lambda \right),\omega \left(\lambda \right)=\frac{\Lambda }{2}\left(\frac{1}{\lambda }+\lambda \right),$$

$$Q\left(x,z,\lambda \right)=$$

Full text

I Introduction

In the absence of external applied magnetic or electric fields and in free space, a cholesteric liquid crystal twists uniformly about a single axis: rod-like chiral molecules self-assemble into a helical arrangement along a single direction. Under confinement, the natural twisted structures of the chiral nematic liquid crystals are often incommensurate with the geometry, dimension and surface anchoring of the confining regions. Due to this geometric frustration, the confined chiral liquid crystals struggle to maintain their helical pitch and the preferred directors at the confining walls. A similar frustration mechanism takes place for a cholesteric liquid crystal in an electric or magnetic field: the field alignment is incompatible with the cholesteric twist. In both cases of frustration, the cholesteric is subject to an anisotropic environment: the competition between favored twist and anisotropy leads to the distortion and the partial or total unwinding of the helical nematic texture and to the formation of intriguing frustrated director configurations with topological defects 1 ; 2 . The simultaneous interplay of the intrinsic molecular chirality, the external fields and the boundary conditions, generates new structures, either localised or extended. These structures can be elongated string-like objects called cholesteric fingers (or threads) or helicoids 4 ; baudry2 ; oswald1 ; oswald2 ; 3 ; 8 , or they may be localised objects called cholesteric bubbles or spherulitic domains tako ; chandra ; nawa ; coche . Recently, it has been recognized that cholesteric bubbles can have the remarkable properties of skyrmions. Skyrmions, originally proposed in the field of nuclear physics skyrme , are localised structures in which the magnitude of the order parameter (the nematic director $\mathbf{n}$ in this case) remains constant, but the orientation continuously varies in a complex texture that cannot be annealed away. Moreover, very recently, Smalyukh et al. generated in confined chiral liquid crystals a type of defects in the director field configurations, called triple-twist torons, by geometrical frustration and by using Laguerre-Gaussian vortex laser beams. They also showed numerically their existence from a theoretical point of view acke1 ; acke2 ; acke3 .

All such configurations, i.e. helicoids, skyrmions, torons, and other specific solitonic textures, e. g. hopfions and merons acke1 ; acke2 ; acke3 ; pandey , are stabilised by some type of topological and/or non-topological conservation laws 13 and, at least in some approximate setting, they can be described in terms of integrable nonlinear equations Spheruliti . Moreover, these novel structures have been studied through a range of techniques, including experiments and numerical simulations in several confinement geometries such as thin layers, wedge cells, cylindrical cavities and spherical droplets lansac ; lagemaat ; fukudazumer ; eun . Relevant variational and numerical calculations for skyrmions have been done in two dimensions by Bogdanov et al. bogdanov1 ; bogdanov2 , and more recently in three dimensions by Leonov et al. leonov2014 . In the latter paper authors studied the Frank-Oseen free energy in the confined geometry

[TABLE]

$L$ being the nematic cell gap, under the action of an external field and with weak anchoring conditions, i.e. the Rapini-Popoular anchoring surface energy stewart_book . They calculated the director texture for skyrmions. The boundary conditions translate into conditions on the partial derivatives of the director on the confining surfaces. Later, in 17 , the authors performed similar variational and numerical calculations in the same confining geometry but for the case of strong anchoring ($\mathbf{n}$ perpendicular at the confining surfaces) and in the absence of external fields. In particular the authors of 17 worked out analytical and numerical solutions for a single helicoid and for a helicoid lattice. By using similar calculations, they investigated isolated skyrmions and skyrmion lattices.

In this paper, we consider a chiral liquid crystal confined within $\mathcal{B}$ as in (1) with strong homeotropic anchoring conditions and in the presence of an external magnetic field.

As opposed to 17 , the inclusion of an external field leads to a nonlinear partial differential equation (PDE) for the director orientation angle, more precisely the elliptic sine-Gordon equation on the strip, possibly with discontinuous boundary conditions. New static chiral states are found and recognised to be similar to cholesteric fingers (helicoids) with defects of disclination type 4 ; 8 . We classify and describe these configurations by analytical and numerical methods. In particular, we discuss a new type of solutions called $2\pi-$helicoids, where the director field $\mathbf{n}(\mathbf{r})$, independent of $y$, rotates by $2\pi$ over the strip, that is, the projection of $\mathcal{B}$ onto the (x, z) plane. In addition to the $2\pi-$helicoids, we also find $\pi-$helicoids, where the director $\mathbf{n}(\mathbf{r})$ only twists once over the strip. All these solutions contain disclination-type singularities. Accordingly, the evaluation of the static free energy leads to the introduction of a phenomenological cut-off, which, in turn, determines a critical parameter for the transition to the nematic uniform state. Transitions from the uniform nematic state to $2\pi-$ and $\pi-$helicoids are represented in the phase space of the spontaneous chiral twist and magnetic strengths.

The paper is organised as follows. In Sec. II, we introduce the mathematical model. In Sec. III and IV we solve the nonlinear PDE problem and find the solutions in a closed analytical form. In Sec. V, we build up a phase diagram of the solutions by energy comparison. Finally, in Sec. VI, we draw our conclusions and plan future work.

II The model

We consider a static CLC layer, confined in between two identical planar surfaces placed at $z=\pm\frac{L}{2}$, and extending to infinity in the orthogonal directions $\left(x,y\right)$ in a suitable cartesian reference system $\left(O,x,y,z\right)$, with orthonormal basis vectors $\mathbf{x},\mathbf{y},\mathbf{z}$. The system is described by the uni-modular director field $\mathbf{n}\left(\mathbf{r}\right)\in\mathbb{RP}^{2}$

[TABLE]

due to the $Z_{2}$ symmetry of the microscopic model deGennes , and governed by the Frank-Oseen free energy density

[TABLE]

where $q_{0}$ is the spontaneous chirality constant of the cholesteric phase, while the positive reals $K_{i}$ denote the splay, twist and bend Frank elastic constants, respectively, for which we use the simplifying one constant approximation $K=K_{1}=K_{2}=K_{3}.$ The last term represents the interaction energy density with an external static magnetic field $\mathbf{H}=\text{H}\;\mathbf{z}$, which is assumed to be uniform.

At the bounding surfaces, we impose strong homeotropic anchoring conditions, i.e. $\mathbf{n}\left(x,y,z=\pm\frac{L}{2}\right)=\mathbf{z}.$

Both $\mathbf{H}$ and the confinement break the general rotational and translational symmetry along the $\mathbf{z}$ direction of the fundamental cholesteric helices, i.e. ${\bf n}\left({\bf r}\right)=\left(0,\sin q_{0}x,\cos q_{0}x\right)$. Hence, cholesteric helices are deformed, possibly leading to extended structures called helicoids, or to localized cholesteric bubble domains, called spherulites, which have been considered in Spheruliti ; PRE .

The special symmetry reduction (constant $\phi=-\pi/2$ and $y$ invariance)

[TABLE]

for the director field leads to the cholesteric finger phase with its axis along the $\mathbf{x}$ direction and it simplifies the Frank-Oseen energy (3) to

[TABLE]

where the constant magnetic contribution of the nematic phase has been subtracted.

The corresponding equilibrium equation is the elliptic sine-Gordon,

[TABLE]

where $\Lambda$ is the reciprocal of the magnetic coherence length. Please note that this equation does not depend on $q_{0}$ as the term $Kq_{0}\partial_{x}\theta$ in (5) is actually a null Lagrangian.

Large classes of solutions in the plane to equation (6) are well known in the literature Borisov3 (and references therein), but here we are dealing with different boundary conditions. Our reduction is the most natural extension of the problem considered in 17 and studied also in Lee92 in a linear setting. Here we study the nonlinear problem. Strong homeotropic anchoring conditions require $\Theta\left(x,\pm\frac{L}{2}\right)=2k\pi$ and $\Theta\left(x,\pm\frac{L}{2}\right)=2k^{\prime}\pi$, with $k,\,k^{\prime}\in\mathbb{Z}$, for negative and positive $x$ respectively. Thus, any non constant solution must have at least a jump singularity on the boundaries.

III $\bm{2\pi}$ - helicoids

In order to embody the above requirements, a known group-theoretical procedure Kaptsov suggests a simple ansatz for the solution, corresponding to the form

[TABLE]

which separately depends on the $x$ and $z$ variables (see also Burylov ).

Plugging this expression into (6) and recalling the trigonometric identity $\sin(\arctan{x})=\frac{x}{\sqrt{1+x^{2}}}$, we get

[TABLE]

where $(\cdot)^{{}^{\prime}}$ and $(\cdot)^{{}^{\prime\prime}}$ indicate the first and the second derivatives with respect to the argument, respectively.

The integration procedure of this equation is well known Burylov , but for self-consistency of the article and for convenience of the readers, we just sketch the main steps.

Then, performing the mixed derivative ${\partial}^{2}_{xz}$ to both sides of eq. (8), the r.h.s. vanishes and it becomes

[TABLE]

Last equation can be separated in

[TABLE]

[TABLE]

where $a$ is a separation constant.

Integrating eq. (10) and (11) with respect to $x$ and $z$, one obtains

[TABLE]

[TABLE]

which can be further integrated (i.e. multiplying the equation for $Z(z)$ in (13) by ($Z^{\prime}/Z^{3}$)) to get

[TABLE]

[TABLE]

Equations (14) and (15) contain five constants (the separation constant $a$ and the integration constants $b,c,d,g$), but only three of them are independent. Indeed, substituting them in eq. (8), one can easily obtain the relations

[TABLE]

Thus, equations (14) and (15) can be written in the final form

[TABLE]

[TABLE]

which are solvable in terms of Jacobi elliptic functions.

The step-like conditions on the boundaries force to have $Z\left(\pm\frac{L}{2}\right)=0$ and to look for functions $X\left(x\right)$ monotonic and unbounded, also at some finite point. Then, one obtains

[TABLE]

where $a>0$ and $\Lambda^{2}>d$.

As for $Z(z)$, the compatibility of equation (17) with the sine-Gordon (6) and the above boundary conditions, set $d=-\frac{\pi^{2}}{L^{2}}$, in order to have real $\theta$ with the semi-period of the corresponding $Z(z)$ exactly equal to the thickness $L$ of the sample. Thus, we are led to the general expression

[TABLE]

for $n\in\mathbb{N}$, and where the effective scale $\ell=\Lambda\,L$ and the modal factor $c_{n}=\left[1+\frac{\left(1+2n\right)^{2}\pi^{2}}{\ell^{2}}\right]^{{\scriptstyle{\frac{1}{2}}}}$ have been introduced. The asymptotic behaviour of this family of solutions is $\theta\stackrel{{\scriptstyle x\to\pm\infty}}{{\longrightarrow}}{}\mp\pi$, that is, the director field $\mathbf{n}(\mathbf{r})$ rotates by $2\pi$ over the strip, as depicted in Fig. 1.

We first observe that the solution only depends on the structural parameters $L$ and $\Lambda$ via their product $\ell$. Second, there exists an entire spectrum of excitations indexed by $n\geq 0$. Since for $n>0$ solutions are discontinuous at $x=0$ and $\cos(\pi(1+2n)z/L)\leq 0$, the most physically meaningful solution corresponds to $\theta_{0}\left(x,z\right)$ ( Fig. 1), which is continuous and differentiable at all points of the strip but at the special points $\left(x=0,z=\pm\frac{L}{2}\right)$. Here it exhibits a discontinuity of $2\pi$, indicating the presence of a disclination along the $\mathbf{y}$ direction.

Whilst the dependency on $z$ is determined uniquely by $L$, on the $x$ variable the typical scale is $\Lambda^{-1}\left({1+\frac{\pi^{2}}{\ell^{2}}}\right)^{-{\scriptstyle{\frac{1}{2}}}}$, which comes from the tendency of the external magnetic field to align the molecules along its direction, further enhanced by the anchoring of the bounding surfaces.

We also notice that solution (LABEL:sol2pihelic) can be obtained by applying the nonlinear Bianchi superposition formula Bianchi for two 1-kink solutions of the elliptic sine-Gordon (6). Such a theorem states that if $\Theta$, $\Theta_{1}$ and $\Theta_{2}$ are three solutions related by the overdetermined first order system, the so-called Bäcklund transformation,

[TABLE]

($\beta_{i}\in\mathbb{C}$, $i=1,2$), then a fourth solution is given by

[TABLE]

As an application, one solves the system (20), (21) corresponding to the trivial solution $\Theta=0$. This leads to the kink solutions $\Theta_{i}=4\arctan\left(e^{a_{i}+\frac{\beta_{i}\Lambda\zeta}{2}+\frac{\Lambda\bar{\zeta}}{2\beta_{i}}}\right)$, with $a_{i}$ integration constants for $i=1,2$. Combining these two expressions into formula (22) and requiring real solutions for real variables, one is led to the constraints $\beta_{2}=\beta_{1}^{-1}$ and $a_{1}=0,\;a_{2}=\imath\pi$. Furthermore, imposing the boundary condition one arrives at the expression (LABEL:sol2pihelic) with $\beta_{1}=c_{n}+\sqrt{c_{n}^{2}-1}$.

This result suggests that multiple $2\pi$-helicoid solutions may be built up by iterating the process. Unfortunately, the boundary conditions we impose are too rigid and they lead to singular solutions, which will be excluded from further discussions.

Finally, let us observe that since the total change of $\theta$ is $2\pi$, we can distinguish two regions where the director rapidly passes through $\frac{\pi}{2}$, sandwiched by three different regions in which the liquid crystal is close to the uniform state. Around the disclination, $\theta_{0}$ has a conformally invariant behaviour which, in cylindrical coordinates $x=\rho\cos t,\;z=-\frac{L}{2}+\rho\sin t$, can be expressed as a series with respect to $\rho$, namely,

[TABLE]

Thus, $\theta_{0}$ is almost independent of the distance $\rho$ from the disclination as shown by the contour plot in Fig. 2. The corresponding singularity in the energy signals a loss of order, and therefore, physically, in the vicinity of the disclination, the liquid crystal is melted from the ordered to the disordered isotropic phase. To avoid the singularity, we single out a small semi-disk, of radius $a\ll\frac{L}{2}$, around the disclination and replace it with a region in the isotropic phase, where it is common practice to uniformly impose an energy density cut-off ${\mathcal{E}}_{max}\approx\frac{K}{a^{2}}$ with $a\approx 10^{-2}-10^{-3}L$ Lee92 .

By using analytical solution (LABEL:sol2pihelic), we can write down an expression for the difference of energy with respect to the nematic phase in the limit of small $a/L$, which, up to second order, reads

[TABLE]

where $G$ is the Catalan constant. Please note that $E_{2\pi}$ corresponds to the energy of the $2\pi$-helicoids while $E_{\rm nem}$ stands for the energy of the nematic phase, which equals zero as can be seen from eq. (5).

At thermal equilibrium, stable 2$\pi$-helicoids may exist if $\Delta E_{2\pi}<0$; from the latter condition, one can infer a critical transition value for the chiral twist parameter as a function of the external field $\Lambda$. Notice that due to the presence of the external field, expansion (24) for small $a/L$ holds as far as $\Lambda L$ is not very high, which is true in the range of interest. By neglecting the $(a/L)^{2}$ contributions in (24) and introducing the cut-off energy, we arrive at

[TABLE]

with

[TABLE]

the critical chiral twist. Since the disclination radius has an important effect in locating the phase transition, $q_{H}$ might be used to estimate the effective size of the defects appearing in real liquid crystal samples.

IV $\bm{\pi}$ - helicoids

In a “more elementary” class of helicoids the director ${\bf n}\bm{(r)}$ twists only once, when going from the vacuum state $\theta=\pi$ at $x\to-\infty$, to the nearest $\mathbb{Z}_{2}$ equivalent one, namely $\theta=0$ for $x\to+\infty$. Moreover, because of the homeotropic anchoring, the boundary conditions on the confining surfaces will be

[TABLE]

Then, in such a kind of solutions, which we call $\pi$-helicoids, a couple of disclinations appear on the boundaries, parallel to the $y$ axis.

The analysis of the full nonlinear problem (6) with boundary values (27) is quite involved, even if one would profit from the integrability properties of the sine-Gordon equation through its Lax pair formulation, as developed in AKNS and references therein (see Appendix for more details). Actually, the Lax pair problem is a matrix linear representation of the overdetemined Bäcklund system (20) and (21), and it is equivalent to finding a matrix integration factor $\Psi\left(x,z;\lambda\right)$, depending on an auxiliary complex parameter $\lambda$, which makes the following 1-form exact

[TABLE]

with ${\hat{\sigma}_{3}}\cdot=\left[\sigma_{3},\cdot\right]$, $\sigma_{i}$ being the Pauli matrices, $\Omega(\lambda)$ and $\omega(\lambda)$ the dispersion functions on $\mathbb{C}_{\lambda}/\{0\}$ defined as

[TABLE]

and $Q(x,z,\lambda)$ the potential matrix

[TABLE]

Then, one can check that the closure condition of $W$ implies that $\Theta$ has to satisfy the sine-Gordon equation (6). However, this type of construction also requires the knowledge of the derivatives normal to the boundaries, possibly containing singularities not included in (27). These difficulties are not alleviated in the limit of the linear approximation, which can be obtained by letting $\Theta=2\theta\;\approx\;0$ and $\Psi\;\approx\;\mathbf{1}_{2}$, namely

[TABLE]

To avoid such an addition of extra information about the derivatives, a long procedure was introduced in FLP . Such a work led to a unifying approach in solving both linear and nonlinear boundary value problems, via a suitable integral mapping of the boundary data into the solution Unified .

The closure condition of $W_{\rm lin}$ yields the modified Helmholtz equation which can be solved by different types of transform methods, thus bypassing the above mentioned problem of the derivatives on the boundary. By exploiting the discrete mirror symmetry along the x-axes of our problem (6) and (27), the modified Helmholtz boundary value problem on the semi-strip arises

[TABLE]

In principle, this system could be solved by standard Fourier series methods Morse . Because of the linearity of the equation and the separability of the independent variables, one can look for solutions of (32) in the form $\theta_{+}\left(x,z\right)=\sum_{k=0}^{\infty}X_{k}\left(x\right)Z_{k}\left(z\right)$, where $Z_{k}$ belong to an orthonormal basis of the $L^{2}_{\left[-\frac{L}{2},+\frac{L}{2}\right]}$ space. Taking into account the vanishing of the series at $|z|=\frac{L}{2}$, one sets $Z_{k}=\sqrt{\frac{2}{L}}\cos\left(\frac{\pi\left(2k+1\right)z}{L}\right)$. Assuming that the series is absolutely and uniformly convergent together with its derivatives, its substitution into (32) can be computed by exchanging the sum with the laplacian operator, which will act term by term. Using the ortho-normalization of the $Z_{k}$, for the generic coefficient one gets the equation

[TABLE]

The corresponding solutions are real exponentials, from which we only keep those vanishing as $x$ tends to infinity, namely $X_{k}=\xi_{k}\exp\left[-\frac{x\sqrt{\pi^{2}(2k+1)^{2}+\Lambda^{2}\,L^{2}}}{L}\right]$. The amplitudes $\xi_{k}$ are simply determined by letting $x{\rightarrow}0$, and using again the ortho-normalization of the $Z_{k}$ for the series $\sum_{k=0}^{\infty}\xi_{k}Z_{k}=\frac{\pi}{2}$, we are led to the solution

[TABLE]

This latter cannot be cast in a factorized form as in (LABEL:sol2pihelic), unless when $\Lambda=0$. In this case, the series can be summed up to the closed form

[TABLE]

which, upon using the trigonometric identity $\arctan u+\arctan v=\arctan\left(\frac{u+v}{1-uv}\right)$, matches Eq. (10) in Ref. 17 up to a global sign due to the equivalence $\theta\rightarrow\pi-\theta$ in order to comply with their far-field values.

Using expression (32b) one can continuously complete the solution also for $x<0$, except at the points $\left(x=0,z=\pm\frac{L}{2}\right)$. The exponential decaying in $x$ has the characteristic length $\Delta x\approx\frac{2L}{\sqrt{\ell^{2}+\pi^{2}}}$.

However, close to $x=0$ the linear approximation to the sine-Gordon might be unsatisfactory, first because $\sin\frac{\pi}{2}\neq\frac{\pi}{2}$, and second because there the series (34) is not uniformly and absolutely convergent, implying all well known convergence problems at the discontinuity points. Moreover, the partial derivatives ${\partial}_{x}\theta\left(x,z\right)$ and ${\partial}_{z}\theta\left(x,z\right)$ are diverging for $x\to 0^{\pm}$ and $z\to\pm\frac{L}{2}$. So the differentiability of $\theta$ on the segment $\{x=0,|z|<L/2\}$ is lost.

By finding an integration factor associated to the 1-form $W_{\rm lin}$ in eq. (31), we can also adopt a different representation of the same solution (see Appendix for details) as follows AF

[TABLE]

where

[TABLE]

[TABLE]

encode the information about the boundary conditions. In particular, the function $G_{1}\left(\lambda\right)$ has the poles ( $\ell=\Lambda L$)

[TABLE]

As expected, the corresponding residues lead to expansion (34), this bringing both approaches together. On the other hand, by looking at the analytic properties of the integrands in the lower half-plane $\mathbb{C}_{\lambda}$ , one may suitably deform the integration contour in formula (36), in order to extract more information about the solution near the

singularities. This alternative approach might be the starting point to study the full nonlinear problem, which still presents difficulties to get an explicit formula FLP (in particular, the singular behaviour of the solution around the disclination).

In Fig. 3, $\mathbf{n}\left(x,z\right)$ is represented by using formula (34) up to $k=10$. At $\left(0,\pm\frac{L}{2}\right)$ an overlap of directional ellipsoids occurs, thus giving rise to disclinations extended in the orthogonal direction $y$. However, the actual configuration is sensitive to the involved parameters, which control also how good the linear approximation is. To have more insight in the $\pi$-helicoids, a numerical boundary value problem solver has been implemented for the full nonlinear system. We employed standard numerical techniques, namely the central finite difference scheme accompanied by appropriate application of Newton’s iterative method for the calculation of the function $\theta(x,z)$ over a suitable grid Press2007 ; LeVeque2007 . The problem can be coded in almost any programming language. However, we used MATLAB© by Mathworks matlab because of the natural way it operates with large and sparse matrices.

In Fig. 4 a comparison between the numerical solutions and the analytical series up to 7000 terms is presented. Although the similarity of the profiles is remarkable at this accuracy, a study of the Laplacians shows the expected differences associated to the different equations solved (either modified Helmholtz or sine-Gordon).

V Phase transition diagram

In order to assess the relative stability of the solutions we have found above, i.e., $\pi$- and $2\pi$-helicoids, we need now to perform an energy analysis and to compare the solutions also with the uniform nematic configurations to see which one is energetically favoured in terms of the physical parameters involved, namely $q_{0}$ and $\Lambda$. To this purpose we use solutions (LABEL:sol2pihelic) for $n=0$ and (34) with its extension for negative $x$, where in the latter we truncated the series at $k_{\max}=7000$.

Recalling that the chiral strength $q_{0}$ is a function of the temperature (linear under certain conditions), one could interpret the phase diagram as the result of thermal-magnetic competitive effects in the formation of the helicoids. When the energy of both solutions is greater than zero, the homeotropic nematic phase is the favoured one. It is important to note that the transitions between the three different configurations occur along two different curves (see Fig. 5). There are two different thresholds in chirality strength, for fixed magnetic field, also depending on the disclination size $a/L$. In Fig. 5, we have set $a/L=10^{-2}$, although for smaller values of $a/L$ the shape of the diagram does not change. In particular, if a study of a pair of $\pi$-helicoids as a function of the distance between them is realised, one can see the energy increasing when they approach. This implies the existence of a barrier preventing the two helicoids from sitting on top of each other, unless some energy is introduced into the system from outside. This underlines the stability of the $2\pi$-helicoids with disclinations located at the boundaries with respect to the decay into two $\pi$-helicoids with disclinations placed at the same point, which correspond to completely different configurations. Conversely, this is not the case when $\Lambda=0$ where the two configurations are equivalent.

Moreover, the curve between the nematic and $\pi$-helicoid phases should be the analogous of the straight line $\Lambda_{0}=\frac{\pi|q_{0}|}{2}$ in the bulk model ABC . Deviations from such a behaviour are related to the anchoring, possibly leading to significant variations in the value of the critical external field. Moreover, in Lee92 it was found a semi-empirical coexistence curve between the homeotropic and the cholesteric phases, providing the critical external field $\Lambda_{c}L$ in terms of the cholesteric twist $q_{0}L$. Using the present notation, such a relation reads

[TABLE]

with $\gamma=K_{2}/K$ and

[TABLE]

where $W_{k}$ is the $k$-th branch of the Lambert $W$ function DLMF . A fitting procedure to our numerical data leads to $\gamma=1.054$, in excellent agreement with our assumed one constant approximation (i.e. $\gamma=1$). The corresponding best fit curve is displayed in Fig. 5 (black dashed line).

VI Conclusions

Summarizing, we have analytically found $2\pi$-helicoids which, to our knowledge, are novel configurations in bounded CLCs with homeotropic anchoring, allowed by the nonlinearity arising from the presence of an external field. Analogously, we studied a different type of helicoids, the $\pi$ ones. If the former configuration can be derived in a closed form, for instance by the Bäcklund transformations, for the latter we have obtained approximated expressions and numerical solutions. Both classes of configurations are characterised by the presence of disclinations, located at the boundaries. The disclinations imply energy density divergences, which may be overcome by introducing a phenomenological energy cut-off, corresponding to excluded regions of melted crystal.

Numerically, we provided a phase diagram in the parameters $q_{0}L$ and $\Lambda L$, which established the energetically favoured configurations among them and the uniform nematic phase and the corresponding transitions. We showed that $\pi$-helicoids switch to $2\pi$-helicoids under certain circumstances dictated by the parameters of the problem. In particular, for fixed value of the external field, we found the sequence uniform-to-$\pi$-to-2$\pi$-helicoids as the value of the chirality increases, as expected.

We would like to specify that other configurations, like the 3-dimensional objects discussed in the Introduction and helicoid lattices, might modify the phase diagram presented in Figure 5. However, it is not necessary to introduce them in order to characterize the cholesteric-nematic transition we found, as shown by the agreement of our results with those obtained in Lee92 , and the nature of the disclinations at the boundaries. This approach is similar to previous studies of analytical solutions limited to single helicoids with disclinations but in the absence of external fields and under the assumption of y-invariance (see for instance goosens ; scheffer and also 17 ; press1976 ). Preceding detailed energy computations are now extended to the case of external fields, in order to draw a theoretical phase diagram which can be considered equivalent to the experimental or semi-empirical ones presented in Lee92 ; oswald1 ; oswald90 .

This work provides a step forward in the analytical construction of the fully comprehensive phase transition diagram for geometrically frustrated chiral nematics, which, to the best of our knowledge, is still a crucial and formidable problem for liquid crystals theory. In this spirit, the next step to be taken would be the study of the existence of lattices of helicoids, their interactions and the corresponding phase diagram.

VII Acknowledgments

The authors are thankful to the Referees for their comments and suggestions. CN is supported by the INFN grant 19292/2017 Integrable Models and Their Applications to Classical and Quantum Problems. LM has been partly supported by INFN through the MMNLP project by CSN4.

Appendix A

In this Appendix we give details about the derivation of the integral solution (36) to the modified Helmoltz boundary value problem and its relation to the series representation (34).

The sine-Gordon Eq. (6) is well known to be solvable via the Inverse Spectral Transform (IST) AKNS , for initial data given at a fixed value of one of the independent variables. The IST is an integration procedure, which consists in studying the analytical properties of the wave-functions of a specific linear differential operator (sometimes called the principal Lax operator) in one independent variable. Then, analytical deformations in the other independent variable of those wave-functions are obtained by the action of a second suitable differential linear operator. In our case, these operators, called Lax pair, possess the key property to have the sine-Gordon as their compatibility (commutativity) condition. This formulation introduces a complex spectral parameter $\lambda$, analogous to the momentum/frequency in the Fourier Transform method. In fact, one may encode information of the solution (like initial data and/or boundary values) into certain specific functions of $\lambda$: the so-called scattering data. The second Lax operator leads to linear differential equations determining the deformations of the scattering data. This new data is used to reconstruct the solution of the sine-Gordon for any value of the pair of independent variables, by solving the so-called inverse linear spectral problem. Moreover, the scattering data is determined by the singularities of the principal Lax operator wave-functions in the complex $\lambda$-plane. Under suitable asymptotic assumptions on the sine-Gordon solutions, independent wave-functions are sectionally holomorphic in the $\lambda$-plane and the corresponding jump functions, along the curves separating the holomorphicity regions, provide the scattering data. Now, in classical complex function theory, solving the Riemann-Hilbert problem means reconstructing these sectionally holomorphic wave-functions from given jump functions AblowitzClarkson . This is an equivalent way to formulate the inverse spectral problem. This idea can be readily applied also in the elliptic case, with rather general boundary conditions, provided that suitable jump functions are assured Unified ; AshF . For the sine-Gordon on the semi-strip this approach was almost fully developed in FLP (and references therein), where a unified method both for linear and nonlinear integrable equations has been developed. Thus, this includes the treatment of the linear (modified) Helmholtz equation, which is exactly the linear approximation to our original model (6). Thus, in the present appendix we will apply such methods to extract interesting information in the linear case given by Eq. (32a). This simplification is supported by numerical calculations (see discussion at the end of section IV). In fact, it turns out that linear and nonlinear treatment show quite small differences in the quantitative behaviour of the solutions near the disclination and for values of $\theta\approx\frac{\pi}{2}$, where we are far from linear scheme.

As mentioned in the main text, in the unified approach to the study of the sine-Gordon and the modified Helmholtz equation, the traditional Lax pair formulation is equivalent to finding a matrix integration factor $\Psi\left(x,z;\lambda\right)$ (the wave-function), depending on $\lambda$, which makes the 1-form $W$ in Eq. (28) exact.

The first observation is that $W_{lin}$ in Eq. (31) remains closed by adding a suitable exact 1-form $d\left(e^{\Omega\,x+\omega\,z}\kappa\left(x,z\right)\right)$. In particular, one can choose the function $\kappa\left(x,z\right)$ in such a way to cancel the ${\partial}_{z}\,\theta$ (${\partial}_{x}\,\theta$) terms from the $dx$ ($dz$) component. A special case is given by

[TABLE]

where $\Omega$ and $\omega$ are as in Eq. (29).

Suppose now that the boundary corresponds to a semi-strip where the value for $\theta$ is assigned on its three sides. Then, in analogy with the equation (28), we look for three functions $\Psi_{j},\;j=1,2,3$ such that

[TABLE]

By integrating along the paths shown in Fig 6 (as many as the sides of the semi-strip) with the initial conditions

[TABLE]

for any $\lambda$ one finds the functions $\Psi_{j}$, which are related by

[TABLE]

The consistency of these equations requires the global condition

[TABLE]

If the analyticity properties of the functions $\Psi_{j}$ are determined and the r.h.s. of (45)-(47) are known, such a system will define a Riemann-Hilbert problem associated to the original modified Helmholtz BVP.

Proceeding in this way, one first integrates equation (43) along the path 1, obtaining

[TABLE]

Hence the constants appearing in (45) and (47) are given by

[TABLE]

Since $x^{\prime}\geq 0$ the above integrals are convergent for $\Im\lambda>0$.

Now, in (46) the quantity $\Psi_{3}\left(0,-\frac{L}{2};\lambda\right)$ can be computed again by integrating (43), for $j=3$, along the segment $\left[\left(0,\frac{L}{2}\right),\left(0,-\frac{L}{2}\right)\right]$. It yields

[TABLE]

Integrating by parts the last double integral

[TABLE]

the previous expression reads

[TABLE]

for any $\lambda\in\mathbb{C}$.

Although the constants introduced in (45)-(47) contain information about the boundary value problem, they still require unknown data. Hence, further manipulations are needed in order to suppress them.

Turning again our attention to the original problem (32), Eq. (50) and Eq. (A) become

[TABLE]

[TABLE]

both still containing unknown functions.

To have more restrictions on them, we look at their symmetries. First, the $z$-mirror symmetry

[TABLE]

implies

[TABLE]

for (54), which, used in the global symmetry (48), leads to

[TABLE]

Now, since $\Psi_{1}\left(0,-\frac{L}{2};\lambda\right)$ depends on $\lambda$ only through $\Omega\left(\lambda\right)$, it will enjoy the same inversion symmetry $\lambda\to-\frac{1}{\lambda}$, namely

[TABLE]

On the other hand, $\Psi_{3}\left(0,-\frac{L}{2};\lambda\right)$ in (A) depends on $\lambda$ through both $\omega$ and $\Omega^{2}$, thus being invariant under $\lambda\to\frac{1}{\lambda}$

[TABLE]

Applying the transformation $\lambda\to\frac{1}{\lambda}$ into (58) and (59), substituting into (60) and rearranging the various terms we are led to the equation

[TABLE]

which, because of the convergence region for $\Psi_{1}$, only makes sense for $\lambda\in\mathbb{R}$ with

[TABLE]

Analogously, using again (58) and the above symmetries, one obtains

[TABLE]

which holds on $\lambda\in\mathbb{C}^{-}$ and where

[TABLE]

Both $G_{1}$ and $G_{2}$ now only encode information about the boundary conditions, but in order to use them one has to suitably modify the relations (45)-(47). Precisely, defining the new $\tilde{\Psi}_{j}$ (coherently with the definition (43))

[TABLE]

the relations (45)-(47) read

[TABLE]

The system (66) is a Riemann-Hilbert problem defined on three branches, where the jumps are completely known functions. Furthermore, going back to (49) and estimating its asymptotic behaviour for $\lambda\to\infty$, one obtains a reference value for $\tilde{\Psi}_{1}$

[TABLE]

Thus, from the jump conditions (66) and the asymptotic value (67), $\tilde{\Psi}_{1}$ is given by

[TABLE]

which holds for $\lambda\in\mathbb{C}^{+}$. On the other hand, from the component $dx$ of the linear problem (43), or from the solution (49), in the limit $\lambda\to 0$ one obtains

[TABLE]

Comparing the last two relations, one obtains the solution of the modified Helmholtz BVP (32) by the final formula

[TABLE]

One can check that the solution (70) certainly satisfies the boundary conditions in (32).

Now, if one sets $z=-\frac{L}{2}+\epsilon$ with $0<\epsilon<\frac{L}{2}$, it can be verified that the integrands in (36) are bounded and analytic functions in $C^{IV}$ . Thus, performing as above, the change of variable $\lambda\to-\frac{1}{\lambda}$ on the first term, we are led to an integrand split into an analytic part plus a meromorphic contribution $-\imath\frac{e^{\Omega\left(\lambda\right)x}}{2\lambda}\frac{\left(1-\lambda^{2}\right)}{\left(1+\lambda^{2}\right)}\frac{\left(1-e^{-\omega\left(\lambda\right)\epsilon}\right)}{\left(e^{L\omega\left(\lambda\right)}+1\right)}$, containing all poles in the family $P_{G}^{-}=\left\{\lambda_{\in}P_{G},\;\Im\lambda_{n}<0\right\}$. Thus, the solution is given in terms of the series of its residues. Noting that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

one can first sum up the contributions coming from the poles $\lambda_{n}$ and $\lambda_{-\left(n+1\right)}$, collecting the exponential $x$-dependence and the trigonometric $z$-dependence. This manipulation leads directly to the formula (34).