Liquid crystal droplets under extreme confinement probed by a multiscale simulation approach

ABSTRACT In this work, we computationally investigate liquid crystal (LC) droplets in the size range 0.03–1 μm, confined within shells of combined anchoring conditions. Two different types of surface were defined to promote homeotropic and planar degenerate anchoring, respectively. We identified the LC behaviour within the nanoscale droplets using a bespoke multiscale simulation approach. To study 30 nm droplets, we used coarse grained simulations within the dissipative particle dynamics formalism; to study 0.1 μm and larger droplets, we used a finite element method based on the Landau–de Gennes theory. Good agreement between the two methods was observed in our prior analysis and was confirmed in the present work. We explicitly study droplets of size 0.1 and 1 μm by using continuum mechanics calculations. Our results for the largest droplet are consistent with those available in the literature, suggesting that the extension to smaller droplets presented here is realistic, and therefore can be helpful for innovations in which device intensification could be achieved using LC nanodroplets. Graphical abstract


Introduction
Electro-optic and temperature-dependent properties of liquid crystals (LCs) have been exploited for a long time in applications that range between traditional displays, high-resolution devices for communications, microwave and terahertz devices to biological and chemical sensors [1][2][3][4][5][6][7][8]. Either in bulk, or geometrically confined, LC mesogens exhibit interesting chemical and physical properties that could lead to several innovations in devices and processes. For future developments in sensor applications, it is observed that controlling the LCs confinement in the form of droplets is one of the major challenges [9,10]. Due to the nature and strength of the confinement, LC droplets and nanodroplets should be treated differently than their bulk counterparts. Conclusions achieved by investigating LC mesogens, for example, in a typical LC display cell, most probably will not apply in a spherical geometry where, for example, inevitable defects, such as boojums, occur in planar anchoring conditions. In addition, the smaller the droplets, the fuzzier the distinction between surface and bulk region within the droplets themselves. This is the challenge faced by the development of applications such as polymer-dispersed liquid crystals [11][12][13] as well as others, in which LC droplets are controlled by electric fields, and in sensors applications where LC droplets react to temperature and chemical stimuli. Reading of highly sensitive biosensors, for example, depends on the quantification of how LC droplets and nanodroplets respond to other chemicals [14][15][16][17].
As LC droplets are used in many applications including displays and sensors, it is necessary to understand well the behaviour of LCs under these confinement conditions [18][19][20][21]. Various experimental studies have identified the sensitivity level for different biological and chemical compounds [22,23], in terms of nematic-to-isotropic transition. Because the small changes in LC droplets' internal structures, including the appearance of defects, can at present only be captured by computational studies for nanometre-size systems [24][25][26], we previously studied the interactions of surfactants and nanoparticles with nanometre-size droplets [27,28]. However, atomistic and coarse-grained simulations quickly become computationally prohibitive when the LC droplet size increases above 40-50 nm. At much larger length scales (e.g. mm), continuum mechanics calculations yield results which positively correlate with experiments [29][30][31]. The mesoscale region in between remains difficult to explore computationally, although we have recently demonstrated a phenomenological correspondence between coarsegrained simulations and continuum studies which could prove useful to develop accurate multiscale studies for LC droplets [32].
As shown by previous works in the field, LC modelling either qualitatively supports selected experimental finding at the large scale using the Q-tensor calculations or provides molecular insights related to experimental observations via the implementation of atomistic or Gay-Berne coarse-grained simulations [33][34][35][36][37][38][39]. The few available reports that use both Q-tensor and molecular simulations focus on bulk behaviour or, at most, cylindrical confinement [40][41][42][43][44]. Besides our recent contribution [32], a comprehensive multiscale analysis of LC droplets that combines and contrasts molecular simulations with Q-tensor theory results is missing. The aim of this study is to expand on our initial findings to fill this knowledge gap. The results could provide an initial roadmap for the interpretation of experimental data obtained for LC droplets and nanodroplets.
In this work, we investigated spherical LC droplets under extreme confinement where both strong planar degenerate and strong homeotropic anchoring were applied on the droplet surface with different distributions. We conducted simulations via the coarse-grained dissipative particle dynamics (DPD) formalism, where we identified the molecular behaviour and mapped defect locations. These simulations are complemented by continuum calculations conducted within the Q-tensor approach, via Landau-de Gennes theory. The combination of the two approaches allowed us to map the defects generated by surface alteration and quantify the system energy as a function of the constraints imposed on the LC droplet surfaces. We identified how elastic (splay, twist, bend) and thermotropic energy densities change in response to confinement. We chose a very small size range to compare findings achieved by simulations (30 nm diameter) and by Q-tensor calculations (0.1 μm diameter). We kept the anchoring strong, so interfacial interactions dominate the ordering within the droplet. We compared some of our findings with those obtained for droplets of 1 μm diameter.
The remaining of the paper is organised as follows: in the Methods and algorithms section we describe the algorithms implemented, the molecular models and parameters. In the Results section we discuss orientational order as obtained by the DPD calculations, and we compare it against the correspondent Q-tensor calculations. For each droplet, we discuss energy profiles as well as elastic, thermotropic and anchoring contributions to the overall system energy for each droplet. We conclude proposing possible applications and a future roadmap.

Spherical confinement
The droplet shells used in our DPD simulations are shown in Figure 1. The shells are differentiated based on their capability to enforce anchoring. In Figure 1(a), strong homeotropic anchoring occupies regions of decreasing size centred on the shell equator (homeotropic centre, HC). In Figure 1(b), strong planar degenerate anchoring occupies regions of decreasing size, also centred on the shell equator (planar centre, PC). In Figure 1(c), shells enforce planar degenerate and homeotropic anchoring from the two poles (two-sided droplets, TS).

Liquid crystal models
For the DPD simulations, each of the shells shown in Figure 1 contains 20,000 beads arranged in a hollow sphere formation. Inside each shell, 5000 inert beads and 16,000 rods were inserted. Figure 2(a) shows a representative snapshot (PC-type droplet with 40% planar degenerate anchoring) for the initial configurations used in the DPD simulations. The tips of the rod-like structures have a different affinity towards the surrounding shell compared to the core to achieve the desired anchoring. A harmonic spring is defined between consecutive LC beads in a mesogen to ensure its molecular rigidity [45]. Bond length E bond ¼ k bond � r À r 0 ð Þ 2 and bond bending potentials E angle ¼ k angle � θ À θ 0 ð Þ 2 are used, where r 0 is 0.6 and θ 0 is 180°. Both k bond and k angle equal 100 k B T/r 2 . For all other DPD parameters, we refer to our previous work [32].

Simulation algorithms
The coarse-grained simulations were conducted using the DPD framework [46], as implemented in the LAMMPS software package [47]. These simulations were conducted for ~0.30 μs, with the last 0.15 μs used for data analysis. Because the 1000 data points collected in the last 0.15 μs showed that enthalpy and the potential energy of the systems did not change, we considered that our DPD simulations had reached equilibrium. To conduct the continuum-level calculations, we utilised a method developed by Fernandez and coworkers [48,49]. This method minimises the total free energy in order to find the Q-tensor distribution within a given volume, based on the Landau-de Gennes theory. Continuum elastic theories are computationally less expensive compared to molecular simulations because they consider local average properties of LC molecules. To execute Q-tensor simulations, we built mesh structures with surface properties mimicking the shells shown in Figure 1. We defined two surface  conditions that impose strong homeotropic anchoring and strong planar degenerate anchoring, respectively. A representative snapshot of a mesh structure can be seen in Figure 2(b). The Q-tensor calculations are iterated until the largest change in Q-tensor component values was below five significant digits, which satisfied the minimum energy criteria as in our previous work [32]. To compare the two simulation approaches on the same scale, we built with a spherical droplet of 0.1 μm size, although Q-tensor calculations could be used for systems as large as tens of micrometres. We then picked three different conditions of anchoring and analysed the results obtained for 1 μm droplets with the same anchoring conditions as shown in Figure 1.

Quantification of orientational order in DPD simulations
From the definition of the uniaxial order tensor given by [50,51], a global counterpart appropriate for DPD calculations that considers an ensemble of molecular positions and orientations was proposed in [52]: In Equation (1), N is the number of mesogens in the ensemble and û i is the unit vector representing the molecule's direction. By analogy with the scalar uniaxial order parameter associated to Q-tensor [50,51], a cumulative uniaxial order parameter, S U , can be extracted from the dominant eigenvalue, λ 1 , of Q d , using the relation Perfect axial alignment of LC mesogens would yield an S U value equal to one. However, for bipolar droplets, boojums will prevent reaching such value, although results will be very close to unity. The S U value reflects how smooth the bipolar alignment is in a droplet. Similarly, a radial global order parameter, S R , which quantifies the degree of achievement of a perfect radial order inside an LC droplet, can be defined as [52] In Equation (3), r i is the vector represents the location of mesogen i in the LC droplet.

Calculating global order in Q-tensor simulations
The Q-tensor calculations were conducted using finite element meshes of about 90,000 tetrahedra with adaptive refinement; the programs were run until a tolerance of 10 −5 was satisfied for all Q-tensor components. The LC modelling program used for the Q-tensor calculations was also employed for post-processing the results and extract the global order parameter. Colourcoded isosurfaces from Q-tensor simulations represent the local scalar order parameters within the droplet volume [32,50,51]. To obtain the visuals for defect structures, we adjusted the visible volume for S = 0.3 and lower. The global order was calculated as the volume average of the scalar uniaxial order parameter distribution over the droplet. To make this calculation easier, after the results for the Q-tensor were obtained, they were interpolated onto a coarser regular cubical grid. Since in the regular grid the cubes have the same volume, the global order parameter is obtained as the average of the values in each of the cubes contained by the droplet. For validation purposes, different grid sizes were used. It was found that regular grids with 30 divisions along the droplet diameter were sufficient to obtain results that do not vary by more than 0.25% with respect to finer grids; based on this analysis, regular grids of 30 � 30 � 30 were used to calculate the averages reported in this manuscript.

Visualisation of molecular contribution to order parameter determination
By definition of the order tensor in Equations (1), its calculation and that of SU need all molecules in the system or in a sufficiently large group of molecules. However, for visualisation purposes, once the SU for the system was determined, we post-processed the contribution of each molecule to the dominant eigenvalue. In contrast, radial parameters, as they are defined relative to a known direction, can be defined and calculated for each individual molecule in the system. Colour-coding protocol for SR was applied by using Equation.
(3) for every molecule and without taking the system average. For that reason, while SU colour-coding scheme represents a contribution between 0 and 1, colour-coding scheme for SR maps the actual scalar order of every molecule between -0.5 and 1.

Orientational order analysis
We describe the results starting from droplets on which homeotropic anchoring is imposed on the surface equator, and planar degenerate anchoring at the poles (i.e. Figure 1(a)). For ease, these droplets are referred to as HC. In our previous study, we showed that DPD simulations do not yield a clear picture of defect structures, but they clearly map the defect regions [32]. Due to the discrete nature of individual molecules in a smaller scale, the tendency to disorder is much higher in molecular simulations than in continuum mechanics calculations, where an adaptive mesh method provides refined defect structures. Therefore, the ratio of defect volume to the entire droplet volume is higher in DPD simulations, yielding higher range for orientational order analysis. Building on those observations, we identified the orientational order of the droplets with respect to the uniaxial and radial order parameters. We also found that the energy profiles predicted by DPD align with those from Q-tensor calculations.
As shown in Figure 3(a), DPD simulations reveal that S U increases as the portion of the droplet surface on which planar degenerate anchoring is imposed increases; correspondingly, S R decreases. This was an expected result as planar degenerate anchoring more strongly promotes axial ordering; on the contrary, homeotropic anchoring strongly promotes the radial orientation of the mesogens within the LC droplet. The results show that HC droplets with 10% planar degenerate anchoring show high S U and low S R , while HC droplets with 90% planar degenerate anchoring show low S U and high S R . Example snapshots provided in Figure 3  chosen represent HC droplets whose surface anchoring is 50% planar degenerate and 50% homeotropic. To visualise the meaning of the radial order parameter S R , we colour-coded the mesogens individually as described in Methods and algorithms section. A contribution factor, which ranges from 0 to 1, is employed to visualise how each individual mesogen increments, or deflates the cumulative S U order for the whole droplet. Molecules in Figure 3(a) were colour-coded accordingly to ease visualisation. Additional snapshots to further illustrate the LC droplet structures are provided in Figure 4. These snapshots show that the ring defect in the core of the droplet is only seen for 10% HC. As the planar degenerate surface anchoring ratio increases, the ring in the centre disappears, giving rise to two +½ line defects connecting two ring defects that form around the caps. While in a droplet with 100% planar degenerate surface anchoring, two boojums diametrically opposed would appear, in the droplets considered here, two +½ ring defects appear at the boundary between the caps and the rest of the droplet surface, where the boundary anchoring conditions change.
Another rather interesting characteristic of HC-type droplets is the change in curvature of the line-defects within the bulk as the planar degenerate surface anchoring ratio increases. Between 20% and 40% HC, the defect lines were curved outwards, while for larger ratios the defect lines curve inwards. This is the direct consequence of the different surface conditions imposing different alignments on the LC molecules. Strong homeotropic anchoring induces splay deformation in the bulk, and planar degenerate anchoring favours bend deformation from the surface to the core, the curvature of defect lines significantly depends on proportion of these two anchoring types. However, in DPD simulations the splay and bend-type deformations are not distinguishable. In fact, after 70% HC, LC molecules align as if they were not confined in a droplet. Such parallel alignment created a sudden increase in the value of the S U parameter.
The second set of results was obtained for the planar centre (PC) droplets, which are shown in Figure 1(b). We expect that as the ratio for planar degenerate anchoring to the homeotropic anchoring (calculated in terms of surface area) increases, the orientation must be almost perfectly uniaxial due to homeotropic anchoring on the boojums. As can be seen in Figure 3(b), S U results are consistent with this expectation. The S U order parameter assumes the value of ~0.2 for low planar degenerate anchoring ratio, and it jumps to ~0.6 when at least 60% of the droplet surface is forced to assume planar degenerate anchoring. The snapshots in Figure 4 alongside the results in Figure 3(b) show that the +½ rings formed at the boundary between different surfaces contribute to the change in LC behaviour. In bulk, nematic LCs are expected to align parallel to each other, whereas in bipolar droplets the orientation is expected to be degenerate, hence the boojums. Rings formed due to surface alteration, however, reduce the confinement effects, which was also observed in some of the HCtype droplets. Another point observed on both HC and PC-type droplets is that when the surface alteration on the tips is very small (10% HC and 90% PC), the effect is local. In these circumstances, the overall orientation of the LC was not changed, and the defects in bulk were not affected. Also, 10% HC preserved its +½ ring defect in the core, and 90% PC kept the boojums, alongside the local surface defects created by the imposed homeotropic anchoring. The jump in the parameter value occurs because the rings at the interface of two different surface types become smaller, enabling more mesogens to align with the polar axis of the droplet. In addition, our simulations show that decreasing the surface area on which homeotropic anchoring is imposed affects the boojums, and eventually they disappear. The example snapshot provided in Figure 3(b) reveals that already at 60% planar degenerate anchoring, the tips of droplets endorsed overall axial ordering for the most part and that the point defects disappeared. Instead, boojums are composed of two +½ ring defects. Additional snapshots representing all the simulated droplets are provided in Figure 4. Profiles in Figure 3(b), for S R , are similar to that of Figure 3(a). As the homeotropic anchoring decreases in terms of surface ratio, S R decreases for the same reasons discussed for Figure 3(a).
The last data set is for the two-sided (TS) droplets, which are shown in Figure 1(c). The results in Figure 3 (c) complement those shown elsewhere for LC droplets with planar degenerate anchoring over 5%, 10%, 20%, 80%, 90% and 95% of their surface [32]. The results in Figure 3(c) show that the S R order parameter obtained from DPD simulations decreases as the surface ratio of planar degenerate anchoring increases. We note that the results for the S U order parameter are somewhat unexpected: they increase regularly from 0.23 to 0.65 as the planar degenerate anchoring increases from 10% to 40%, but then they become erratic in the range 0.3-0.8 at higher values of the planar degenerate anchoring. Analysis of our results shows that local mesogen ordering is responsible for these effects: as we force the mesogens in a controlled volume to align with an imposed boundary condition (i.e. homeotropic anchoring), the locally achieved alignment may or may not be consistent with the global order of the droplet, leading to the observed S U behaviour. Snapshots of TS-type droplets provided in Figure 4 show that TS droplets possess defect structures similar to HC droplets when the ratio of planar degenerate surface anchoring is low; when the ratio is high, they became more PC-like. It is expected that due to physical similarities between these droplets, 10-40% TS droplets resemble HC, and 50% and higher resemble PC. Both quantitative and qualitative properties of TS droplets shift from HC to PC with respect to increasing planar degenerate surface anchoring ratio.
Q-tensor calculations with adaptively refined mesh structures revealed that within that level of analysis, the global order within the droplets (i.e. averaged within the volume of the droplet) did not depend on the choice of surface anchoring. For all cases in Figure 3(a-c), the global order parameter increases from 0.48 to 0.52 as the ratio of planar degenerate surface anchoring increases. This is likely a result of the fact that the volumetric ratio between defects and the whole droplet is small for each of the droplets considered, even when defect types and locations change with the surface anchoring conditions. Although results obtained from Q-tensor calculations are somewhat in line with the results of S U parameter for DPD simulations, the change in order is much less pronounced in the former method, almost negligible. While DPD was effective at mapping the defect regions within a droplet, the Q-tensor approach was essential to elucidate the defects themselves. This synergism between the two computational approaches was harnessed to quantitatively analyse the orientational order parameters. In our DPD simulations, the volumes affected by defects were a larger proportion of the droplet volumes, while in Q-tensor calculations defects were confined in small volumes and did not affect the cumulative order for the droplets. This complementarity yields a difference between order parameter calculations conducted via the two computational approaches [32]. The order parameter, whether radial or uniaxial, is dependent on individual molecules in DPD simulations. On the other hand, the volume average of order in Q-tensor calculations is dependent of a complete volume, as explained in Methods and algorithms section. Additionally, the conditions of homeotropic and planar degenerate anchoring was also set differently in these two algorithms. In DPD simulations, anchoring type and strength was satisfied by the strong repulsion coefficients between different type of beads. In Q-tensor calculations, however, it is defined via different elastic constants on different volumes and boundaries. It is reassuring that although the two methods have fundamentally different approaches, the results show the same tendency towards cumulative order and local defect regions. Perhaps, the implications of these phenomena could be better analysed and appreciated with the help of advanced experiments. Nevertheless, we stress that both DPD simulations and Q-tensor calculations yield results that are visually similar in terms of the localisation of the defect regions within the LC droplets.
We point out that DPD simulations also yield the total energy and the enthalpy for each simulated LC droplet. As shown in Figure 3(d), all droplet types (HC, PC and TS) yield the same qualitative results: both the total droplet energy and the enthalpy increase monotonically as the planar degenerate anchoring surface ratio increases. The difference in slope between the two curves suggests that the pressure within the LC droplet also increases with the planar degenerate anchoring surface ratio. The increase in pressure could only be observed for DPD simulations as it was not included in our Q-tensor calculations. Such change in pressure was somewhat expected. Because increase in pressure as planar degenerate surface ratio increases means that LCs (particularly LC tips) are exposed to highly repulsive surface mostly, and less to attractive. As a result, within an enclosed boundary conditions, the same amount of LC molecules are exposed to higher pressure, which is reflective in the change in increase between enthalpy and total energy of droplets in DPD simulations.

Energy comparison and energy density analysis
The Q-tensor analysis allows us to extract the elastic and thermotropic energy contributions to the total energy of the LC droplets considered. The results are shown in Figure 5. Figure 5(a) shows the elastic energy density of all LC droplets simulated. For all cases, the elastic energy increases as the planar degenerate surface anchoring ratio increases, regardless of the LC droplet type. For the PC droplets, the elastic energy increase flattens as fully planar degenerate anchoring on the surface is approached; the opposite holds for HC droplets. For TS droplets, on the other hand, the results show similarity with HC droplets at one end of the graph, and with PC droplets at the other end, due to the structural similarities of the shells of the three droplet types. In fact, TS droplets have similar morphological properties than the HC droplets when the planar degenerate surface area ratio is low, and similar to the PC droplets when the planar degenerate surface area ratio is high. Similar results were obtained for the thermotropic energy density profiles, as shown in Figure 5(b).
Results so far suggest that the relation between DPD and Q-tensor simulations is qualitative for defect structures within the droplets, and quantitative in terms of total energy. The same conclusion was reached in our previous work, which may validate findings presented in this paper [32]. Motivated by such correspondence, we fitted the energy densities of DPD and Q-tensor calculations in Figure 5(c), obtaining the correlation y ¼ 2:22 � 10 9 x À 1:11 � 10 6 , where x represents the energy value in reduced DPD units. We find that comparing energy densities rather than total energy values is more accurate as it allows us to take into consideration the different sizes of the droplets considered in DPD simulations (30 nm) and Q-tensor calculations (0.1 μm).
From the results in Figure 5(c), we observe that DPD simulations yield similar results for all LC droplet types considered (HC, PC and TS), whereas Q-tensor calculations predict that PC droplets have higher energy and HC ones have lower energy at a given surface area ratio with planar degenerate anchoring. These results suggest that in DPD simulations changes in surface morphology do not reflect on energetics, consistent with the results shown in Figure 3(d). Our results show that energy values predicted by the DPD simulations are lower than those predicted by the Q-tensor approach for PC droplets, higher for HC droplets and for the most part comparable for TS droplets. The data sets from both approaches are provided in Table 1.

Effect of droplet size on Q-tensor calculations
While it would be very computationally demanding to simulate LC droplets of diameter larger than ~30 nm using DPD simulations, our preliminary analysis suggests that Q-tensor calculations could yield large errors for droplets smaller than 0.1 μm. Because these size limits do not overlap, we did not discuss the effect of droplet size on the shape of defects and on the droplet energy. However, using the Q-tensor calculations we can quantify droplet size effects for LC droplets larger than 0.1 μm. For three LC droplets (60% HC, 60% PC and 60% TS, respectively), we increased the diameter to 1 μm while maintaining similar surface properties.   Figure 6(a) shows the total energy densities of droplets of diameter 0.1 and 1 μm. The energy density of 0.1 μm droplets is almost one order of magnitude higher than that for the 1 μm droplets. For both diameters, PC droplets have the highest energy density (1.40 × 10 4 and 9.44 × 10 4 J/m 3 ), followed by TS (1.13 × 10 4 and 7.63 × 10 4 J/m 3 ) and HC droplets (7.42 × 10 3 and 5.82 × 10 4 J/m 3 ). In Figure 6(b), HC droplets of diameter 0.1 and 1 μm with 60% planar degenerate anchoring on the surface are presented. Note that these simulations were conducted only at the Q-tensor level because of computing-power limitations. The major difference between the two droplets is the alignment on the surface that imposes planar degenerate anchoring. Inside the 0.1 μm droplet, there is a layer of mesogens aligned parallel to the planar degenerate anchoring interface. For the 1 μm droplet, however, only a thin surface layer of mesogens is aligned parallel to the surface, with the rest of the mesogens perpendicular to the surface; the latter configuration yields a radial alignment at the core of the droplet. When large, the direction within droplet gets diffused and defects are less pronounced. In the small droplets the higher degree of confinement restricts LCs from freely moving, and high distortion occurs.
Releasing the degree of LCs confinement yield a change from the more defined defect structures, which is a finding consistent with literature [20,21]. Because of the strong surface anchoring, the equilibrium order distribution within droplets of different diameters changes due to different surface-to-volume ratios. For the 0.1 μm droplet, the strong anchoring dominated the alignment within the droplet while for the 1 μm droplet the bulk energy overcame the anchoring conditions and limited their effect to a thin interfacial layer.
The discussed difference caused by surface-to-volume ratio is visible also in Figure 6(c), where a PC-type droplet with 60% planar degenerate surface anchoring is shown. The effect on the global order due to the planar degenerate anchoring becomes less effective for the larger droplets, and the LCs mesogens form two +½ defects at the centre of the anchoring surface, caused by the homeotropic alignment forced by the bulk of the droplet. A similar behaviour is observed in Figure 6(d) for TS droplets. Different from 0.1 μm size droplets, we visually detect twisted structures in the 1 μm droplets, although these structures are not quantitatively identified. We presume that, because the surface energy is less dominant as the droplet volume increases, the bulk energy dictates the LC mesogens alignment, yielding twisted structures in the LC droplets of diameter 1 μm.
It is reassuring that the change in orientation due to the droplet size has been reported previously, although our results differ somewhat from those reported in the literature. Gupta et al. reported that for the same interfacial chemistry around the droplets the change in ordering as size decreases occurs from bipolar (8 μm), to pre-radial (1 μm), and radial (0.7 μm) [53]. In another work, Miller and Abbott reported that in a system with multiple droplets the number of droplets that exhibit radial configuration sharply decreases as the size of these droplets increases from 2 μm to 10 μm [16]. It should be noted that we refer to experimental literature results, in which surface anchoring was modified via the addition of surface-active chemicals. According to our Q-tensor calculations with different surface conditions presented in Figure 6, increasing the droplet diameter from 0.1 μm to 1 μm yield radial alignment inside the droplet, a change towards a predominantly radial configuration was observed as the droplet size increased. What we observed for 1 μm droplets is in line with experimental observations for 6.5 μm droplets [16]. We however highlight that, because of the different scales being probed, the increase in homeotropic alignment observed in our work corresponds to the decrease observed experimentally upon an increase in LC droplet size. The overall differences are due to the contribution of the surface energy, which becomes dominant when the droplet size reduces to 0.1 μm and smaller. For droplets in this size range, the surface energy dictates the alignment within the LC droplet, a prediction which should be verified experimentally, once suitable tools are designed to conduct defect analysis within nano-scale droplets.
Overall, this study showed that droplets under extreme confinement yield complex results, which differ compared to results observed in micrometre-sized droplets. This is a significant step towards understanding nanodroplets for intensified devices.

Conclusions
In this work, we compared synergistic simulation techniques to elucidate the properties of LC droplets in the meso-scale range. DPD simulations were used to quantify molecular phenomena for 30 nm droplets, whereas Q-tensor calculations provided insights for droplets of size 0.1 μm and above. The two methodologies are complementary, with DPD providing molecular-level insights and Q-tensor defect formation and morphology. Because the LC droplets considered in the Q-tensor calculations are rather large compared to the size of one LC mesogen, the volume of the defects is small compared to that of the LC droplet, and the order parameter shows that all LC droplet considered here of diameter at least 0.1 μm maintain nematic order for all surface anchoring conditions chosen. Because DPD simulations interrogate much smaller LC droplets, whose size is comparable to the size of a few LC mesogens, the order parameters obtained strongly depend on the anchoring conditions. While the defects are rather confined in the results of our Q-tensor calculations, they occupy broad volumes within the LC droplets simulated with the DPD formalism.
The profiles for the total LC droplet energy provide the foundation for a meso-scale model that allows us to link DPD with Q-tensor calculations. The model results depend on surface anchoring, especially for DPD simulations. Our analysis revealed that DPD cannot differentiate the morphological difference of the shells with the same planar degenerate anchoring ratio, in terms of the total energy calculations.
Finally, once the correlation between nano-and mesoscale LC droplets was established, we used the Q-tensor approach to compare defects and energies of six LC droplets whose diameter was increased from 0.1 μm to 1 μm. The results revealed that the energy density in 0.1 μm droplet is one order of magnitude higher than that in 1 μm droplets, which shows how surface energy becomes dominant, as opposed to bulk energy, as the surface-tovolume ratio increases for LC droplets. Because our results for 1 μm LC droplets are consistent with literature reports, the approach presented here could offer a reliable methodology for sampling the properties of LC nanodroplets, which is expected to be important particularly for device intensification.

Acknowledgments
Funding for this research was provided, in part, by the Department of Chemical Engineering at UCL. The DPD simulations presented were conducted at the University College London Research Computing Platforms Support (Grace, Myriad). ZS would like to thank Dezhi Shen for fruitful discussions on finite-element modelling and support in the implementation of the related software. AS acknowledges funding from EPSRC via grant number EP/ T004282/1.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
This work was supported by the Engineering and Physical Sciences Research Council [EP/T004282/1].