Neural networks determination of material elastic constants and structures in nematic complex fluids

The developed neural networks-based method for determining elastic constants is based on the combined modeling of (i) liquid crystal effective dynamics, (ii) light transmission, and (iii) supervised machine learning, as then applicable both to experimental or modeling data. The method starts by calculating a large number of time-dependent light beam transmittance functions I(t) during the relaxation of the NLC sample that correspond to different elastic constants, which we then use to train the neural network that recognizes elastic constants from a simulated measurement. Later, the well-trained neural network is used to predict elastic constants also from signals measured from real samples in the laboratory. A general overview of the method, which could, in principle, be used for any other experimentally relevant setup and determining other material parameters, is shown in Fig. 1. Simulation of LC profiles is explained in Methods.

Neural networks based method for determining elastic constants

The orientational dynamics of a nematic liquid crystal depends on the initial state, rotational viscosity \(\gamma _1\), cell thickness D and also on Frank elastic constants. Therefore, if the director field is initially deformed by short pulses of the electric or magnetic field, it will reconfigure back to the equilibrium state after the fields are turned off. Consequently, the intensity of the transmitted light through the liquid crystal cell between the crossed polarizers also varies during the relaxation and its time-dependence therefore indirectly carries information about the elastic constants. So the idea of this method is to use a neural network to identify elastic constants from the time-dependent intensities of transmitted light, regardless of the initial state of the director. However, to perform this task, the neural network must be trained on data, i.e. on many examples of time-dependent intensities I(t) and associated elastic constants. In the setup for the determination of the nematic elastic constants, we assume a specific cell geometry where the liquid crystal is confined in a thin cell of thickness \(D\sim 10\ \upmu \text {m}\) with strong anchoring which is uniform in x and y directions on both boundaries, so that the director, varies only along the z axis and in time t, \(\textbf{n}= \textbf{n}(z,t)=(n_x(z,t),n_y(z,t),n_z(z,t))=(\cos \phi (z,t)\cos \theta (z,t),\sin \phi (z,t)\cos \theta (z,t),\sin \theta (z,t))\), where \(\theta \in [-\pi /2,\pi /2]\) and \(\phi \in [-\pi ,\pi ]\) are spherical angles. For determining only splay and bend constants, \(K_{11}\) and \(K_{33}\), even more simple 2D director geometry (i.e. director profile variability), \(\textbf{n}(z,t)=(\cos \theta (z,t),0,\sin \theta (z,t))\), proves sufficient.

Determining \(K_{11}\) and \(K_{33}\)

The scheme of the method for determining \(K_{11}\) and \(K_{33}\) is shown in Fig. 2. To create a training set, we start with generating a pair of random elastic constants \(K_{11}\) and \(K_{33}\) from the uniform distribution in the interval of possible expected values. For example, to predict the elastic constants of 5CB liquid crystal at approx 23 °C, we choose elastic constants from the uniform distribution \(U(2\text { pN}, \ 18\text { pN})\) as we expect the actual constants somewhere within this interval^29,35,36. Next, the initial non-equilibrium state of the director \(\textbf{n}(z,t=0)=(\cos \theta (z,t=0),0,\sin \theta (z,t=0))\) is set by generating a random, not necessarily physically meaningful smooth 1D function \(\theta (z,t=0)\) for \(z\in [0,D]\) by quadratic interpolation between a random number of points at different random positions within the interval [0, D]. Having the initial state \(\textbf{n}(z,t=0)\), a pair of elastic constants and the rotational viscosity \(\gamma _1\) (e.g. \(0.098 \text { Pa s}\) for 5CB at the room temperature^37,38), the dynamics \(\textbf{n}(z,t)\) is numericallly simulated (see “Methods”). The boundary conditions \(\textbf{n}(z=0,t)\) and \(\textbf{n}(z=D,t)\) are determined by the selected anchoring type on each boundary. Different combinations of planar (\(\theta =0\)) and homeotropic (\(\theta =\pi /2\)) anchoring are tested. From the configuration of the director at each time step, the transmission of light through crossed polarizers and the sample of reconfiguring liquid crystal between them are calculated using the Jones matrix formalism^39,40. We comment that more advanced light propagation methods could also be used, such as finite difference time domain (FDTD)⁴¹, but give no qualitative difference for the considered nematic geometries. The values of the ordinary and extraordinary refractive indices, the thickness of the cell and the spectrum of the light source are required to be known precisely for particular liquid crystal material and experimental setting. In this way, the time dependence of the transmittance I(t) is calculated and discretized at, for example, 500 time steps for each pair of corresponding random elastic constants and the initial state of the director. Repeating this 200,000 times, a set of 200,000 pairs of input vectors \({\textbf{X}}_i=[I_i(t=0),I_i(t=\Delta t), …,I_i(t=T=499\Delta t)]\) and expected (true) output vectors \({\textbf{T}}_i=\left[ \widetilde{K_{11}^i},\widetilde{K_{33}^i}\right]\), packed as \(\{(\textbf{X}_1,\textbf{T}_1), …, (\textbf{X}_{200000},\textbf{T}_{200000}) \}\), is obtained, which is later split into a training and a validation data set of lengths 185,000 and 15,000, respectively. The time of the interval T depends on the geometry and rotational viscosity \(\gamma _1\), and we set it to such a value that the intensity I(t) saturates and effectively becomes constant. For training, the data is – as usual for neural network training⁴² – linearly scaled onto the [0, 1] interval. While the relative intensities I(t) are already limited to this interval by themselves, we transform elastic constants by the “MinMaxScaler” as \(\widetilde{K_{11}^i}=\left( K_{11}^i – \min _j\left( K_{11}^j\right) \right) /\left( \max _j\left( K_{11}^j\right) -\min _j\left( K_{11}^j\right) \right)\) and analogously for \(\widetilde{K_{33}}\) . The training data set is used to train the weights and biases of a dense sequential neural network so that the difference between the predicted output \({\textbf{Y}}\) and the expected output \({\textbf{T}}\) iteratively becomes as small as possible. To quantify the difference, the mean absolute error was used as a loss function. The validation set is used to test the neural network’s performance for data that was not used for training. Training via Adam optimization algorithm⁴³ was performed using Tensorflow Keras software^44,45 with batch size 25 and learning rate \(\eta =0.0003\). A neural network with an input layer of 500 neurons and four hidden layers of 500, 400, 250, and 100 neurons and rectified linear unit (ReLU) activation functions and an output layer of two neurons with sigmoid activation functions was used. It turns out that the network’s architecture can be substantially varied as long as the total number of parameters (connections between neurons) is large enough (for our study, larger than order-of-magnitude \(\sim 10^{5}\)). As observed, adding more layers or increasing the number of neurons in layers in the described model does not significantly improve its accuracy.

A neural network can be trained to recognize both elastic constants \(K_{11}\) and \(K_{33}\) from time-dependent transmittances I(t) through samples with \(\textbf{n}(z,t)=(\cos \theta (z,t),0,\sin \theta (z,t))\) director geometry only if the planar (\(\theta (z=0)=0\)) and the homeotropic (\(\theta (z=D)=\pi /2\)) anchoring are used at the opposite cell surfaces. Planar-planar anchoring geometry allows us to determine \(K_{11}\) only, while in cells with homeotropic anchoring on both surfaces, only \(K_{33}\) can be determined. This is shown in Fig. 3, where 2D histograms show the number count of points in particular sections of the 2D space of true versus predicted constants for examples from the validation set. For a single type of anchoring on both plates, the equilibrium orientation of molecules is always constant along z-axis, either \(\theta (z,t\rightarrow \infty )=0\) or \(\theta (z,t\rightarrow \infty )=\pi /2\), regardless of the elastic constants. Consequently, the transmittance after relaxation, \(I(t\rightarrow \infty )\), does not depend on them and therefore the information about elastic constants is efffectively embedded in the effective dynamics of the time-dependent I(t). In the planar-planar geometry, the dynamics in the last part of the relaxation, when the deformations are small, and the director is only slightly different from its equilibrium configuration \(\textbf{n}(z,t\rightarrow \infty )=(1,0,0)\), is mainly governed by \(K_{11}\). This follows from the fact that in a small-deformations regime, \((\partial n_z/\partial z)^2\gg (\partial n_x/\partial z)^2\) and therefore the Frank-Oseen elastic free energy density (“Methods”: Eq. 1) simplifies to \(f_{FO}\approx K_{11/2}(\partial n_z/\partial z)^2\). This is the reason why only \(K_{11}\) can be determined in such geometry. In homeotropic-homeotropic geometry, it is just the opposite. The dynamics is governed by \(K_{33}\) and therefore only \(K_{33}\) can be determined. However, the equilibrium configuration of the director in a cell with planar-homeotropic anchoring geometry is described by \(\theta (z,t\rightarrow \infty )=\pi z/2D\) if \(K_{11}=K_{33}\) and by a convex or a concave function \(\theta (z,t\rightarrow \infty )\) when \(K_{11}>K_{33}\) or \(K_{11}<K_{33}\), respectively. This results also in the dependence of \(I(t\rightarrow \infty )\) on both elastic constants and for this reason, I(t) carries more information in such geometry and it is possible to determine both elastic constants at the same time. Consequently, such geometry has been chosen to be studied in more detail. As depicted by the red curve in Fig. 6, the mean absolute error of the predicted elastic constants can decrease down to \(\overline{\sigma _{ K_{ii}}}\approx 0.5 \text { pN}\) after 60 epochs of training and batch size 25 using the training set corresponding to the parameters, described in the caption of Fig. 4.

Of course, the prediction depends on the accuracy of other parameters (refractive indices, \(n_o\), \(n_e\), thickness of the cell, D, rotational viscosity, \(\gamma _1\)) that should therefore be known as precisely as possible. For example, the cell thickness and the refractive indices explicitly determine the phase difference between ordinary and extraordinary polarization and thereby the magnitude of the light intensity, whereas the rotational viscosity directly scales with the elastic constants and in turn with the I(t) curves in time. In Fig. 4, we show the distributions of the elastic constants determined by the neural network that was trained on the data set corresponding to \(D= 15.0 \ \upmu \text {m}\), \(n_o=1.523\), \(n_e=1.744\), \(\gamma _1=0.20\ \text {Pa s}\) and tested for time-dependent transmittances that were simulated in systems with elastic constants \(K_{11}=11.0 \text { pN}\), \(K_{33}=17.0 \text { pN}\), but slightly modified parameters D, \(n_o\), \(n_e\), \(\gamma _1\). Training a neural network with data corresponding to the inaccurate thickness of the cell D or refracitve indices \(n_o\), \(n_e\) causes the predicted elastic constants to be in the wrong ratio, while the inaccuracy of the rotational viscosity \(\gamma _1\) causes that both predicted constants are shifted by the same factor.

Once the neural network is well-trained from numerically generated data pairs (\({\textbf{X}},{\textbf{T}}\)) calculated with the same thickness of the cell D, refractive indices \(n_o\), \(n_e\) and the rotational viscosity \(\gamma _1\) as used in the experimental setup, the trained network can be utilized to predict elastic constants of a real sample from an experimentally measured I(t) as well.

Method validation: prediction of \(K_{11}\) and \(K_{33}\) from experimental data

Experimentally, the nematic liquid crystal samples of 5CB and E7 were used in \(D=10.0\ \upmu \text {m}\) and \(D=15.0\ \upmu \text {m}\) thick cells, respectively, with strong uniform planar anchoring at the bottom (\(z=0\)) and homeotropic anchoring at the top surface (\(z=D\)) in both cases. The sample of the NLC between crossed polarizers was illuminated by a lamp with a known spectrum and the transmitted light was measured. First, using transparent electrodes positioned at the boundaries of the cell, the director was randomly deformed within the x–z plane by the pulses of the electric field at times \(t<0\) and after the voltage was turned off at \(t=0\), the director freely relaxed from its randomly deformed initial state to the minimum free energy state. 30 experimental measurements I(t) were done for each material and then normalized to the interval [0, 1] and interpolated to 500 time steps within the same time domain as numerically simulated intensities that were used for training, \([0,T=499\Delta t]\). Since the performance of our method does not depend on the initial state of the director – due to the variety of the director initial states used for the training data set – the transmittance I(t) could be measured after the sample was already partially relaxed. This allows us that the measured intensities I(t) can be interpolated to many random intervals \([\Lambda ,\Lambda +T]\), where \(\text {max}(\Lambda )=T/10\), to increase the number of inputs \(\textbf{X}\) and consequently get more results \({\textbf{Y}}=[K_{11}\), \(K_{33}]\) to achieve better statistics of the predicted values.

The measured time-dependent intensities and distributions of predicted elastic constants for the 5CB and the E7 liquid crystal samples at room temperature (approx. 23 °C) are shown in Fig. 5. In each histogram, there are five different distributions of predictions for \(K_{11}\) and \(K_{33}\) made by five differently trained neural networks. The weights and biases of the networks are always initially set to random values⁴² and since the number of them is very large, their final values are different after each training, but the predictions are generally still similar. The transposed matrices of the values of weights (of size \(100\times 2\)) between the last hidden layer (of size 100) and the output layer (of size 2) from the networks predicting elastic constants of the E7 sample are shown in the panel (g) of Fig. 5. The elastic constants determined by neural networks are compared to the values published in the literature in Table 1. As illustrated in Fig. 4, the accuracy of the predicted elastic constants depends on the accuracy of different material and geometrical parameters, such as D, \(n_o\), \(n_e\), and \(\gamma _1\). In the used (experimental) setup, the thickness of the cell was measured to \(\pm 0.1\, \upmu\)m accuracy, the refractive indices to \(\pm 0.001\) accuracy, in agreement with typical values reported in the literature^46,47,48. Rotational viscosity \(\gamma _1\) was taken from the literature^37,38,49.

By increasing the number of training epochs, the mean absolute error of predicted constants for numerical examples from the validation set decreases, as shown by the red curve in Fig. 6, but the predictions from experimental measurements of the transmittance I(t) become increasingly inaccurate, compared to the expected values⁵¹, \(K_{11}\approx 11.1\ \text {pN}\), \(K_{33}\approx 17.1\ \text {pN}\), as shown by blue and green histograms in Fig. 6. Besides that, multiple peaks can emerge in the distributions of elastic constants determined by neural networks. We speculate that the reason is the following. Differences between the measured and the simulated time-dependent intensities are inevitable: There are always inaccuracies in the parameters (D, \(n_o\), \(n_e\), \(\gamma _1\)) that are used to generate the training data, the backflow and the light scattering are neglected in simulations, there is also some noise in the experimental measurements. Therefore, we can only approximately simulate the transmittances I(t). Besides that, it is known that neural networks can become overfitted^42,53 and therefore less general and more sensitive to details and noise during training with the increasing number of training epochs. We observe that as long as the model is general enough (after less than \(\sim 10\) epochs), the distributions of the predicted constants are broader but centered near the actual constants, while after many (\(\gtrsim 20\)) epochs, the distributions of the predictions get narrower but less accurate. For this reason, to achieve optimal predictions, one should not over-train the network when the aim is to predict material parameters from experimentally measured signals that are, in details, always slightly different from the simulated ones. Therefore, to avoid overfitting, we stopped training after 10 or even fewer epochs of training. In Fig. 5 and in Table 1, the predicted elastic constants were determined by models that were trained in 5 epochs.

Determining \(K_{11}\), \(K_{22}\) and \(K_{33}\)

To extract all three elastic constants \(K_{11}\), \(K_{22}\) and \(K_{33}\) from I(t), it proves that more complex deformations of the director that include twist are needed to emerge in the measuring cell or geometry. To keep the cell geometry, therefore, the director parametrization \(\textbf{n}(z,t)=(n_x(z,t),n_y(z,t),n_z(z,t))=(\cos \phi (z,t)\cos \theta (z,t),\sin \phi (z,t)\cos \theta (z,t),\sin \theta (z,t))\) with infinitely strong anchoring that gives fixed boundary conditions \(\phi (z=0,t)\), \(\theta (z=0,t)\), \(\phi (z=D,t)\), \(\theta (z=D,t)\) is assumed. The initial state of the director in numerical simulations is then determined by two random functions \(\theta (z,t=0)\) and \(\phi (z,t=0)\). The transmittance is again simulated similarly using the Jones formalism, and the training of neural networks is almost the same as for the determination of two constants except for the output vector that is, in this case, 3-dimensional (for the three elastic constants), \({\textbf{T}}_{i}=[K_{11}^i,K_{22}^i,K_{33}^i]\).

As it is shown in Fig. 7, we have found out that in hybrid aligned nematic (HAN) cells (panel (a)) and in twisted nematic (TN) cells (panel (b)), determination of \(K_{11}\), \(K_{22}\) and \(K_{33}\) at once is not possible, while in a tilted twist nematic (TTN) cell with planar (\(\theta (z=0)=0\), \(\phi (z=0)=0\)) and tilted (\(\theta (z=D)=\pi /3\), \(\phi (z=D)=\pi /2\)) anchoring (panel (c)), the I(t) apparently carries information about all three constants, that can therefore be determined by a properly trained neural network. The experimental realization of such cells is beyond the scope of this work, but is clearly realizable⁵⁴.

Determination of initial director configurations

As an alternative use of our method, if all material parameters — including elastic constants — are known, our neural networks methodology can also be used to predict the initial director field \(\textbf{n}(z,t=0)\), from the time-dependent transmittance I(t). Below, we show results for effective 2D director geometry with \(\textbf{n}(z,t)=(n_x(z,t),0,n_z(z,t))=(\cos \theta (z,t),0,\sin \theta (z,t))\).

Using a data set of 390,000 pairs \((\textbf{X}, \textbf{T})=\left( [I



Source link