Lily-like twist distribution in toroidal nematics
Andrea Pedrini, Marco Piastra, and Epifanio G. Virga
Toroidal nematics are droplets of nematic liquid crystals in the form of a circular torus. When the nematic director is subject to planar degenerate boundary conditions, the bend-only director field with vector lines along the parallels of all internal torodial shells is an equilibrium solution for all values of the elastic constants. Local stability analyses have shown that an instability is expected to occur for sufficiently small values of the twist elastic constant. It is natural to conjecture that in this regime the global equilibrium would be characterized by a maximum twist deflection on the boundary of the torus, with a twist distribution over the torus’ cross-section represented by a fennel-like surface. We prove that surprisingly the stable twist distribution is instead represented by a lily-like surface. Thus the overall maximum twist deflection falls well within the torus. To cope with the complexity of the elastic free-energy functional in the fully non-linear setting, we developed an ad hoc deep-learning optimization method, which here is also duly validated and documented for it promises to be applicable to other similar problems, equally intractable analytically.

footnotetext:  Dipartimento di Matematica, Università di Pavia, via Ferrata 5, 27100 Pavia, Italy. E-mail: andrea.pedrini@unipv.itfootnotetext:  Dipartimento di Ingegneria Industriale e dell’Informazione, Università di Pavia, via Ferrata 5, 27100 Pavia, Italy. E-mail: marco.piastra@unipv.itfootnotetext:  Dipartimento di Matematica, Università di Pavia, via Ferrata 5, 27100 Pavia, Italy. E-mail: eg.virga@unipv.itfootnotetext: ‡ All authors have contributed equally to this work.

1 Introduction

Toroidal nematics are typically droplets of nematic liquid crystals in the shape of a circular torus surrounded by an isotropic fluid that enforces degenerate planar anchoring conditions for the nematic director . We imagine that the droplet’s shape is somehow prescribed and inquire about the most stable equilibrium nematic structure inside it. This problem was studied by Koning et al.1 who envisaged a first destabilization mechanism by twist-injection from the surface into a bend-only nematic texture.

The elastic free energy stored in is given by the functional


where is the volume element and , , , and are Frank’s elastic constants, weighting the costs of splay, twist, bend, and saddle-splay distortions, respectively (as explained in a number of textbooks, 2, 3, 4 which the reader is referred to for further details). Frank’s energy functional in (1) measures the elastic cost involved in distorting a uniform state where . As first remarked by Ericksen, 5 this functional, which vanishes on any uniform state, is positive semi-definite only if the elastic constants obey the inequalities


which are also referred to as Ericksen’s inequalities. Here, for convenience, we shall scale all elastic constants to , assumed to be positive,


so that Ericksen’s inequalities become


As is well known, the energy in (1) is a null-Lagrangian and can be transformed into a surface energy, so that takes the equivalent form


where is the area element and is the outer unit normal to .

We want to find the director field that minimizes subject to the boundary condition


As stated, this is a formidable variational problem. Luckily, a number of partial results are known from the literature. First, there is a special, bend-only director field whose vector lines are parallel circles with centers on the straight axis of the torus, also called the parallels of . This field is a universal solution6, 7, 8 for the Euler-Lagrange equation associated with , for any choice of the elastic constants. This means that , which complies with (6), is also an equilibrium configuration for our variational problem. It remains to be seen whether is the absolute minimizer of as well.

The local stability of for toroidal nematics was first studied by Koning et al.1 and then further refined.9 In brief, for any given , becomes unstable for sufficiently large, that is, when the cost of bending is sufficiently greater than the cost of twisting. The critical values of increases with decreasing , thus suggesting that this is a surface-dominated transition, as (5) identifies as a surface-like elastic constant. When becomes unstable a certain degree of twist is likely to be injected into the torus, but nothing is known about how the newly acquired twist is distributed inside the torus: the local stability analysis cannot answer this question. One other thing we know from symmetry: the chiral degeneracy of the system entails that the energy minimizers will come in pairs, with equal energy but opposite chirality.1, 9

In this paper, we solve the fully non-linear problem and find the optimal twist distribution in toroidal nematics when is unstable. If the interpretation of as being surface-driven were to be supported by the fully non-linear results, one would expect the minimizers of to exhibit the maximum twist deflection at the boundary of the torus. On the contrary, we shall see that this is not the case, as falls generally in the interior of : more precisely, on every circular cross-section of , the twist angle distribution has typically a lily-like appearance, as opposed to the fennel-like surface that would describe the twist distribution had fallen on .

In other words, in toroidal nematics, the twist favoured by surface elasticity, does not extend freely up to the boundary, as it would happen in a cell bounded by planar plates, but it remains self-confined in the interior, giving rise to a characteristic twist structure which has the potential of being detectable optically. No doubt all this is conjured by the toroidal geometry, whether it is typical of that we cannot tell.

The paper is organized in the following way. Section 2 illustrates the special class of nematic director fields within which we minimize ; in this class, which includes as well as quite different fields, we write a reduced energy functional. In Section 3, we summarize the deep-learning method applied to this minimization problem, highlighting the mathematical infrastructure behind the code that we devised for this problem. In Section 4, we illustrate the results obtained through our method and comment on their physical significance, pausing on the details of the generic twist structure characterizing the states with minimum energy. Section 5 is devoted to a number of validation tests for our computations, both numerical and analytical in character. In particular, we determine numerically the curve in parameter space marking the instability onset for the bend-only solution and contrast it against the known analytical estimates.1, 9 Finally, in Section 6, we draw the conclusions of this work and we discuss possible extensions to other soft matter systems of the deep-learning method adopted here. The paper is closed by an Appendix, where we outline the features of the accompanying code.

2 Reduced energy functional

A toroidal nematic is described as a circular torus . We denote by and the radius of the inner circular axis of and the radius of its meridian cross-sections (see Fig. 1). We employ toroidal coordinates to describe points in : are polar coordinates in the meridian cross-section of with a plane making the angle with a given reference plane, taken as the plane of a Cartesian frame (see Fig. 1). Thus ranges in the interval and both and in . The coordinates are orthogonal and the associated frame is positively oriented, that is, . In these coordinates, the volume and area elements are given by

Fig.  1: (Colour online) Representation of a circular torus (left) obtained by revolving a disk of radius (bottom-right) around the -axis. The centreline of has radius . Each point of the torus is described by the toroidal coordinates and the nematic vector field at is represented in the associated toroidal frame (top-right) through the angle , as in Equation (8).

The special class of director fields that we shall consider in are represented by


where . These fields are characterized by being independent of the azimuth and having zero radial component. Moreover, we shall assume that to avoid a disclination of along the centreline of (which would cost infinite energy in Frank’s model). It follows from (8) that for and that reversing the sign of reverses the chirality of , that is, its sense of winding around the torus.1

To make our variables dimensionless, we shall scale lengths to , so that


and we shall call the ratio


Following Koning et al. 1, we further restrict the admissible class of director fields by requiring that , which amounts to the following differential equation for ,


which can be easily integrated, delivering


where is a real function of only. For (12) to obey the inequality , must be such that for all , which is satisfied whenever


Thus, eventually, the energy functional can be reduced to a functional in the single scalar function subject to (13) and


Making repeated use of Appendices A2 and A3 of Pedrini and Virga9, and exercising patience, we finally arrive from (5) to the following formula for ,


where .

Clearly, , showing that minimizers of come in pairs of functions differing in sign. We shall systematically consider only one member of such a pair, keeping in mind that it also represents its companion with opposite sign (and chirality).

The (dimensionless) energy of is readily obtained from (2) as


The excess energy associated with a generic will then be given by


Formula (2) has the advantage of being valid for all director fields in the class (8), no matter how they differ from , but it is intractable analytically. We devised a numerical scheme to minimize , which is outlined in the following section.

3 Deep-learning method

In general, deep-learning methods are constructed to find smooth functions (typically, of class ) that minimize a given cost functional .333All definitions and results presented in this section are more generally valid for all spaces . Formally, we shall write that


In general, denotes the set of functions in where attains its minimum. Here, however, we shall adopt the simplifying assumption (corroborated by physical intuition) that, possibly apart from a geometric degeneracy (such as the sign parity of both and ), has a unique minimizer.

Deep-learning methods are based on a nested (or deep) representation of such as


where , are vectors and are matrices, collectively represented by the set of parameters . Here is the overall number of scalars contained in and is the depth of the representation. The function is applied element-wise to its vector and matrix arguments.10

Representation (19), which is also called deep neural network (DNN), is often said to be universal. Informally, this means that (19) is a versatile and adaptable representation of a function through a finite number of scalar parameters. More technically, it has been proved11 that for all , all compact subset of , and all , there exist and such that


where and denotes the derivative of .

Adopting the above representation, we change the problem of minimizing over into that of finding


where, to ease the numerics, the cost functional is evaluated on a discretised function defined over a finite dataset of input points. Thus, we reduce a minimization problem in a function space into an optimization problem in a finite, -dimensional parameter space. Clearly, by construction, this method is only capable of identifying in an approximate fashion the smooth minimizers of the cost functional ; these latter are assumed both to exist and to bear physical meaning.

The deep-learning method outlined above is here applied to the problem at hand by taking in (17) as . Accordingly, is replaced by the discretized version , where is subject to the constraint


and evaluated over the dataset , with and . The optimization problem then becomes: find




so that the boundary condition (22) is automatically satisfied.

We leave out a number of interesting mathematical questions, such as the degree of regularity of the minimizers of and the convergence of both and as both and are increased. These issues, though neither trivial nor irrelevant, exceed the scope of our present work, which is more concerned with the physical relevance of the results and treats the mathematical tools employed here as learning experiments for an automatic resolution of our problem. This is why our code is made available in a workable fashion alongside this paper, see the Appendix for further explanatory notes. Here we are content to give some additional details on how the code was implemented to arrive at the results described in the following section.

3.1 Numerical optimization:

TensorFlow***TensorFlow, the TensorFlow logo and any related marks are trademarks of Google Inc. is a recent software tool12, 13 for the numerical solution to minimization problems such as (23). In particular, TensorFlow allows translating functionals like into a flow graph in which each node represents an elementary operation and each edge represents a dependency of an operator on its operands. In addition, TensorFlow allows computing the partial derivative of a flow graph with respect to any of its nodes via a technique called automatic differentiation. 14 The result of any such differentiation is yet another flow graph that shares the relevant parts of the flow graph being derived. In passing, we note that the derivative in Equation (2) was also computed in this way.

Flow graphs can also be used to compute the numerical values of functions and functionals by specializing variables and parameters to definite tensorial values. With TensorFlow, such computations can be carried on transparently on either standard CPUs or GPU-accelerated hardware.

The optimization process that, in the case at hand, aims to obtain minimizing values of the parameters is performed by applying specific optimization operators, also provided by TensorFlow, which simplify the process of taking derivatives with respect to and applying iterative, quasi-second order methods.15 Processes of the above kind can only lead to local minimizers for hence are dependent on the assigned initial values.

TensorFlow is an open source software package that runs on several different hardware platforms and operating systems.

3.2 Experiments

We adopted an approximator with and , . The discretized functional was obtained by using a dataset with and . Among several possible choices, we selected Adam16 as the quasi-second order method needed for the optimization. At the beginning of each run was assigned small random values. For deterministic repeatability, specific seeding policies were adopted for pseudo-random number generators. The most effective values for the hyperparameters , , and were determined experimentally.

All our experiments were carried out on a workstation Xeon(R) CPU E3-1240 v3 @ 3.40GHz x 8 with 16 GB of RAM, equipped with Quadro K2200 GPU with 4 GB of DRAM.

In the rest of the paper, will always denote the (presumably unique) minimizer found as a result of the numerical optimization process.

4 Results and discussion

We performed many different optimizations of for several choices of the ratio and the normalized elastic constants and ; we illustrate here only a few representative experiments to highlight the main properties of the minimizers .

Fig.  2: (Colour online) The minimizer of for , , and different choices of .
Fig.  3: (Colour online) The minimizer of for , , and different choices of .
Fig.  4: (Colour online) The minimizer of for , and different choices of .

As shown in Figures 2, 3 and 4, when Ericksen’s inequalities hold, is a concave function, nearly linear for small values of , attaining in general its maximum at some . The value of depends on the choice of , , and : it decreases as or increase (see Figs. 2, 3), and as decreases (see Fig. 4). As a consequence, the corresponding function (obtained from Equation (12)), which represents the twist angle in the toroidal frame , can attain its maximum within the torus .

Fig.  5: (Colour online) Lily-shaped surface that represents the plot of corresponding to through Equation (12). The dark (blue) line describes the maximum twist deflection along every fixed radial direction in the cross-section of the torus. In this experiment: , , and .
Fig.  6: (Colour online) The line of radial maximum twist, that is, the projection on the cross-section of the torus of the dark (blue) line of Figure 5, for , and .

An example in which is located in the interior of is provided by choosing , , and , as in Figure 5, where a lily-shaped surface represents the twist angle , for and varying in the cross-section of the torus. The dark line marked on this surface describes how the maximum twist deflection along every fixed radial direction varies with the angle . The projection of this curve on the torus’ circular cross-section designates the line of radial maximum twist (rmt); it encircles the centre of the cross-section, running through the whole interior of the torus, from the outer radius (at ) to the inner radius (at ), being closer to the latter end than to the former, see Figure 6.

Fig.  7: (Colour online) Contour plots of on the plane of the scaled elastic constants, for (left) and (right) and for (top), (center) and (bottom). The white bottom-left region in each panel is where the bend-only director field is stable. Darker shades of blue correspond to higher values of , and so the line of radial maximum twist deflection is closer to the boundary of the torus.

By increasing or the line rmt opens up and touches the boundary on the inner side of the torus, generally remaining well inside on the other side. In polar coordinates , the line rmt is represented by a function , which changes with the geometric parameter and the scaled elastic constants and . The contour graphs in Figure 7 show how depends on , for and ; for , we have chosen its three most significant values. The farther is from unity, the farther is the line rmt from the boundary of the torus: visually, these cases correspond to the lighter shades of blue in Fig. 7, whereas the the darker shades correspond to lines rmt closer to the boundary. Accordingly, the white regions on the bottom-left corners of the plots in Fig.7 represent values of for which is zero and the bend-only director field is stable. For a fixed , the shades of blue darken uniformly in the plane as the angle increases from to , in perfect accord with Figure 6: the line rmt is off-centre, closer to the inner equator than to the outer equator. The intermediate value is considered here both for completeness and for being somewhat special, as there, by Equation (12), . The countour plots of Figure 7 also reveal that for high values of the overall maximum twist deflection tends to fall on the boundary , suggesting that the inner twist is surface-driven. However, this tendency is significantly reduced upon increasing : then the inner twist would be surface-driven only for very small values of .

Finally, by both Equations (8) and (12), the knowledge of also allows us to describe the vector field that minimizes . Figure 8 shows the vector lines of for , and (same parameters as in Figure 5) on the internal toroidal shell at : they approach the parallels on the outer side and bend slightly towards the meridians on the inner side.

Fig.  8: (Colour online) Vector lines of the director field for the same as in Figure 5 (, , and ), wrapped around the internal toroidal shell at .

5 Validation

As a validation of the optimization method presented in Section 3 we tested it on the very special case . For small values of , indeed, the torus can locally be approximated by a cylinder, and so we can compare with the analytical solution , existing for only, given by Davidson et al.17

Fig.  9: (Colour online) In the limit the torus becomes a cylinder and converges to the expected in Equation (25). Here is an example for and : for the error is .
Fig.  10: (Colour online) In the limit the torus becomes a cylinder and converges to the expected in Equation (25). Here is another example for and : for the error is .

Figures 9 and 10 show the results of two series of experiments for : the former with and , and the latter with and . For , the error (where denotes the standard sup norm in ) is negligible in both cases.

Fig.  11: (Colour online) Comparison between instability conditions for (top), (center) and (bottom). The green straight line delimits the instability region found by Koning et al.1, while the red line delimits the instability region obtained by Pedrini and Virga9. The light blue area represents the refined instability region found by performing a grid search in the plane with strides and , respectively, and by marking as instable any point for which , where is the result of the optimization method described in Section 3.

A further validation of our method is shown in Figure 11. Here, we compare the instability region for the pure-bend director field found numerically by our optimization method with the instability conditions given analytically by Koning et al.1 and by Pedrini and Virga9, with different choices of admissible test functions. For , , and , we performed a grid search in the plane with strides in and in . Each point of the grid is marked as unstable whenever the corresponding has . For all values of , the instability region obtained in this way is an improvement on both regions predicted analytically. The graphs in Fig. 11 also show that the instability analysis of Pedrini and Virga9 was quite on the mark: it effectively provides a reliable stability criterion.

Finally, our optimization method proves effective even if less guided. Further experiments showed that suppressing the boundary condition imposed by Equation (24) (thus leaving the representation in Equation (19) unconstrained) does not significatively change the final result: remains essentially the same function, still satisfying to a very good approximation.

6 Conclusions

We studied toroidal nematics subject to planar degenerate anchoring for the director by solving numerically the fully non-linear problem for the minimisers of Frank’s elastic free energy. When the pure-bend director field (with vector lines along the parallels of all internal toroidal shells) becomes unstable, two twisted director fields (with opposite handiness) take its place. It is quite natural to expect that such an instability is surface-driven, and so the maximum twist injected in the director texture should appear at the boundary of the torus. Remarkably, however, this turned out not to be the case. Typically, the twist distribution inside the torus is instead well described by a lily-like surface. The radial maximum twist deflection takes place along an inner, off-centre ring of the torus’ circular cross-section, away from which the twist angle (radially) decreases in both directions: to zero on the centre of the cross-section and to a positive value on the cross-section’s periphery. We think that such a characteristic twist structure could be observed experimentally and might become the optical signature of toroidal nematics.

The elastic free-energy functional that describes toroidal nematics in the fully non-linear case is rather involved, especially when all elastic constants are considered as free parameters. To cope with such a complexity, we developed an ad hoc deep-learning optimization method that proved quite versatile and also reliable, in light of two independent validation tests, in which the outcomes of our method were confronted with analytical predictions available in the literature. We plan to exploit further this method, which is duly documented and made available on the web, in the study of similar problems, among which is primarily the stability analysis of toroidal nematics subject to strong anchoring conditions for the director. Our objective will be to show under what condition on the smallness of the twist elastic constant the characteristic lily-like distribution of twist found here also arises there (if it ever does). We would expect that the answer to this question could be significant for chromonic liquid crystals,18 for which is by far the smallest among all four elastic constants.

Conflicts of interest

There are no conflicts to declare.


The work of A.P. has been supported by the University of Pavia under the FRG initiative, meant to foster research among young postdoctoral fellows.

Appendix A Deep-learning code

The code for the optimization method is available at the following repository: https://bitbucket.org/AndreaPedriniUniPV/lily-like_twist_distribution_in_toroidal_nematics/

The package contains the main file run.py and a lib folder with all the ausiliary code; the code is written in Python444Python Software Foundation. Python Language Reference, version 3.x. Available at http://www.python.org and all further requirements are listed in requirements.txt.

A README.md file, also available on the web, contains all information and instructions needed to run properly the code, to change the values of parameters , , and , and to modify some advanced configurations of the program. It also describes four examples, chosen from the experiments run to produce the figures presented in the main body of the paper.


  • Koning et al. 2014 V. Koning, B. C. van Zuiden, R. D. Kamien and V. Vitelli, Soft Matter, 2014, 10, 4192–4198.
  • de Gennes and Prost 1993 P. G. de Gennes and J. Prost, The Physics of Liquid Crystals, Clarendon Press, Oxford, 2nd edn., 1993.
  • Virga 1994 E. G. Virga, Variational Theories for Liquid Crystals, Chapman & Hall, London, 1994.
  • Stewart 2004 I. W. Stewart, The Static and Dynamic Continuum Theory of Liquid Crystals, Taylor & Francis, London, 2004.
  • Ericksen 1966 J. L. Ericksen, Phys. Fluids, 1966, 9, 1205–1207.
  • Ericksen 1967 J. L. Ericksen, Trans. Soc. Rheol., 1967, 11, 5–14.
  • Marris 1978 A. W. Marris, Arch. Rational Mech. Anal., 1978, 67, 251–303.
  • Marris 1979 A. W. Marris, Arch. Rational Mech. Anal., 1979, 69, 323–333.
  • Pedrini and Virga 2018 A. Pedrini and E. G. Virga, Liq. Cryst., 2018, DOI:10.1080/02678292.2018.1495771.
  • Goodfellow et al. 2016 I. Goodfellow, Y. Bengio and A. Courville, Deep Learning, MIT Press, 2016.
  • Hornik 1991 K. Hornik, Neural Networks, 1991, 4, 251–257.
  • Abadi et al. 2015 M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irvin, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyal, P. Warden, M. Wattenberg, M. Wicke, Y. Yu and X. Zheng, TensorFlow: large-scale machine learning on heterogeneous systems, 2015, https://www.tensorflow.org/, Software available from tensorflow.org.
  • Abadi et al. 2016 M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, M. Kudlur, J. Levenberg, R. Monga, S. Moore, D. G. Murray, B. Steiner, P. A. Tucker, V. Vasudevan, P. Warden, M. Wicke, Y. Yu and X. Zhang, TensorFlow: a system for large-scale machine learning., 2016, http://arxiv.org/abs/1605.08695.
  • Bartholomew-Biggs et al. 2000 M. Bartholomew-Biggs, S. Brown, B. Christianson and L. Dixon, J. Comput. Appl. Math., 2000, 124, 171 – 190.
  • Ruder 2016 S. Ruder, An overview of gradient descent optimization algorithms, 2016, http://arxiv.org/abs/1609.04747.
  • Kingma and Ba 2014 D. P. Kingma and J. Ba, Adam: a method for stochastic optimization, 2014, http://arxiv.org/abs/1412.6980.
  • Davidson et al. 2015 Z. S. Davidson, L. Kang, J. Jeong, T. Still, P. J. Collings, T. C. Lubensky and A. G. Yodh, Phys. Rev. E, 2015, 91, 050501.
  • Masters 2016 A. Masters, Liq. Cryst. Today, 2016, 25, 30–37.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description