1 Introduction

Nematic liquid crystals (NLCs) are classical examples of partially ordered materials or viscoelastic anisotropic materials, with long-range orientational ordering de Gennes and Prost (1995). The NLC molecules are typically asymmetric in shape, e.g. rod-like, and they tend to align along locally preferred directions, referred to as nematic directors in the literature de Gennes and Prost (1995). The optical, mechanical and rheological NLC responses are direction-dependent, with distinguished responses along the nematic directors. Indeed, the anisotropic NLC responses to external stimuli, interfaces and boundaries make them soft, self-organising and the cornerstone of several NLC applications in science and engineering Phillips and Rey (2011). Nematics in confinement have attracted substantial scientific interest in the academic and industrial sectors Lagerwall and Scalia (2012). In fact, nematics inside a planar cell are the building block of the celebrated twisted nematic display, and contemporary work has focused on the tremendous potential of NLCs for multistable systems, i.e. NLC systems with multiple observable states without applied fields, such that each observable state offers a distinct mode of functionality Spencer et al. (2010). Nematic defects play a key role in multistability, where a nematic defect is a point/line/surface where the nematic director cannot be uniquely defined de Gennes and Prost (1995). Nematic defects often organise the space of stable states in multistable systems in the sense that we can classify the stable states in terms of the nature, multiplicity and locations of defects. Additionally, defects have distinct optical signatures and can act as binding sites or “special” sites in materials design or applications. The delicate interplay between geometric frustration, boundary effects, material properties and defects in multistability leads to a suite of challenging mathematical questions at the interface of applied topology, nonlinear partial differential equations and scientific computation (to name a few). Equally, it gives new inroads into engineered soft materials, topological materials or meta-materials which could find new applications in photonics, robotics and artificial intelligence Geng et al. (2021).

In this paper, we focus on the mathematical modelling and numerical computation of NLC equilibria outside regular two-dimensional (2D) polygons with homeotropic boundary conditions, i.e. the nematic director is constrained to be normal to the polygon boundary. We work in a reduced 2D Landau–de Gennes framework, details of which are given in the next section, and this can capture the nematic directors and the nematic defects in a 2D setting, along with informative insights into how multistability can be tailored by the shape and size of the polygon. This toy mathematical problem models a single colloidal particle, in the shape of a 2D polygon, suspended in an extended NLC medium, which is of both experimental and theoretical interest Muševič and Škarabot (2008); Evans et al. (2011); Yuan et al. (2019); Gupta and Rey (2005); Phillips and Rey (2011). In Muševič and Škarabot (2008); Evans et al. (2011); Yuan et al. (2019), the authors fabricate almost 2D platelets of different polygonal shapes suspended in a NLC medium. Using advanced microfabrication techniques and optical methods, they can track the director profiles and the associated defects. The authors observe multiple types of defects: dipoles, Saturn rings, and various linked or entangled defect loop lines intertwining arrays of embedded inclusions. In fact, in Muševič and Škarabot (2008), the authors use optical tweezers to manipulate the defect lines, to link them or disentangle them, and in doing so, create various exotic knotted defect patterns. In all cases, the observed states and their defect patterns strongly depend on the geometry and orientation of colloidal particle(s) and their boundary effects, offering excellent examples of organic self-assembled structures in NLC media.

In Han et al. (2020), we study the “interior” problem of multistability for NLCs confined to a regular 2D polygon with tangent boundary conditions, i.e. the nematic director is tangent to the polygon edges, in a reduced Landau–de Gennes (LdG) framework. We use conformal mapping techniques and methods from elliptic partial differential equations to study nematic equilibria in two distinguished limits, phrased in terms of a dimensionless variable \(\lambda \). The variable, \(\lambda \), is the ratio of two length scales - the ratio of the polygon edge length to the nematic correlation length, which is a material-dependent length scale related to typical defect sizes. In the \(\lambda \rightarrow 0\) limit, the limiting profile is unique with a single isolated point defect at the centre of the polygon, coined as the Ring solution. The only exceptions are the triangle and the square, which are dealt with separately. In fact, the limiting profile for a triangle has an isolated fractional point defect at the centre whereas the unique limiting profile for a square domain is labelled as the Well Order Reconstruction Solution (WORS), first reported in Kralj and Majumdar (2014). The WORS is a distinctive profile with two orthogonal defect lines along the square diagonals (also see Kralj and Majumdar (2014)). In the \(\lambda \rightarrow \infty \) limit, we use combinatorial arguments to demonstrate multistability, i.e. there are at least \(\left( {\begin{array}{c}K\\ 2\end{array}}\right) \) distinct NLC equilibria on a K-polygon with K edges. In related papers, we comprehensively study NLC solution landscapes on 2D polygons with tangent boundary conditions Han et al. (2021, 2021).

We study the complementary “exterior” problem in this paper: asymptotic and numerical studies of NLC equilibria outside a regular polygonal hole, immersed in \(\mathbb {R}^2\). This problem, and related problems, have received some analytical interest although systematic studies are missing. For example, in Phillips and Rey (2011), the authors study stable NLC profiles outside a square hole with homeotropic boundary conditions. They numerically observe string defects (line defects) pinned to square edges, defects at square vertices and interior point defects, depending on the temperature and square size. In Wang et al. (2017, 2018), the authors numerically investigate the possible structures of NLCs with one, two and multiple spherical inclusions. In Gupta and Rey (2005), the authors model NLCs with multiple spherical inclusions and numerically investigate how the defect set depends on the spatial organisation and properties of the spherical particles, along with those of the ambient NLC media. In Alama et al. (2016), the authors rigorously analyse NLC equilibria outside a spherical particle with homeotropic boundary conditions, in a three-dimensional (3D) LdG framework. They obtain elegant limiting profiles in the small particle and large particle limit, and in fact, produce an analytic expression for the celebrated Saturn ring solution with a distinct defect loop around the spherical particle.

In this paper, we focus on the effects of shape and size of the polygonal hole on the corresponding NLC equilibria, in a reduced LdG model. The methodology follows that of Han et al. (2020), the key difference being the extra degree of freedom rendered by the far-field boundary conditions, away from the polygon boundary. As with the interior problem, we compute limiting profiles for the stable NLC equilibria in the \(\lambda \rightarrow 0\) and \(\lambda \rightarrow \infty \) limits, where \(\lambda \) has the same interpretation as in Han et al. (2020), accompanied by supplementary numerical results. The stable equilibria are modelled by local or global energy minimisers of the reduced LdG free energy, which in turn are solutions of the associated Euler-Lagrange equations that are a system of two coupled nonlinear partial differential equations. In the \(\lambda \rightarrow 0\) limit, there is a unique NLC equilibrium or equivalently, a unique minimiser of the reduced LdG free energy. However, the limiting profiles are more varied compared to the interior problem. For a square hole, we can observe either line defects or point defects at the square vertices, depending on the far-field condition, in the \(\lambda \rightarrow 0\) limit. There is qualitative agreement with the numerical results in Phillips and Rey (2011). In general, the locations of the defects for the unique limiting profile depend on the far-field condition. Using a result from Ginzburg–Landau theory in Bauman et al. (1993), we show that there are exactly two interior point defects for a generic polygonal hole, \(E_K\) with K edges and \(K>4\), and the location of these defects can be tuned with the far-field condition. In the \(\lambda \rightarrow \infty \) limit, we provide a simple estimate for the number of stable NLC equilibria using combinatorial arguments, and multistability is enhanced compared to the interior problem. This is because the exterior problem has lesser symmetry than the interior problem, due to the far-field conditions. For example, on the interior of a square domain, there are two rotationally equivalent diagonal solutions for which the NLC profile is approximately aligned along the square diagonal. We lose this equivalence for the exterior problem since the far-field condition breaks the symmetry.

Applied mathematics focuses on the development of new mathematical methods, and equally elegant applications of known methods to new settings. Our work falls into the second category, where we largely build on previous work, and use techniques from complex analysis, Ginzburg–Landau theory for superconductivity, symmetry results and combinatorial arguments to analyse limiting profiles, complemented by numerical studies to support the theory. In doing so, we demonstrate how geometric frustration and nematic defects go hand in hand for tailored multistability, and this is a good forward step for rigorous mathematical studies of NLC solution landscapes in complex geometries with voids, mixed boundary conditions and in some cases, multiple order parameters Han et al. (2021).

The paper is organised as follows. In Sect. 2, we review the reduced 2D LdG framework for modelling NLCs in 2D scenarios. In Sect. 3, we focus on the \(\lambda \rightarrow 0\) limit of minimisers of the reduced Landau–de Gennes free energy. We use the Schwarz–Christoffel mapping to define an associated boundary-value problem on the unit disc, for each regular polygonal hole, and this boundary-value problem is explicitly solved. The defect set is tracked analytically, along with its dependence on K and the far-field condition. In Sect. 4, we shift focus to the \(\lambda \rightarrow \infty \) limit and the limiting problem is captured by the Laplace problem for an angle in the plane, with Dirichlet boundary conditions. This angle models the 2D nematic director. The Dirichlet boundary conditions for the angle are not uniquely defined, and this leads to multistability in this limit. We present illustrative numerical examples for a square and a hexagon and conclude in Sect. 5 with a summary and some perspectives.

2 Theoretical Framework

The Landau–de Gennes (LdG) theory is a celebrated continuum theory for nematic liquid crystals de Gennes and Prost (1995) and was indeed one of the reasons for awarding the Nobel Prize for physics to Pierre-Gilles de Gennes in 1991. The LdG theory is a phenomenological theory that assumes macroscopic quantities of interest vary slowly on the molecular length scale and is based on the crucial concept of a LdG order parameter, which is a macroscopic measure of the NLC anisotropy or orientational ordering. The LdG order parameter, known as the \(\textbf{Q}\)-tensor order parameter is a symmetric, traceless \(3\times 3\) matrix with five degrees of freedom. The nematic director is often interpreted as the eigenvector of the LdG \(\textbf{Q}\)-tensor with the largest positive eigenvalue Han et al. (2020). Consider a three-dimensional (3D) domain, with a 2D cross section \(\Omega \), and height, h. Then, the 3D LdG energy is given by :

$$\begin{aligned} I[\textbf{Q}]: = \int _{\Omega \times [0,h]} \frac{L}{2}|\nabla \textbf{Q}|^2 + F_b(\textbf{Q})~dV \end{aligned}$$
(2.1)

where \(F_b(\textbf{Q}): = \frac{A}{2}\text {tr}\textbf{Q}^2 - \frac{B}{3} \text {tr}\textbf{Q}^3 + \frac{C}{4} (\text {tr}\textbf{Q}^2)^2\), A is a re-scaled temperature, B and C are positive material-dependent constants Mottram and Newton (2014). The constant \(L>0\) is a material-dependent elastic constant, with units of Newtons. Typical values of L are \(10^{-12}\) Newtons Priestly (2012). We have adopted the Dirichlet elastic energy density, \(w(\nabla \textbf{Q}) \propto |\nabla \textbf{Q}|^2 = \sum _{i,j,k=1}^{3} Q_{ij,k}^2\), where \(Q_{ij,k} = \frac{\partial Q_{ij}}{\partial x_k}\), based on the assumption that all elastic deformations, e.g. splay, twist and bend deformations are energetically degenerate. The physically observable stable equilibria are modelled by minimisers of (2.1) in an appropriately defined admissible space.

For thin 3D systems, for which h is much smaller than the dimensions of \(\Omega \), it suffices to work with the reduced Landau–de Gennes (rLdG) model, based on the assumption that the nematic director is in the cross-section plane, and structural details are invariant across the height of the system Golovaty et al. (2017). The rLdG model has been successfully applied for capturing the qualitative properties of physically relevant solutions and for probing into defect cores Brodin et al. (2010); Gupta and Rey (2005); Muševič and Škarabot (2008); Muševič et al. (2006). In the rLdG model, the nematic state in the 2D cross section/2D domain is described by a reduced order parameter: a symmetric, traceless \(2\times 2\) matrix, \(\textbf{P}\), as given below

$$\begin{aligned} \textbf{P} = \left( \begin{array}{cc} P_{11}&{}P_{12}\\ P_{12}&{}-P_{11}\\ \end{array}\right) . \end{aligned}$$

We work at a special temperature, \(A = -B^2/3C\), where B and C are as before de Gennes and Prost (1995), and with the following rLdG free energy:

$$\begin{aligned} F[\textbf{P}]: = \int _{\Omega }\frac{L}{2}|\nabla \textbf{P}|^2+f_b(\textbf{P}) \textrm{d A}. \end{aligned}$$
(2.2)

where \(\Omega \) is the 2D cross section of our 3D domain. In the remainder of this manuscript, \(\Omega \) is the complement of a regular 2D polygon with K edges, \(E^C_K\), in \(\mathbb {R}^2\). Here, we have employed the Dirichlet elastic energy density as in (2.1) and the bulk energy density as

$$\begin{aligned} f_b{(\textbf{P})} = -\frac{B^2}{4C}\text {tr}\textbf{P}^2 + \frac{C}{4}(\text {tr}\textbf{P}^2)^2. \end{aligned}$$

We choose this temperature, partly for comparison with previous work in Han et al. (2020) and Fang et al. (2020), and partly because for this special temperature, the critical points of the rLdG model exist as critical points of the full 3D LdG free energy in 3D settings, for suitably defined boundary conditions. This is generally not true for arbitrary \(A<0\); see Canevari et al. (2020) for more details. More precisely, at \(A=-B^2/3C\), given a critical point \(\textbf{P}_c\) of (2.2), there exists a critical points \(\textbf{Q}_c\) of (2.1) such that

$$\begin{aligned} \textbf{Q}_c = \textbf{P}_3 - \frac{B}{6C}\left( 2\hat{\textbf{z}}\otimes \hat{\textbf{z}}- \hat{\textbf{x}}\otimes \hat{\textbf{x}}- \hat{\textbf{y}}\otimes \hat{\textbf{y}}\right) , \end{aligned}$$

where \(\hat{\textbf{x}}, \hat{\textbf{y}}. \hat{\textbf{z}}\) are coordinate unit-vectors. The matrix, \(\textbf{P}_3\) is a \(3\times 3\) symmetric traceless matrix, such that \( (\textbf{P}_3)_{ij} = (\textbf{P}_{c})_{ij}\) for \(i,j=1,2\) and all remaining matrix entries are set to zero.

The energy (2.2) is non-dimensionalised with \(\overline{x} = x/\bar{\lambda }\), where \(\bar{\lambda }\) is the edge-length of the polygonal hole.

$$\begin{aligned} F[\textbf{P}]: = \int _{E_K^C}\frac{1}{2}|{\bar{\nabla }} \textbf{P}|^2+\frac{\bar{\lambda }^2}{L}f_b(\textbf{P}) \textrm{d A}. \end{aligned}$$
(2.3)

The working domain is now \(E_K^C\), the complement of a 2D re-scaled polygon \(E_K\) with K edges of unit length, centred at the origin, with vertices

$$\begin{aligned} w_k = \left( \cos \left( 2\pi \left( k-1\right) /K\right) ,\sin \left( 2\pi \left( k-1\right) /K\right) \right) ,\ k = 1,\ldots ,K. \end{aligned}$$

We label the edges counterclockwise as \(C_1,\ldots , C_K\), starting from \(\left( 1,0\right) \). See figure 1. In the following, we drop the bar over \(\nabla \) for brevity.

Fig. 1
figure 1

The normalised domain \(E_K^C\), with \(K = 6\), as an illustrative example

We can also write \(\textbf{P}\) in terms of an order parameter, s, and an angle \(\gamma \) as shown below

$$\begin{aligned} \textbf{P} = 2s\left( \textbf{n}\otimes \textbf{n}-\frac{1}{2}\textbf{I}_2\right) , \end{aligned}$$
(2.4)

where \(\textbf{n} = \left( \cos \gamma ,\sin \gamma \right) ^T\) is the nematic director in the plane and \(\textbf{I}_2\) is the \(2\times 2\) identity matrix, so that

$$\begin{aligned} P_{11} = s\cos \left( 2\gamma \right) ,\ P_{12} = s\sin \left( 2\gamma \right) . \end{aligned}$$

The defect set is simply identified with the nodal set of s, also see Han et al. (2020), interpreted as the set of no planar order in \(E_K^C\).

We impose homeotropic boundary conditions on \(\partial E_K\), which requires \(\textbf{n}\) in (2.4) to be homeotropic/normal to the edges of \(E_K\). However, there is a necessary mismatch at the corners/vertices. We impose a continuous Dirichlet boundary condition, \(\textbf{P} = \textbf{P}_b\), on \(\partial E_K\) as defined below:

$$\begin{aligned} \begin{aligned}&P_{11b}\left( w\right) = \alpha _k = {\left\{ \begin{array}{ll} &{}\frac{B}{2C}\left( (1-\frac{a}{2})\cos \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) \right. \\ &{}\quad \left. + \frac{a}{2}\cos \left( \frac{\left( 2k+1\right) 2\pi }{K}\right) \right) , ||w-w_k||\ge ||w-w_{k+1}||,\\ &{}\frac{B}{2C}\left( (1-\frac{a}{2})\cos \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) \right. \\ &{}\quad \left. +\frac{a}{2}\cos \left( \frac{\left( 2k-3\right) 2\pi }{K}\right) \right) , ||w-w_k||\le ||w-w_{k+1}||, \end{array}\right. }\\&P_{12b}\left( w\right) = \beta _k = {\left\{ \begin{array}{ll} &{}\frac{B}{2C}\left( (1-\frac{a}{2})\sin \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) \right. \\ &{}\quad \left. +\frac{a}{2}\sin \left( \frac{\left( 2k+1\right) 2\pi }{K}\right) \right) , ||w-w_k||\ge ||w-w_{k+1}||,\\ &{}\frac{B}{2C}\left( (1-\frac{a}{2})\sin \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) \right. \\ &{}\quad \left. +\frac{a}{2}\sin \left( \frac{\left( 2k-3\right) 2\pi }{K}\right) \right) , ||w-w_k||\le ||w-w_{k+1}||, \end{array}\right. } \end{aligned} \end{aligned}$$
(2.5)

where \(a(w,\sigma )\in [0,1]\), \(w\in C_k\), \(k= 1,\cdots , K\) is an interpolation function, which is continuous and strictly monotonic about \(||w-(w_k+w_{k+1})/2||\), increasing from 0 at the centre of \(C_k\), to 1 at the two vertices of \(C_k\). Further, we take \(0<\sigma \le 1/2\) and \(a(w,\sigma )\rightarrow 0\) as \(\sigma \rightarrow 0\). At the vertices \(w = w_k\), we set \(\textbf{P}_b\) to be equal to the average of the two constant values on the two intersecting edges, and at the edge mid-point, i.e. \(w = (w_k+w_{k+1})/2\), we have strictly homeotropic conditions. As \(\sigma \rightarrow 0\), the Dirichlet boundary conditions are piece-wise constant and compatible with strict homeotropic conditions

$$\begin{aligned} P_{11b}(w) = \hat{\alpha }_k= & {} \frac{B}{2C}\cos \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) , \nonumber \\ P_{12b}(w) = \hat{\beta }_k= & {} \frac{B}{2C}\sin \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) ,\ w\in C_k. \end{aligned}$$
(2.6)

Further, the domain \(E_K^C\) is unbounded and we impose uniform/constant boundary conditions at infinity as given below:

$$\begin{aligned} \lim _{|w|\rightarrow \infty } \textbf{P}= \textbf{P}^*:= {\frac{B}{C}}\left( \textbf{n}^*\otimes \textbf{n}^* - \frac{1}{2}\textbf{I}_2\right) , \end{aligned}$$
(2.7)

where \(\textbf{n}^* = (\cos (\gamma ^*),\sin (\gamma ^*))\) is the constant nematic director at infinity. To avoid confusion, we reiterate that \(\gamma \) is the director angle and \(\gamma ^*\) is associated with the boundary condition at infinity.

The corresponding Euler–Lagrange equations are

$$\begin{aligned} \begin{aligned} \Delta P_{11}&= \frac{2C\bar{\lambda }^2}{L}\left( P_{11}^2+P_{12}^2-\frac{B^2}{4C^2}\right) P_{11},\\ \Delta P_{12}&= \frac{2C\bar{\lambda }^2}{L}\left( P_{11}^2+P_{12}^2-\frac{B^2}{4C^2}\right) P_{12}. \end{aligned} \end{aligned}$$
(2.8)

For the purposes of this paper, the crucial dimensionless parameter is

$$\begin{aligned} \lambda ^2 = \frac{ 2 C}{L} \bar{\lambda }^2 \end{aligned}$$

where the nematic correlation length, \(\xi _n \propto \sqrt{\frac{L}{C}}\) at the fixed temperature \(A = -\frac{B^2}{3C}\). Recall that L has units of Newton (N) and C has units of \(N m^{-2}\), so that \(\xi _n\) has the units of length.

The admissible space is

$$\begin{aligned}&\mathcal {H}_{\infty }:=\textbf{P}^* + \mathcal {H}, \end{aligned}$$
(2.9)
$$\begin{aligned}&\mathcal {H}:=\left\{ \textbf{H}\in H^1_{loc}(E_K^C;\mathcal {S}_0):\int _{E_K^C}|\nabla \textbf{H}|^2 + \int _{E_K^C}\frac{|\textbf{H}|^2}{|w|^2}<\infty \right\} \end{aligned}$$
(2.10)

where

$$\begin{aligned} \mathcal {S}_0:=\{\textbf{P}\in M_2(\mathbb {R}): P_{ij} = P_{ji},tr(\textbf{P}) = 0\}. \end{aligned}$$

The free energy (2.3) is not always finite, since the potential term \(f_b\ge 0\) may very well not be integrable in \(E_K^C\). However, since \(f_b(\textbf{P}^*) = 0\), we may find a compactly supported \(\textbf{H}\in \mathcal {H}\), for which \(\textbf{P} = \textbf{P}^* + \textbf{H}\) satisfies the boundary condition (2.5). The existence of solution of (2.8) in \(\mathcal {H}_{\infty }\) has been proven in Proposition 3 of Alama et al. (2016).

We study two distinguished limits in what follows—the \(\lambda \rightarrow 0\) limit which is relevant for polygon holes \(E_K\) with edge length comparable to \(\xi _n\) which is typically on the nanometre scale, and the \(\lambda \rightarrow \infty \) limit, which is the macroscopic limit relevant for micron-scale or larger polygonal holes. In the following sections, we study the limiting problems and the limiting minimiser profiles, their defect sets and multistability in the \(\lambda \rightarrow \infty \) limit.

3 The \(\lambda \rightarrow 0\) Limit

In Theorem 1 of Alama et al. (2016), as \(\lambda \rightarrow 0\), the solution of (2.8) converges to the unique solution of (3.1) below, with Dirichlet boundary conditions.

$$\begin{aligned} \begin{aligned}&\Delta P_{11}^0 = 0,\ \Delta P_{12}^0 = 0, on\ E_K^C,\\&P_{11}^0 = P_{11b},\ P_{12}^0 = P_{12b},\ on\ \partial E_K,\\&\textbf{P} = \textbf{P}^*,\ |x|\rightarrow \infty . \end{aligned} \end{aligned}$$
(3.1)

In other words, the limiting problem is a boundary-value problem for the Laplace equation for \(\textbf{P}\) on \(E^C_K\), in the \(\lambda \rightarrow 0\) limit. This problem is explicitly solvable and in the following sections, we exploit the symmetries of the Laplace equation, the symmetries of the polygon and boundary conditions to illustrate how the limiting profile depends on K—the number of polygon edges, and \(\gamma ^*\)—the director angle at infinity. In fact, these two parameters tune the existence, location and dimensionality of defects in this limit, amenable to experimental verification in due course.

3.1 Defect Patterns Outside a 2D Disc

As an illustrative example, we first consider the limiting problem (3.1) on the complement of a disc. As \(K\rightarrow \infty \), the domain, \(E_K^C\), converges to the exterior of a disc, \(D^C\). The conformal mapping from unit disc, D, to exterior of disc, \(D^C\), is given by

$$\begin{aligned} w = f(z) = \frac{1}{z}. \end{aligned}$$

Under the mapping \(f^{-1}:D^C\rightarrow D\), the limiting problem for \(\lambda = 0\) is given by:

$$\begin{aligned}&\Delta p_{11} = 0, \ \Delta p_{12} = 0,\ z\in D \nonumber \\&p_{11} = \frac{B}{2C}\cos (-2\theta ),\ p_{12} = \frac{B}{2C}\sin (-2\theta ),\ z\in \partial D \nonumber \\&p_{11} = \frac{B}{2C}\cos (2\gamma ^*),\ p_{12} = \frac{B}{2C}\sin (2\gamma ^*),\ z=0 \end{aligned}$$
(3.2)

where \(z = r e^{i\theta }\), \(\theta \) is the azimuthal angle and r is the radius. The corresponding solution is

$$\begin{aligned} p_{11}(re^{i\theta },\gamma ^*)&= \frac{B}{2C}\left( r^2\cos (-2\theta ) + \lim _{\epsilon \rightarrow 0}\frac{\cos (2\gamma ^*)ln r}{ln \epsilon }\right) ,\nonumber \\ p_{12}(re^{i\theta },\gamma ^*)&= \frac{B}{2C}\left( r^2\sin (-2\theta ) + \lim _{\epsilon \rightarrow 0}\frac{\sin (2\gamma ^*)ln r}{ln \epsilon }\right) . \end{aligned}$$
(3.3)

This solution has the rotational symmetry property

$$\begin{aligned} (p_{11},p_{12})(re^{i\theta },\gamma ^*)= & {} (p_{11}(re^{i\theta +i\gamma ^*},0)\cos 2\gamma ^* - p_{12}(re^{i\theta +i\gamma ^*},0)\sin 2\gamma ^*,\nonumber \\{} & {} p_{12}(re^{i\theta +i\gamma ^*},0)\cos 2\gamma ^* + p_{11}(re^{i\theta +i\gamma ^*},0)\sin 2\gamma ^*),\nonumber \\ \end{aligned}$$
(3.4)

so that it suffices to assume \(\gamma ^* = 0\). With \(\gamma ^* = 0\), we can check that \(p_{11}=p_{12}=0\) at exactly two points, located at \(\theta =\frac{\pi }{2}\) and \(\theta = \frac{3\pi }{2}\), respectively. The corresponding limiting solution on \(D^C\) is:

$$\begin{aligned} \textbf{P}(w{, \gamma ^*})= & {} \textbf{p}(f^{-1}(w){, \gamma ^*})\nonumber \\= & {} s_+\left( \frac{1}{\rho ^2}(\textbf{e}_{\rho }\otimes \textbf{e}_{\rho }-\textbf{I}_2/2 )+\lim _{\epsilon \rightarrow 0}\frac{ln \rho }{ln 1/\epsilon }(\textbf{n}^*\otimes \textbf{n}^* -\textbf{I}_2/2)\right) \end{aligned}$$
(3.5)

where \(w = \rho e^{i\psi }\),\(\textbf{e}_{\rho }= (\cos \psi ,\sin \psi )\) is the unit radial vector and the constraint at infinity is \(\textbf{n}^* = \left( \cos \gamma ^*,\sin \gamma ^*\right) \).

This method can be easily generalised to piecewise constant boundary conditions on segments of \(\partial D\), relevant for solving the limiting problem (3.1) on \(E^C_K\), as we shall see in the Section 3.3. Consider the following boundary-value problem on the unit disc D,

$$\begin{aligned} \Delta u&= 0,\ z\in D, \nonumber \\&u = u_k,\ on\ z\in D_k,\ k = 1,\cdots , K,\nonumber \\&u = u_0,\ on\ z = 0. \end{aligned}$$
(3.6)

where

$$\begin{aligned} D_k = \{e^{i\theta },\theta \in (-2\pi k/K, -2\pi k/K + 2\pi /K)\},\ k = 1,\ldots ,K, \end{aligned}$$
(3.7)

are the segments of \(\partial D\). The solution, u, can be written as, \(u = u_a + u_b\), where \(u_a\) and \(u_b\) are defined by the following boundary-value problems:

$$\begin{aligned}&\Delta u_a = 0,\ z\in D,\nonumber \\&u_a = u_k,\ on\ z\in D_k,\ k = 1,\cdots , K, \nonumber \\&u_a = \frac{1}{2\pi }\sum _{k = 1}^K \int _{2\pi (K-k)/K}^{2\pi (K-k+1)/K}u_k d\phi \, on\ z=0, \end{aligned}$$
(3.8)

and

$$\begin{aligned}&\Delta u_b = 0,\ z\in D, \nonumber \\&u_b = 0,\ on\ z\in D_k,\ k = 1,\cdots , K, \nonumber \\&u_b = u_0-\frac{1}{2\pi }\sum _{k = 1}^K \int _{2\pi (K-k)/K}^{2\pi (K-k+1)/K}u_k d\phi \, on\ z=0. \end{aligned}$$
(3.9)

Using the Poisson integral, and standard separation of variables method for the Laplace equation on an annulus, with the radius of inner ring (denoted by \(\epsilon \)) approaching zero, the solution of (3.6) is given by

$$\begin{aligned} u(re^{i\theta })= & {} \frac{1}{2\pi }\sum _{k = 1}^K\int _{2\pi (K-k)/K}^{2\pi (K-k+1)/K}u_k\frac{1-r^2}{1-2r\cos (\phi -\theta )+r^2}d\phi \nonumber \\{} & {} + \lim _{\epsilon \rightarrow 0}\frac{(u_0- \frac{1}{2\pi }\sum _{k = 1}^K \int _{2\pi (K-k)/K}^{2\pi (K-k+1)/K}u_k d\phi )ln r}{ln \epsilon }. \end{aligned}$$
(3.10)

3.2 Schwarz–Christoffel Mapping

The example of \(D^C\) gives useful insights into the computation of the limiting profile for the exterior of a regular polygon, \(E^C_K\). The conformal mapping from D to \(D^C\) is straightforward to construct. The analogous conformal mapping from D to \(E_K^C\) is the following Schwarz–Christoffel mapping:

The mapping from D to \(E_K^C\) (exterior of a regular polygon with K sides of unit length since we have re-scaled the length variable) is defined by Driscoll and Trefethen (2002)

$$\begin{aligned} w = f(z) = A - C\int ^z x^{-2}\prod _{k = 1}^K\left( 1-\frac{x}{w_k} \right) ^{1-\alpha _k} dx,\ \forall z\in D. \end{aligned}$$
(3.11)

Here \(\alpha _k\pi \) is the exterior angle of \(E_K^C\) (interior angle of polygon \(E_K\)) and \(\omega _k\) are the polygon vertices. In particular, for the mapping from the unit disc to \(E^C_K\), \(\alpha _k = 1-2/K\) \(k = 1,\cdots ,K\), when the first vertex is located at \(w_1 = (1,0)\) and the Schwarz–Christoffel mapping is uniquely given by:

$$\begin{aligned} w = f(z) = C(K)\int ^z x^{-2}\left( 1-x^K \right) ^{2/K}dx, \end{aligned}$$
(3.12)

see Fig. 2. The leading term of f in (3.12) is \(-C(K)z^{-1}\) and hence, as \(z\rightarrow 0\), \(f(z) \rightarrow \infty \) and f is single-valued near the origin Driscoll and Trefethen (2002). The pre-factor, C(K) is real and |C(K)| is the capacity or transfinite diameter of the region \(E_K^C\)Driscoll and Trefethen (2002). One can check that f maps the circle, \(\partial D\), onto the polygon boundary, \(\partial E_K = f(\partial D)\) and the segments of \(\partial D\) to the corresponding segments of \(\partial E_K^C\), i.e.

$$\begin{aligned} f(D_k) = C_k, \end{aligned}$$
(3.13)

where \(D_k\), \(k = 1,\cdots ,K\) is defined in (3.7), and

$$\begin{aligned} f(0) = \infty \end{aligned}$$
(3.14)

The mapping, f, has the following properties

$$\begin{aligned} f(\overline{z})&= C(K)\int ^{\overline{z}} x^{-2}(1-x^K)^{2/K}dx = C(K)\int ^z \overline{x}^{-2}(1-\overline{x}^K)^{2/K}dx = \overline{f(z)}, \end{aligned}$$
(3.15)
$$\begin{aligned} f(ze^{2\pi n i/K})&= C(K)\int ^{ze^{2\pi n i/K}} x^{-2}(1-x^K)^{2/K}dx \end{aligned}$$
(3.16)
$$\begin{aligned}&= C(K)\int ^{z} e^{-4\pi n i/K}x^{-2}(1-x^K)^{2/K} e^{2\pi n i/K}dx = e^{-2\pi n i/K} f(z), \end{aligned}$$
(3.17)

which are useful when studying the properties of solutions of (3.1) on \(E^C_K\).

Fig. 2
figure 2

Schwarz–Christoffel mapping from unit disc to exterior of hexagon

3.3 The Limiting Problem for \(E^C_K\) in Terms of a Boundary-Value Problem on D

Under the Schwarz–Christoffel mapping \(f^{-1}:E_K^C\rightarrow D\), the limiting problem (3.1) can be equivalently written in terms of a rLdG tensor, \(\textbf{p}\) defined on D, as shown below:

$$\begin{aligned} \begin{aligned}&\Delta p_{11}^0 = 0,\ \Delta p_{12}^0 = 0, on\ D,\\&p_{11}^0 = p_{11b},\ p_{12}^0 = p_{12b},\ on\ D_k,\\&p_{11}(0,0) = \frac{B}{2C}\cos (2\gamma ^*),\ p_{12}(0,0) = \frac{B}{2C}\sin (2\gamma ^*),\\ \end{aligned} \end{aligned}$$
(3.18)

The Dirichlet boundary conditions (2.5) translate to boundary conditions, on \(D_k\), \(k = 1,\cdots , K\) defined in (3.7), which are segments of \(\partial D\); see below:

$$\begin{aligned} \begin{aligned}&p_{11b} = \overline{\alpha }_k = {\left\{ \begin{array}{ll} &{}\frac{B}{2C}\left( (1-\frac{\overline{a}}{2})\cos \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) \right. \\ &{}\quad \left. +\frac{\overline{a}}{2}\cos \left( \frac{\left( 2k+1\right) 2\pi }{K}\right) \right) ,\ \theta \in [\frac{2\pi (K-k)}{K},\frac{2\pi (K-k+1)}{K}-\frac{\pi }{K}],\\ &{}\frac{B}{2C}\left( (1-\frac{\overline{a}}{2})\cos \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) \right. \\ &{}\quad \left. +\frac{\overline{a}}{2}\cos \left( \frac{\left( 2k-3\right) 2\pi }{K}\right) \right) ,\ \theta \in [\frac{2\pi (K-k+1)}{K}-\frac{\pi }{K},\frac{2\pi (K-k+1)}{K}], \end{array}\right. }\\&p_{12b} = \overline{\beta }_k = {\left\{ \begin{array}{ll} &{}\frac{B}{2C}\left( (1-\frac{\overline{a}}{2})\sin \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) \right. \\ &{}\quad \left. +\frac{\overline{a}}{2}\sin \left( \frac{\left( 2k+1\right) 2\pi }{K}\right) \right) ,\ \theta \in [\frac{2\pi (K-k)}{K},\frac{2\pi (K-k+1)}{K}-\frac{\pi }{K}],\\ &{}\frac{B}{2C}\left( (1-\frac{\overline{a}}{2})\sin \left( \frac{\left( 2k-1\right) 2\pi }{K}\right) \right. \\ &{}\quad \left. +\frac{\overline{a}}{2}\sin \left( \frac{\left( 2k-3\right) 2\pi }{K}\right) \right) ,\ \theta \in [\frac{2\pi (K-k+1)}{K}-\frac{\pi }{K},\frac{2\pi (K-k+1)}{K}], \end{array}\right. } \end{aligned} \nonumber \\ \end{aligned}$$
(3.19)

where the interpolation function, \(\overline{a}(z,\sigma )\), \(\forall z\in \partial D\) satisfies \(\overline{a}(e^{i(2K-2k+1)\pi /K},\sigma ) = a((w_k+w_{k+1})/2,\sigma ) = 0\) and \(\overline{a}(e^{i2(K-k)\pi /K},\sigma ) = a(w_{k+1},\sigma ) = 1\). We assume \(\overline{a}\in C^2\) and \(\frac{\partial ^2\overline{a}}{\partial \theta ^2}\ge 0\). As \(\sigma \rightarrow 0\), \(\overline{a}\rightarrow 0\), the Dirichlet boundary conditions approach \(p_{11b} = \hat{\alpha }_k\), \(p_{12b} = \hat{\beta }_k\) on \(D_k\) uniformly, where \(\hat{\alpha }_k\) and \(\hat{\beta }_k\) are given in (2.6).

For convenience, we extend the definition of \(\overline{\alpha }_k\), \(\overline{\beta }_k\), \(k = 1,\cdots ,K\), to \(k\in \mathbb {Z}\) and use the periodicity of \(\tan \), \(\cos \) and \(\sin \) to define

$$\begin{aligned} \overline{\alpha }_{k+nK} = \overline{\alpha }_k,\overline{\beta }_{k+nK} = \overline{\beta }_k,\ n\in \mathbb {Z}. \end{aligned}$$
(3.20)

One can check the following relations between \(\overline{\alpha }_k\) and \(\overline{\alpha }_{n+k}\)(\(\overline{\alpha }_{K-k+1}\)), \(\overline{\beta }_k\) and \(\overline{\beta }_{n+k}\)(\(\overline{\beta }_{K-k+1}\)):

$$\begin{aligned} \overline{\alpha }_{n+k}&= \cos (4\pi n/K)\overline{\alpha }_k-\sin (4\pi n/K)\overline{\beta }_k, \end{aligned}$$
(3.21)
$$\begin{aligned} \overline{\beta }_{n+k}&= \sin (4\pi n/K)\overline{\alpha }_k+\cos (4\pi n/K)\overline{\beta }_k, \end{aligned}$$
(3.22)
$$\begin{aligned} \overline{\alpha }_{K-k+1}&= \overline{\alpha }_k, \end{aligned}$$
(3.23)
$$\begin{aligned} \overline{\beta }_{K-k+1}&= -\overline{\beta }_k. \end{aligned}$$
(3.24)

The constants \(\hat{\alpha }_k\) and \(\hat{\beta }_k\) have similar properties. From (3.24), we have

$$\begin{aligned} \int _{0}^{2\pi } p_{12b} d\theta&= \sum _{k = 1}^K \overline{\beta }_k = \frac{1}{2}\sum _{k = 1}^K\overline{\beta }_k+\frac{1}{2}\sum _{k = 1}^{K}\overline{\beta }_{K-k+1} = 0. \end{aligned}$$
(3.25)

From (3.22) and (3.25), we have

$$\begin{aligned} \int _{0}^{2\pi } p_{11b} d\theta&= \sum _{k = 1}^K \overline{\alpha }_k = \sum _{k = 1}^K \frac{\overline{\beta }_{k+1}-cos(4\pi /K)\overline{\beta }_k}{\sin (4\pi /K)} \nonumber \\&= \frac{1}{\sin (4\pi /K)}\sum _{k = 1}^K \overline{\beta }_{k+1}-\frac{1}{tan(4\pi /K)}\sum _{k = 1}^K \overline{\beta }_{k} = 0,\ for\ K\ne 4. \end{aligned}$$
(3.26)

Additionally, for \(K =4\), \(\overline{\alpha }_k = 0\), \(k = 1,\cdots ,4\), i.e. \(\int _{0}^{2\pi } p_{11b} d\theta = 0\).

Hence, following the solution of (3.6) with Dirichlet boundary conditions on \(\partial D\), and constraint at origin in (3.10), and the results in (3.25)–(3.26), the solution of (3.18) is given by

$$\begin{aligned} p_{11}(r,\theta ,\gamma ^*)&= \frac{1}{2\pi }\sum _{k = 1}^K\int _{2\pi (K-k)/K}^{2\pi (K-k+1)/K}\overline{\alpha }_{k}\frac{1-r^2}{1-2r\cos (\phi -\theta )+r^2}d\phi \nonumber \\&\quad + \lim _{\epsilon \rightarrow 0}\frac{B}{2C}\frac{\cos (2\gamma ^*)ln r}{ln \epsilon }\nonumber \\&= \frac{1}{2\pi }\sum _{k = K}^1\int _{2\pi (k-1)/K}^{2\pi k/K}\overline{\alpha }_{K-k+1}\frac{1-r^2}{1-2r\cos (\phi -\theta )+r^2}d\phi \nonumber \\&\quad + \lim _{\epsilon \rightarrow 0}\frac{B}{2C}\frac{\cos (2\gamma ^*)ln r}{ln \epsilon }\nonumber \\&= \frac{1}{2\pi }\sum _{k = 1}^K\int _{2\pi (k-1)/K}^{2\pi k/K}\overline{\alpha }_k\frac{1-r^2}{1-2r\cos (\phi -\theta )+r^2}d\phi \nonumber \\&\quad + \lim _{\epsilon \rightarrow 0}\frac{B}{2C}\frac{\cos (2\gamma ^*)ln r}{ln \epsilon } \end{aligned}$$
(3.27)
$$\begin{aligned} p_{12}(r,\theta ,\gamma ^*)&= \frac{1}{2\pi }\sum _{k = 1}^K\int _{2\pi (K-k)/K}^{2\pi (K-k+1)/K}\overline{\beta }_{k}\frac{1-r^2}{1-2r\cos (\phi -\theta )+r^2}d\phi \nonumber \\&\quad + \lim _{\epsilon \rightarrow 0}\frac{B}{2C}\frac{\sin (2\gamma ^*)ln r}{ln \epsilon }\nonumber \\&= \frac{1}{2\pi }\sum _{k = 1}^K\int _{2\pi (k-1)/K}^{2\pi k/K}\overline{\beta }_{K-k+1}\frac{1-r^2}{1-2r\cos (\phi -\theta )+r^2}d\phi \nonumber \\&\quad + \lim _{\epsilon \rightarrow 0}\frac{B}{2C}\frac{\sin (2\gamma ^*)ln r}{ln \epsilon }\nonumber \\&= -\frac{1}{2\pi }\sum _{k = 1}^K\int _{2\pi (k-1)/K}^{2\pi k/K}\overline{\beta }_k\frac{1-r^2}{1-2r\cos (\phi -\theta )+r^2}d\phi \nonumber \\&\quad + \lim _{\epsilon \rightarrow 0}\frac{B}{2C}\frac{\sin (2\gamma ^*)ln r}{ln \epsilon }. \end{aligned}$$
(3.28)

In the above, we use (3.23) and (3.24), along with standard changes of variables, i.e. \(s - 1 = K - k\), and swap the summation variable from s to k etc. The following proposition is crucial for restricting \(\gamma ^*\) to a specified range, for a given K, using rotation and reflection symmetries, the proof of which deferred to an appendix, so as not to detract from the main details. The following corollaries describe the symmetry properties of the limiting solution \(\textbf{P}\) on \(E_K^C\), for any K.

Proposition 3.1

We can restrict \(\gamma ^*\in [0,\frac{\pi }{K}]\), since there are rotation relations between

\((p_{11}, p_{12})\vert _{(re^{i\theta +2\pi ki/K},\gamma ^*-\frac{2\pi k}{K})}\), \(k = 1,\cdots ,K\), and \((p_{11}, p_{12})\vert _{(re^{i\theta },\gamma ^*)}\), and reflection relations between \((p_{11}, p_{12})\vert _{(re^{i\theta },\gamma ^*)}\) and \((p_{11}, p_{12})\vert _{(re^{-i\theta },-\gamma ^*)}\).

Proof

See Appendix. \(\square \)

Corollary 3.1

If K is even, using (6.1) and (6.3) with \(n = K/2\), we have the symmetry property \(\textbf{p}(re^{i\theta +i\pi },\gamma ^*) = \textbf{p}(re^{i\theta },\gamma ^*)\). Due to the property of the SC mapping f in (3.17), \(f(re^{i\theta +i\pi }) = e^{-i\pi }f(re^{i\theta }) = e^{i\pi }f(re^{i\theta })\), the corresponding \(\textbf{P}\), defined on \(E_K^C\), has the symmetry property: \(\textbf{P}(\rho e^{i\psi +i\pi },\gamma ^*) = \textbf{P}(\rho e^{i\psi },\gamma ^*)\).

Corollary 3.2

From (6.4) and (6.6) with \(\gamma ^* = 0\), we have

$$\begin{aligned} p_{11}(re^{-i\theta },0) = p_{11}(re^{i\theta },0),\qquad p_{12}(re^{-i\theta },0) = -p_{12}(re^{i\theta },0). \end{aligned}$$
(3.29)

The SC mapping in (3.12) preserves reflection symmetry, \(f(re^{-i\theta }) = \overline{f(re^{i\theta })}\). We have \((P_{11},P_{12})(\rho e^{\psi i},0) = (p_{11},p_{12})(f^{-1}(\rho e^{\psi i}),0)\), for \(w=\rho e^{\psi i}\in E_K^C\), \((P_{11},P_{12})\) has reflection symmetry about \(\psi = 0\), i.e.

$$\begin{aligned} P_{11}(\rho e^{-i\psi },0) = P_{11}(\rho e^{i\psi },0),\qquad P_{12}(\rho e^{-i\psi },0) = -P_{12}(\rho e^{i\psi },0). \end{aligned}$$
(3.30)

Corollary 3.3

With \(\gamma ^* = \pi /K\), \((P_{11},P_{12})\) has reflection symmetry about \(\psi = \pi /K\), i.e.

$$\begin{aligned} P_{11}(\rho e^{\pi /K i-\psi i},\pi /K)&= P_{11}(\rho e^{\pi /K i+ \psi i},\pi /K)\cos (4\pi /K) \nonumber \\&\quad + P_{12}(\rho e^{\pi /K i+ \psi i},\pi /K)\sin (4\pi /K), \end{aligned}$$
(3.31)
$$\begin{aligned} P_{12}(\rho e^{\pi /K i-\psi i},\pi /K)&= -P_{12}(\rho e^{\pi /K i+ \psi i},\pi /K)\cos (4\pi /K)\nonumber \\&\quad + P_{11}(\rho e^{\pi /K i + \psi i},\pi /K)\sin (4\pi /K). \end{aligned}$$
(3.32)

Proof

See Appendix. \(\square \)

In the next subsections, we apply these results to \(K=3\), \(K=4\) and generic \(E^C_K\) with \(K>4\). These specific examples demonstrate the interplay between K and \(\gamma ^*\), and how this can be exploited to tailor defect sets in reduced 2D problems.

3.3.1 The Limiting Profile for \(E^C_3\) - Exterior of an Equilateral Triangle

We first consider the solution of (3.1) on \(E^C_3\)- the complement of a regular, re-scaled equilateral triangle with homeotropic boundary conditions. Recall the SC mapping from the unit disc D to \(E^C_3\) given by \(f(z) = f(re^{i\theta }) = \rho e^{i \psi } = w\).

From Corollary 3.2, for any K, with \(\gamma ^* = 0\), we have \(p_{12}(r,0) \equiv 0\) on \(\theta =0\) or \(\pi \). The boundary condition in (3.19), \(p_{12b} = \beta _k\le 0\) with \(K=3\) on \(\partial D\cap \{\theta \in [0,\pi ]\}\). From the maximum principle for the Laplace equation on \(D\cap \{\theta \in [0,\pi ]\}\), \(p_{12}= 0\) if and only if \(\theta = 0\) or \(\pi \). Analogously, on \(D\cap \{\theta \in [\pi ,2\pi ]\}\), \(p_{12}=0\) if and only if \(\theta = \pi \) or \(2\pi \). In conclusion, defects can only appear on the diameter \(z = r\), \(r\in (-1,1)\). Using (3.26), we obtain \(p_{11}(0,\gamma ^*)-\frac{B}{2C}\cos (2\gamma ^*)c(r) = \frac{1}{2\pi }\int _0^{2\pi }p_{11b}d\theta = 0\), where \(c(r) = \lim _{\epsilon \rightarrow 0}\frac{ln r}{ln \epsilon }\). For \(K = 3\), with \(\gamma ^* = 0\), as in Lemma 4.5 in Canevari et al. (2017), we show that as r increases from 0 to 1, \(p_{11}(r,0)-\frac{B}{2C}c(r)\) monotonically decreases from 0 to \(\overline{\alpha }_1 = \frac{B}{2C}\left( (1-\overline{a}/2)\cos (2\pi /3) + \overline{a}/2\cos (-2\pi /3)\right) = \frac{B}{2C}\cos (2\pi /3) = -\frac{B}{4C}\). The term, \(\frac{B}{2C}c(r)\), monotonically decreases from \(\frac{B}{2C}\) to 0. Hence, \(p_{11}(r,0)\) monotonically decreases from B/2C to \(-B/4C\). There exists a \(r^*\) such that \((p_{11},p_{12})(r^*,0) = 0\), i.e. there is a defect on \(\theta = 0\). As r increases from 0 to 1, \(p_{11}(-r,0)-\frac{B}{C}c(r)\) is monotonically increasing from 0 to \(\overline{\alpha }_2 = \frac{B}{2C}(1-\overline{a}/2)\cos (\pi ) = \frac{B}{2C}\), since \(\overline{a} = 0\) at the middle of \(D_2\). So with \(\gamma ^*=0\), under the SC mapping \(f^{-1}\), the defect set of the limiting profile consists of a single isolated point, on \(E^C_3\).

Using the rotation-based relations in (6.1), with \(\gamma ^* = 0\), \(\theta =0\), \(n = 1\), and \((p_{11},p_{12})(r^*,0) = 0\), we obtain \((p_{11},p_{12})(r^*e^{i2\pi /3}, \pi /3)=(p_{11},p_{12})(r^*e^{i2\pi /3},-2\pi /3) = 0\), i.e. for \(\gamma ^* = \pi /3\), there is a unique defect on \(\theta = 2\pi /3\), and consequently for \(\psi = -2\pi /3\) in the limiting profile \((P_{11},P_{12})(\rho e^{i\psi },\gamma ^*)\). For \(\gamma ^* \in \left( 0, \pi /3 \right) \), the defect smoothly rotates between \(\psi =0\) and \(\psi = -2 \pi /3\) in the limiting profile (see Fig. 3) and this is a good example of how \(\gamma ^*\) can tune the location of the defect for the limiting profile. In the first row of Fig. 3, we plot the solution of (3.18), for \(K=3\). The domain is a unit disc D, since we non-dimensionalised the rLdG energy in (2.2). In the second row of Fig. 3, we plot the corresponding solution of (3.1) on \(E_3^C\), and again the polygonal hole is re-scaled to have unit length. The computational domain is the same for all values of \(\lambda \), because of the rescaling in (2.3). In Fig. 3 and subsequent figures, the colour bar labels the order parameter, \(s = \sqrt{p_{11}^2 + p_{12}^2}\), and the white lines correspond to the nematic director \(\textbf{n} = (\cos (arctan(p_{12}/p_{11})/2),\sin (arctan(p_{12}/p_{11})/2))\), where \(p_{11}\) and \(p_{12}\) are the two scalar fields in (3.18). In Fig. 4, we plot the difference between the limiting solution of (3.1) and the solution of the Euler-Lagrange equations (2.8), with \(K=3\) and \(\lambda ^2 = 0.01\). The error is proportional to \(\lambda ^2\) and we deduce that the limiting solution is a good approximation to the solutions of (2.8) for small values of \(\lambda ^2\). This can also be formally justified using further asymptotic analysis. The remaining figures in this Section have been computed with the same numerical method, with the same details. The numerical details are given in Section 4.

Fig. 3
figure 3

a The mapped solution of (3.18) on D with \(K = 3\) and b the numerical limiting solution of (3.1) on \([-10,10]^2\backslash E_3\) as an approximation to \(E_3^C\). From left to right, the boundary condition at infinity is \(\gamma ^*=0\), \(\pi /6\), and \(\pi /3\). In this figure and all subsequent figures where the value of B and C are required, we have \(B = 0.64 \times 10^4 N/m^2\) and \(C = 0.35 \times 10^4 N/m^2\). The colour bar is the scalar order parameter, \(s = \sqrt{p_{11}^2 + p_{12}^2}\), and the white lines correspond to the nematic director \(\textbf{n} = (\cos (arctan(p_{12}/p_{11})/2),\sin (arctan(p_{12}/p_{11})/2))\)

Fig. 4
figure 4

From left to right: the nematic order, s, of the limiting solution \((P_{11}^0,P_{12}^0)\) for (3.1) on \([-10,10]^2\backslash E_3\); the solution \((P_{11}^1,P_{12}^1)\) of Euler-Lagrange equations (2.8) with \(\lambda ^2 = 0.01\); comparison between the exact and limiting solutions: \(|\textbf{P}^1|^2/2-|\textbf{P}^0|^2/2 = (s^1)^2-(s^0)^2\); the difference in the nematic director \(P^1_{11}P^0_{12} - P^1_{12}P^0_{11} = s^0s^1\sin (2\gamma ^0 - 2\gamma ^1)\)

3.3.2 The Limiting Profile for \(E^C_4\) - Exterior of a Unit Square

In Kralj and Majumdar (2014), Canevari et al. (2017), the authors show that the unique limiting rLdG minimiser on \(E_4\) with tangent boundary conditions, is the WORS solution in the \(\lambda \rightarrow 0\) limit; the label “WORS” has been proposed by the authors in recognition of the fact that the WORS supports two nodal or defect lines along the square diagonals and is an example of widely studied order reconstruction solutions in LdG theory. The interested reader is referred to these papers, Kralj and Majumdar (2014), Canevari et al. (2017), for more details. We study the analogous problem on \(E_4^C\) with homeotropic boundary conditions, to study whether line defects survive on \(E_4^C\) for different choices of \(\gamma ^*\).

Proposition 3.2

The solution of (3.1) on \(E^C_4\) has two line defects for \(\gamma ^* = \pi /4 + n\pi /2\), \(n\in \mathbb {Z}\), and no interior defects otherwise.

Proof

Substituting \(K=4\) into the Dirichlet boundary condition (3.19), \(\overline{\alpha }_k\equiv 0\), \(k = 1,\cdots , 4\). According to (3.27), the solution of (3.18) on \(E^C_4\) is given by

$$\begin{aligned} p_{11}(re^{i\theta },\gamma ^*) = \frac{B}{2C}\cos (2\gamma ^*)c(r). \end{aligned}$$
(3.33)

As in Proposition 3.1, we can assume that \(\gamma \in [0,\pi /4]\). The zero set of \(p_{11}\) is \(\emptyset \) for any \(\gamma \in [0,\pi /4)\). For \(\gamma ^* = \pi /4\), we have \(p_{11}(re^{i\theta },\pi /4)\equiv 0\).

Substituting \(K = 4\), into (3.28), we have

$$\begin{aligned}{} & {} p_{12}(re^{i\theta },\gamma ^*) - \frac{B}{2C}\sin (2\gamma ^*)c(r)\nonumber \\{} & {} \quad = -\frac{1}{2\pi }\sum _{k =1}^4\int _{\pi (k-1)/2}^{\pi k/2}\overline{\beta }_k\frac{1-r^2}{1-2r\cos (\phi -\theta )+r^2}d\phi . \end{aligned}$$
(3.34)

which only depends on r and \(\theta \). In the following, we prove that with \(\gamma ^* = \pi /4\), \(p_{12}(re^{i\theta },\pi /4)-\frac{B}{2C}c(r)\) is monotonically decreasing in the r-direction. Firstly, we prove that on \(\theta = \pi \), i.e. \(y = 0\) and \(x\le 0\), \(p_{12}(re^{i\pi },\gamma ^*)-\frac{B}{2C}\sin (2\gamma ^*)c(r) \equiv 0\). Using (6.6), we obtain

$$\begin{aligned} p_{12}(re^{i \pi },\gamma ^*)-\frac{B}{2C}\sin (2\gamma ^*)c(r)= -(p_{12}(re^{-i \pi },-\gamma ^*)-\frac{B}{2C}\sin (-2\gamma ^*)c(r)).\nonumber \\ \end{aligned}$$
(3.35)

Since \(p_{12}(re^{i\theta },\gamma ^*) - \frac{B}{2C}\sin (2\gamma ^*)c(r)\) in (3.34) only depends on r and \(\theta \), we have \(p_{12}(re^{i\theta },\gamma ^*) - \frac{B}{2C}\sin (2\gamma ^*)c(r) = p_{12}(re^{i\theta },-\gamma ^*) - \frac{B}{2C}\sin (-2\gamma ^*)c(r), \) and subsequently

$$\begin{aligned} p_{12}(re^{-i \pi },-\gamma ^*)-\frac{B}{2C}\sin (-2\gamma ^*)c(r)= & {} p_{12}(re^{-i \pi },\gamma ^*)\nonumber \\ -\frac{B}{2C}\sin (2\gamma ^*)c(r)= & {} p_{12}(re^{i \pi },\gamma ^*)-\frac{B}{2C}\sin (2\gamma ^*)c(r).\nonumber \\ \end{aligned}$$
(3.36)

Combining (3.35) and (3.36), we get the desired conclusion,

$$\begin{aligned} p_{12}(re^{i \pi },\gamma ^*)-\frac{B}{2C}\sin (2\gamma ^*)c(r)=0. \end{aligned}$$
(3.37)

Similarly, we prove that on \(\theta = 3\pi /2\), i.e. \(x=0\) and \(y\le 0\), \(p_{12}(re^{i3\pi /2},\gamma ^*)-\frac{B}{2C}\sin (2\gamma ^*)c(r) \equiv 0\). Using (6.3) with \(n = 1\), \(K = 4\), \(\theta = \pi \), and \(\gamma = \pi /4\), we have

$$\begin{aligned} p_{12}(re^{i 3\pi /2},-\pi /4) = \cos (\pi )p_{12}(re^{i \pi },\pi /4) = -p_{12}(re^{i \pi },\pi /4). \end{aligned}$$

From (3.37), \(p_{12}(re^{i \pi },\pi /4) = \frac{B}{2C}\sin (\pi /2)c(r)\) and hence,

$$\begin{aligned} p_{12}(re^{i3\pi /2},\gamma ^*)-\frac{B}{2C}\sin (2\gamma ^*)c(r)&=p_{12}(re^{i3\pi /2},-\pi /4)-\frac{B}{2C}\sin (-\pi /2)c(r)\nonumber \\&= -p_{12}(re^{i \pi },\pi /4) + \frac{B}{2C}\sin (\pi /2)c(r) \equiv 0. \end{aligned}$$
(3.38)
Fig. 5
figure 5

The mapped domain D, with \(K = 4\). We show the sign of \(\overline{\beta }_k\), \(k = 1,\cdots , 4\), the boundary condition of \(p_{12}\) on \(D_k\), the segment of \(\partial D\), and \(\overline{\beta }_2, \overline{\beta }_4\le 0\), \(\overline{\beta }_1, \overline{\beta }_3\ge 0\). The quarter of disc \(D_q = D\cap \{x\le 0\}\cap \{y\le 0\}\) is shown in dark blue (Color figure online)

On the quadrant \(\theta \in [\pi ,3\pi /2]\) or the quarter disc, \(D_q = D\cap \{x\le 0\}\cap \{y\le 0\}\) (see Fig. 5), the function, \(p_{12}(r e^{i\theta },\pi /4)-\frac{B}{2C}\sin (\pi /2)c(r)\), is the solution of

$$\begin{aligned}&\Delta u_a = 0,\ on\ D_q,\nonumber \\&u_a = \overline{\beta }_2,\ on\ D_2,\nonumber \\&u_a = 0,\ on\ x = 0\ and\ y = 0. \end{aligned}$$
(3.39)

where \(D_2\) is the segment of \(\partial D\) with \(x\le 0\) and \(y\le 0\), \(\overline{\beta }_2\le 0\) is the corresponding Dirichlet boundary condition of \(p_{12}\) on \(D_2\). Following the arguments from Lemma 4.5 in Canevari et al. (2017), the first step is to show that \(u_a\) is non-increasing in the r direction, i.e.

$$\begin{aligned} u_a(r,\theta )\ge u_a(\tau r, \theta ),\ \forall \tau >1,\ (r,\theta )\in D_q. \end{aligned}$$
(3.40)

We consider Problem (3.41), the analogue of Problem (3.39) on the extended domain, \(D_q^{\tau } = (\tau r,\theta )\) where \((r,\theta )\in D_q\),

$$\begin{aligned}&\Delta u_a^{\tau } = 0,\ on\ D_q^{\tau },\nonumber \\&u_a^{\tau } = \overline{\beta }_2,\ on\ D_2^{\tau },\nonumber \\&u_a^{\tau } = 0,\ on\ x = 0\ and\ y = 0. \end{aligned}$$
(3.41)

where \(D_2^{\tau } = (\tau r,\theta )\), with \((r,\theta )\in D_2\). By scaling invariance, the unique non-negative solution to Problem (3.41) is given by

$$\begin{aligned} u_a^{\tau }(\tau r,\theta ):=u_a(r,\theta )\ for\ any\ (r,\theta )\in D_q. \end{aligned}$$
(3.42)

Moreover, the function

$$\begin{aligned} u_a^s(r,\theta ) = {\left\{ \begin{array}{ll} &{}u_a\ on\ D_q,\\ &{}\overline{\beta }_2(\theta ),\ on\ \{r>1\}\cup \{x\le 0\}\cup \{y\le 0\}\\ \end{array}\right. } \end{aligned}$$
(3.43)

is a subsolution of (3.41), with \(\frac{\partial ^2 \overline{a}}{\partial \theta ^2}\ge 0\), so that

$$\begin{aligned} u_a^{\tau }(\tau r,\theta )\ge u_a^s(\tau r,\theta )\ for\ any\ (r,\theta )\in D_q\ s.t.\ (\tau r,\theta )\in D_q \end{aligned}$$
(3.44)

Recalling (3.42) and using \(u_a^s = u_a\) on \(D_q\), we conclude the proof of (3.40). Let \(v_a = \partial u_a/\partial r\). By (3.40), \(v_a\le 0\) on \(D_q\); we want to prove that the strict inequality holds. We differentiate Equation (3.39) with respect to r, so that \(\nabla v_a = 0\ on\ D_q.\) By the strong maximum principle, we deduce that either \(v_a\equiv 0\) or \(v_a<0\) in \(D_q\). The first possibility is clearly inconsistent with the boundary conditions for \(u_a\) in (3.39), and hence, \(v_a\) must be strictly negative inside \(D_q\).

Therefore, on the quadrant \(\theta \in [\pi ,3\pi /2]\), \(p_{12}(re^{i\theta },\pi /4)-\frac{B}{2C}\sin (\pi /2)c(r)\) is monotonically decreasing in r directions from 0 to \(\overline{\beta }_2\). Using analogous arguments, \(p_{12}(re^{i\theta },\pi /4)-\frac{B}{2C}\sin (\pi /2)c(r)\) is monotonically decreasing in the r-direction, from 0 to \(\overline{\beta }_4\) for on \(\theta \in [0,\pi /2]\), and monotonically increasing in r, from 0 to \(\overline{\beta }_1\) or \(\overline{\beta }_3\), on the quadrants \(\theta \in [\pi /2,\pi ]\cup [3\pi /2,2\pi ]\).

Moreover, \(\frac{B}{2C}\sin (\pi /2)c(r)\) is monotonically decreasing from \(\frac{B}{2C}\) to 0, hence there exists \(r^*(\theta )\) for each \(\theta \in [0,\pi /2]\cup [\pi ,3\pi /2]\) s.t. \(p_{12}(r^*(\theta )e^{i\theta },\pi /4) = 0\), so that there are two line defects in limiting profile for \(E^C_K\) with \(\gamma ^* = \pi /4\). Using (6.1) and (6.3), \((p_{11},p_{12})(re^{i\theta + i 2\pi n/K},\gamma ^*-\frac{2\pi n}{K})\), \(n\in \mathbb {Z}\), can be obtained by rotating \((p_{11},p_{12})(re^{i\theta },\gamma ^*)\). Hence, there are two line defects for the limiting profile \((P_{11},P_{12})(\rho e^{i\psi },\gamma ^*)\) on \(E^C_4\), for \(\gamma ^* = \pi /4 + n\pi /2\), \(n\in \mathbb {Z}\). \(\square \)

We plot the mapped solution of (3.18) on D with \(K = 4\) and corresponding limiting solution of (3.1) on \(E_4^C\) with \(\gamma ^* = 0\) and \(\pi /4\) in Fig. 6, respectively. On \(E^C_4\), with \(\gamma ^* = \pi /4\), there are two line defects; for \(\gamma ^* = 0\), there are no interior defects on \(E^C_4\) although there are defects at the square vertices where the nematic director cannot be uniquely defined.

Fig. 6
figure 6

The mapped solution of (3.18) on D with \(K = 4\) and the corresponding limiting solution for (3.1) on \(E_4^C\) with (a) \(\gamma ^* = 0\) and (b) \(\gamma ^* = \pi /4\)

3.3.3 The Limiting Profile for \(E^C_K\) with \(K>4\)

Lemma 3.3

For \(\gamma ^*\in [0,\pi /K]\), \(K\in \mathbb {Z}\) and \(K>4\), the solution of (3.18) has no defect on \(D_{q1} = D\cup \{2\pi -\frac{\pi }{K}\le \theta \le 2\pi \}\) and \(D_{q2} = D\cup \{\pi -\frac{\pi }{K}\le \theta \le \pi \}\).

Proof

As in Proposition 3.2, we can prove that on \(D_{q1}\) and \(D_{q2}\), \(p_{12}(re^{i\theta },\gamma ^*)-\frac{B}{2C}\sin (2\gamma ^*)c(r)\) increases monotonically in the r-direction from 0 to \(\beta _1>0\), as r increases from 0 to 1. Hence, \(p_{12}(re^{i\theta },\gamma ^*)-\frac{B}{2C}\sin (2\gamma ^*)c(r)\) is positive for \(r>0\) on \(D_{q1}\) and \(D_{q2}\). For \(K>4\) and \(K\in \mathbb {Z}\), \(\frac{B}{2C}\sin (2\gamma ^*)c(r)\) is non-negative, and is zero if and only if \(r=1\), i.e. on \(\partial D\). Hence, \(p_{12}\) is nonzero on \(D_{q1}\) or \(D_{q2}\), i.e. there is no defect on \(D_{q1}\) or \(D_{q2}\). \(\square \)

As \(p_{12}\) and \(p_{11}\) are analytic on D, the angle \(\gamma = arctan(p_{12}/p_{11})/2\) is continuous. We can separate the disc into four parts by two smooth curves (see Fig. 7): on \(curve_1\) (from \(z=e^{i(K-1)\pi /K}\) to \(z= 0\) to \(z=1\)), \(\gamma \) decreases from \(\pi /K\) to \(\gamma ^*\) to 0 monotonically, and on \(curve_2\) (from \(z=-1\) to \(z= 0\) to \(z=e^{-i\pi /K}\)), \(\gamma \) increases from 0 to \(\gamma ^*\), and then to \(\pi /K\) monotonically. Since there is no defect on \(D_{q1}\) and \(D_{q2}\) as proved in Lemma 3.3, there is no defect on \(curve_1\) and \(curve_2\). We define \(D_u\) to be the part of D above \(curve_1\), and \(D_b\) to be the part of D below \(curve_2\). The boundary \(\partial D_u\) (\(\partial D_b\)) consists of a segment of \(\partial D\) and the interior curve \(curve_1\) (\(curve_2\)). In the following, we prove that the director angle \(\gamma \) is strictly monotonic on \(\partial D\) and subsequently, that \(\gamma \) is strictly monotonic on either \(\partial D_u\) or \(\partial D_b\), and rotates by \(\pi \).

Fig. 7
figure 7

The partition of domain D. \(D_u\) and \(D_b\) are surrounded by blue and pink curves, respectively. \(D_{q1} = D\cap \{2\pi -\frac{\pi }{K}\le \theta \le 2\pi \}\) and \(D_{q2} = D\cap \{\pi -\frac{\pi }{K}\le \theta \le \pi \}\) are domains between the two dashed black lines and solid black curves in D. The black two-head arrow in the centre of D represents the direction \(\textbf{n}\) of constraint \((p_{11},p_{12})(0,\gamma ^*) = (\cos 2\gamma ^*,\sin 2\gamma ^*)\). The arcs with different colours in the outer circle represents the segment of boundary \(D_k\), \(k = 1,\cdots ,K\). The two-head arrows on the boundary \(D_k\) represent the direction \(\textbf{n}\) of boundary conditions; here, we take hexagon as an example (Color figure online)

Lemma 3.4

On \(\partial D\), the angle \(\gamma = arctan(p_{12}/p_{11})/2\) is strictly monotonic, and there is no defect.

Proof

On \(D_{1}\cap \theta \in [2\pi -\pi /K,2\pi ]\), the boundary condition is

$$\begin{aligned}&(p_{11b},p_{12b}) = (\overline{\alpha }_1,\overline{\beta }_1)\\&\quad = \frac{B}{2C}\left( (1-\frac{\overline{a}}{2})\cos \left( \frac{2\pi }{K}\right) + \frac{\overline{a}}{2}\cos \left( \frac{-2\pi }{K}\right) \right. ,\left. (1-\frac{\overline{a}}{2})\sin \left( \frac{2\pi }{K}\right) + \frac{\overline{a}}{2}\sin \left( \frac{-2\pi }{K}\right) \right) \\&\quad = \frac{B}{2C}\left( \cos \left( \frac{2\pi }{K}\right) ,(1-\overline{a})\sin \left( \frac{2\pi }{K}\right) \right) . \end{aligned}$$

The angle, \(\gamma = arctan((1-\overline{a})tan(\frac{2\pi }{K}))/2\) monotonically increases from 0 to \(\frac{\pi }{K}\) as \(\overline{a}\) decreases from 1 to 0, i.e. \(\theta \) decreases from \(2\pi \) to \(2\pi -\frac{\pi }{K}\). The angle \(\gamma \) on \(\partial D\) rotates by \(-2\pi k/K\), as \(\theta \) rotates by \(2\pi k/K\) and reflects about axis \(\theta = n\pi /K\), \(n\in \mathbb {Z}\). Hence, \(\gamma \) on \(D_b\cap \partial D\) increases monotonically from \(\pi /K\) to \(\pi \), in a clockwise fashion. Since \(\gamma \) on \(curve_2\) increases monotonically from 0 to \(\pi /K\), then \(\gamma \) on \(\partial D_b\) rotates by \(\pi \) (analogous remarks apply to \(D_u\)).

Consider \(p_{11b}^2 + p_{12b}^2\) on \(D_1\), \(\theta \in [2\pi -\pi /K,2\pi ]\), for \(K>4\), i.e.

$$\begin{aligned} p_{11b}^2 + p_{12b}^2&= \frac{B^2}{4C^2}\left( (1-\frac{\overline{a}}{2})\cos \left( \frac{2\pi }{K}\right) + \frac{\overline{a}}{2}\cos \left( \frac{-2\pi }{K}\right) \right) ^2 \\&\quad +\frac{B^2}{4C^2}\left( (1-\frac{\overline{a}}{2})\sin \left( \frac{2\pi }{K}\right) + \frac{\overline{a}}{2}\sin \left( \frac{-2\pi }{K}\right) \right) ^2\\&\quad = \frac{B^2}{4C^2}\left( (\overline{a}^2-2\overline{a})\sin ^2\frac{2\pi }{K} + 1 \right) >0. \end{aligned}$$

From the rotation-based relations in (3.21) and (3.22), \(p_{11b}^2 + p_{12b}^2\) is nonzero on \(D_k\), \(\theta \in [2\pi (K-k+1)/K - \pi /K,2\pi (K-k+1)/K]\), \(k = 1,\cdots ,K\). From the reflection relations in (3.23) and (3.24), \(p_{11b}^2 + p_{12b}^2\) is nonzero on \(D_k\), \(k = 1,\cdots ,K\). Hence, there are no defects (or zeros of \(\textbf{p}\)) on \(\partial D\), for any regular polygon with \(K>4\) edges. \(\square \)

Theorem 3.1

If \(\textbf{p}\) is a minimiser of \(J = \int _{D_u}|\nabla \textbf{p}|^2 \text {dA}\) in \(\mathcal {A} = \{\textbf{p}\in W^{1,2}(D_u;\mathbb {R}^2):\textbf{p} = \textbf{p}_{bu}\ on\ \partial D_u\}\), then there is a unique point defect in \(D_u\). Similarly, there is a unique point defect in \(D_b\). The solution of (3.18) has two zeroes or point defects and hence, under the SC mapping f (3.12), the limiting profile has two point defects in exterior of regular polygon with \(K>4\) edges.

Proof

Let Y(l) for \(0\le l \le L\) be a one-to-one parametrisation of \(\partial \textrm{D}_u\) with respect to arclength. Considering the Dirichlet data, \((p_{11bu},p_{12bu})\), in \(C^{2,\mu }(\partial \textrm{D}_u;\mathbb {R}^2)\), we have

$$\begin{aligned} (p_{11bu},p_{12bu})(Y(l)) = (s(l) cos2\gamma (l), s(l) sin2\gamma (l)). \end{aligned}$$
(3.45)

As proved in Lemma 3.3 and 3.4 , we have

$$\begin{aligned} s(l) > 0,\ \gamma '(l)\ne 0\ for\ l\in [0,L],\ and\ |2\gamma (L)- 2\gamma (0)| = 2\pi . \end{aligned}$$
(3.46)

For \(\alpha \in \mathbb {R}\) and \(\textbf{p}\), a minimiser for J in \(\mathcal {A}\), set

$$\begin{aligned} w_{\alpha }(X) = -p_{11}(X)\sin (\alpha ) + p_{12}(X)\cos (\alpha ). \end{aligned}$$
(3.47)

Define the nodal set of \(w_{\alpha }\) by

$$\begin{aligned} N_{\alpha }\equiv \{X\in \overline{D}_u:w_{\alpha }(X) = 0\}. \end{aligned}$$
(3.48)

Note that for any pair \(\alpha _1\), \(\alpha _2\) with \(0 \le \alpha _1< \alpha _2 < \pi \), the set of zeros of \(\textbf{p}\) is just \(\Gamma (\textbf{p}) = N_{\alpha _1}\cap N_{\alpha _2}\). Thus as \(\alpha \) varies, \(\Gamma (\textbf{p})\) is the subset of \(N_{\alpha }\) that remains fixed. Analogous to Lemma 2.2 in Bauman et al. (1993), we can prove for each \(\alpha \in [0,2\pi ]\), \(N_{\alpha }\) is a \(C^1\) imbedded curve in \(\overline{D}_u\), which enters and exits \(\overline{D}_u\), at distinct points of \(\partial D_u\). The following proof is identical to Theorem 2.3 in Bauman et al. (1993). Let \(\{P_1, P_2\} = N_{0}\cap \partial \textrm{D}_u\) and assume without loss of generality that \(\sin (2\gamma )\vert _{P_1} = 0\) with \(\cos (2\gamma )\vert _{P_1}<0\) and \(\sin (2\gamma )\vert _{P_2} = 0\) with \(\cos (2\gamma )\vert _{P_2}>0\). The points \(P_1\) and \(P_2\) partition \(\partial D_u\) into two arcs, \((\partial D_u)^-\equiv \{X\in \partial D:\sin (2\gamma )\vert _X\le 0\}\) and \((\partial D_u)^+\equiv \{X\in \partial D:\sin (2\gamma )\vert _X\ge 0\}\). Starting at \(P_1\) and moving along \(N_0\), let \(X_0\) be the first point reached in \(\Gamma (\textbf{p})\). Since \(X_0\in N_{\alpha }\) for all \(\alpha \), each \(N_{\alpha }\) can be parametrised by arclength, \(X = X(\tau , \alpha )\), such that

$$\begin{aligned}{} & {} a(\alpha )< \tau< b(\alpha ),\ a(\alpha ) < 0, X(a(\alpha ),\alpha )\in (\partial \textrm{D}_u)^-,\ b(\alpha )>0,\\{} & {} \quad X(b(\alpha ),\alpha )\in (\partial \textrm{D}_u)^+, \left| \frac{\partial X}{\partial \tau }\right| = 1,\ and\ X(0,\alpha ) = X_0. \end{aligned}$$

Set \(\mathcal {D} = \{(\tau ,\alpha ):a(\alpha )\le \tau \le b(\alpha ),0\le \alpha \le \pi \}\). In part 1 of Theorem 2.3, the authors proved \(X\in C^1(\mathcal {D})\), and \(a(\alpha )\) and \(b(\alpha )\) are in \(C^1([0,\pi ])\). In part 2, they have that for all \(\alpha \) in \([0,\pi ]\) and \(\tau <0\), \(X(\tau ,\alpha )\notin \Gamma (\textbf{p})\). In particular, \(X(\tau ,\pi )\notin \Gamma (\textbf{p})\). In part 3, now \(X(\tau , 0)\) and \(X(\tau , \pi )\) are each parametrisations of \(N_0\) in opposite directions. Since \(X_0\) is the first point in each direction along \(N_0\) that is in \(\Gamma (\textbf{p})\), we conclude that \(\Gamma (\textbf{p}) = {X_0}\), i.e. there is a unique point defect in \(D_u\). Similarly, there is a unique point defect in \(D_b\). Under the SC mapping f (3.12), there are two point defects in exterior of regular polygon with \(K>4\) edges. \(\square \)

Remark

Theorem 3.1 does not work for \(E^C_K\) with \(K = 3\) or 4. For \(K = 3\), since there exists a \(r^*\) s.t. \((p_{11},p_{12})(r^*e^{i2\pi /3}, \pi /3) = (p_{11},p_{12})(r^*, 0) = 0\), we may have a defect on \(D_{q1}\) or \(D_{q2}\) which does not satisfy Lemma 3.3. For \(K = 4\), with the given boundary conditions, we have defects at the four vertices, \((p_{11},p_{12}) = \frac{B}{2C}\left( cos(\frac{2\pi }{K}),(1-\overline{a})\sin (\frac{2\pi }{K})\right) = (0,1-\overline{a}(0)) = (0,0)\), which does not satisfy Lemma 3.3 and 3.4 .

Fig. 8
figure 8

The mapped solution of (3.18) on D with \(K = 5,6,\infty \) and the corresponding limiting solution of (3.1) on \(E_K^C\), the complement or exterior domain of a pentagon, hexagon and disc with (a) \(\gamma ^* = 0\) and (b) \(\gamma ^* = \pi /K\)

The mapped solution of (3.18) with \(K = 3,4,5,6,\infty \), and the corresponding solutions of (3.1) in the exterior of a triangle, square, pentagon, hexagon and disc, with \(\gamma ^* = 0\) and \(\gamma ^*=\pi /K\), are plotted in Fig. 3, Fig. 6 and Fig. 8, respectively. With even K, the solutions have the symmetry property in Corollary 3.1. With \(\gamma ^* = 0\) and \(\pi /K\), the solutions have the symmetry properties in Corollary 3.2 and Corollary 3.3, respectively. The mapped solution has two \(-1/2\) defects on two sides of \(\theta = \gamma ^*+\pi /2\) as in Theorem 3.1, for \(K>4\).

4 The \(\lambda \rightarrow \infty \) Limit

In this section, we use heuristics and numerical methods to study stable rLdG equilibria in the \(\lambda \rightarrow \infty \) limit. This limit is relevant for polygonal holes whose edge length, \(\bar{\lambda }\), is much greater than the nematic correlation length, i.e. micron-scale holes. Consider the rLdG energy in (2.3); as \(\lambda \rightarrow \infty \), minimisers of the rLdG energy will belong to the set of minimisers of the bulk energy, away from defects. Recall that the bulk energy density is given by

$$\begin{aligned} f_b(\textbf{P}) = -\frac{B^2}{4C} \text {tr}\textbf{P}^2 + \frac{C}{4}\left( \text {tr}\textbf{P}^2 \right) ^2. \end{aligned}$$

More precisely, in Theorem 2 of Alama et al. (2016), the authors prove that global minimisers of the rLdG energy in (2.3), converge strongly in \(W^{1,2}\), to

$$\begin{aligned} \textbf{P}^{\infty } = \frac{B}{C}\left( \textbf{n}^{\infty }\otimes \textbf{n}^{\infty }-\frac{1}{2}\textbf{I}_2\right) , \end{aligned}$$

where

$$\begin{aligned} \textbf{n}^{\infty } = \left( \cos \gamma ^*,\sin \gamma ^*\right) , \end{aligned}$$
(4.1)

and \(\gamma ^*\) is the solution of the following equations

$$\begin{aligned} \begin{aligned}&\Delta \gamma = 0,\ on\ E_K^C\\&\gamma = \gamma _b,\ on\ \partial E_K,\\&\gamma = \gamma ^*,\ |x|\rightarrow \infty , \end{aligned} \end{aligned}$$
(4.2)

everywhere away from the vertices of \(E_K\). The convergence is actually stronger but that is not the focus of this paper.

Fig. 9
figure 9

vertex classification according to director profile: a splay and b bend

We use this straightforward characterisation to compute an (under)estimate of the number of stable rLdG equilibria on \(E_C^K\), in the \(\lambda \rightarrow \infty \) limit. Our discussion does not cover local energy minimisers but does serve to illustrate the rich prospects of multistability tailored by two model parameters - K and \(\gamma ^*\). Multistability arises from the fact that for a given \(\gamma ^*\), there are multiple choices of \(\gamma _b\) consistent with the homeotropic boundary conditions in (2.6) as \(\sigma \rightarrow 0\), and each choice of \(\gamma _b\) yields a limiting rLdG minimiser profile in the \(\lambda \rightarrow \infty \) limit. We estimate the number of limiting profiles based on two assumptions - (i) we restrict \(\gamma _b\) so that \(\gamma \) rotates by either \(2\pi /K\) or \(2\pi /K-\pi \) at a vertex (see Fig. 9; referred to as “splay” and “bend” vertices, respectively). This rotation is required to match the homeotropic boundary conditions at the two intersecting edges. (ii) We assume that \(\text {deg}\left( \textbf{n}_b,\partial E_K^C\right) = 0\), where \(\textbf{n}_b = (\cos \gamma _b,\sin \gamma _b)\), so that if \(n_b\) denotes the number of “bend” vertices, then we necessarily have

$$\begin{aligned} \text {deg}\left( \textbf{n}_b,\partial E_K^C\right) = (K-n_b)\left( 2\pi /K\right) +n_b\left( 2\pi /K-\pi \right) = 0, \end{aligned}$$

equivalent to \(n_b = 2\). The two assumptions ensure that \(\textbf{n}_b\) has the minimal rotation around \(\partial E^C_K\) consistent with the homeotropic boundary conditions. We can of course, have scenarios wherein \(\gamma _b\) rotates by a greater amount or when the degree of \(\textbf{n}_b\) is nonzero, where \(\textbf{n}_b\) is as above, but it is reasonable to conjecture that rLdG energy minimisers have the simplest admissible topology. Consider \(E_K^C\); based on the assumptions above, we have 2 bend vertices on \(E_K^C\) for each \(K \ge 3\) and hence, at least \(\left( {\begin{array}{c}K\\ 2\end{array}}\right) \) limiting profiles and [K/2] classes of solutions which are not related by rotation or reflection in this limit. This is identical to the combinatorial estimate for the number of stable equilibria for the interior problem on \(E_K\) in Han et al. (2020), in the \(\lambda \rightarrow \infty \) limit. However, this estimate does not account for the constraint at infinity (2.7), and multistability is enhanced for the exterior problem compared to the interior problem, by virtue of \(\gamma ^*\).

For \(K=4\), there are 6 different choices for the two “bend" vertices, (i) 2 of which correspond to the two pairs of diagonally opposite vertices, (ii) 4 of which correspond to pairs of adjacent vertices. We refer to (i) as diagonal, D states, (ii) as rotated R states, following the nomenclature in Luo et al. (2012), Tsakonas et al. (2007).

The diagonal states are defined by the following choices of \(\gamma _b\) on \(E_4\):

$$\begin{aligned} \gamma _{b1} = {\left\{ \begin{array}{ll} &{}\pi /4,\ on\ C_1,\\ &{} -\pi /4,\ on\ C_2,\\ &{} \pi /4,\ on\ C_3,\\ &{} -\pi /4\ on\ C_4, \end{array}\right. } \gamma _{b2} = {\left\{ \begin{array}{ll} &{}\pi /4,\ on\ C_1,\\ &{} 3\pi /4,\ on\ C_2,\\ &{} \pi /4,\ on\ C_3,\\ &{} 3\pi /4\ on\ C_4. \end{array}\right. } \end{aligned}$$
(4.3)
Fig. 10
figure 10

The limiting profiles for \(E^C_4\), in the \(\lambda \rightarrow \infty \) limit, with \(\gamma ^*= n\pi \) (a) and \(\gamma ^* = n\pi + \pi /4\) (b). The diagonal states are plotted in the first row of a, b. The rotated states are plotted in the second row of a, b. The white lines represent the nematic director \(\textbf{n} = (\cos (\gamma ),\sin (\gamma ))\), and the red colour indicate the nematic order \(s \equiv \frac{B}{2C}\) (Color figure online)

There are four different choices of \(\gamma _b\), corresponding to the four R rotated states, as enumerated below:

$$\begin{aligned} \gamma _b= & {} {\left\{ \begin{array}{ll} &{}\pi /4,\ on\ C_1,\\ &{} -\pi /4,\ on\ C_2,\\ &{} \pi /4,\ on\ C_3,\\ &{} 3\pi /4\ on\ C_4, \end{array}\right. } \gamma _b = {\left\{ \begin{array}{ll} &{}\pi /4,\ on\ C_1,\\ &{} 3\pi /4,\ on\ C_2,\\ &{} \pi /4,\ on\ C_3,\\ &{} -\pi /4\ on\ C_4, \end{array}\right. }\nonumber \\ \gamma _b= & {} {\left\{ \begin{array}{ll} &{}\pi /4,\ on\ C_1,\\ &{} -\pi /4,\ on\ C_2,\\ &{} -3\pi /4,\ on\ C_3,\\ &{} -\pi /4\ on\ C_4, \end{array}\right. } \gamma _b = {\left\{ \begin{array}{ll} &{}\pi /4,\ on\ C_1,\\ &{} 3\pi /4,\ on\ C_2,\\ &{} 5\pi /4,\ on\ C_3,\\ &{} 3\pi /4\ on\ C_4. \end{array}\right. } \end{aligned}$$
(4.4)

For the interior problem on \(E_4\), discarding rotation and reflection symmetries, there is simply one D state and one R state. This is not the case for the exterior problem on \(E_4^C\), due to the choice of the \(\gamma ^*\). In Figure 10, we plot scalar order parameter and planar director of the reduced \(\textbf{P}\)-tensors defined in (4.1) for the different choices of \(\gamma _b\) defined above, for two different choices of \(\gamma ^*\). For \(\gamma ^* = n\pi \), there are three distinct classes of solutions: D1, D2 and R with \(\gamma ^* = n\pi \), and for \(\gamma ^* = \pi /4 + n\pi \), there are four distinct classes of solutions: St, D, R1, and R2. The label St originates from the work in Phillips and Rey (2011), which is the abbreviation of string defect mode.

We can repeat the same combinatorial arguments for \(K=6\), to get simple estimates for the number of stable states on \(E^C_6\), in the \(\lambda \rightarrow \infty \) limit. In Han et al. (2020), the authors study the same problem on \(E_6\) and report three rotationally equivalent classes of solutions - Para, Meta and Ortho distinguished by the locations of the splay vertices. This nomenclature originates from organic chemistry to indicate the position of non-hydrogen substituents on a hydrocarbon ring (benzene derivative), and one can identify these substituents with splay or bend vertices. Following the same recipe as in Figure 10, we compute the limiting profiles for \(\left( {\begin{array}{c}6\\ 2\end{array}}\right) = 15\) different choices of \(\gamma _b\) in (4.1), for \(\gamma ^* = \pi /6+n\pi \). Discarding rotation and reflection symmetries, we recover at least six different classes of limiting profiles (as opposed to three on \(E_6\)), plotted in Fig. 11. On \(E^C_6\), there are two classes of Para solutions - Para1 and Para2 with two diagonally opposite bend vertices. Para2 is distinguished in the sense that the diagonal connecting the bend vertices is orthogonal to the direction \(\textbf{n}^* = (\cos \gamma ^*, \sin \gamma ^*)\). There are two classes of Meta states: Meta1 and Meta2, for which the two bend vertices are separated by one vertex, and two further classes: Ortho1 and Ortho2 with two “adjacent” bend vertices connected by an edge. The unique \(\textbf{P}\)-solution for small enough \(\lambda \), is Para2 with two diagonally opposite bend vertices on the diagonal orthogonal to the fixed director at infinity; see Fig. 11 (b). Para2 is stable for all \(\lambda \).

Fig. 11
figure 11

Left: distinct classes of solution of (4.2) which cannot be recovered from each other by reflection or rotation on the exterior of hexagon domain with \(\gamma ^* = \pi /6+n\pi \), (a1) Para1, (a2) Para2, (b1) Meta1, (b2) Meta2, (c1) Ortho1, (c2) Ortho2. Right: the energy plot of solutions in (2.8) versus \(\lambda ^2\)

Fig. 12
figure 12

Left : plots for D1 with \(\lambda ^2 = 5\), D2 and R with \(\lambda ^2 = 100\), the solutions of the Euler–Lagrange equations in (2.8) on \([-10,10]^2\backslash E_4\) with \(\gamma ^* = n\pi \) (the configuration of D is zoomed in). Right: the plot of the rLdG energy in (2.3) for states on the left versus \(\lambda ^2\)

The limiting profiles in Sects. 3 and 4 give us excellent initial conditions for finding solution branches of the Euler–Lagrange equations (2.8). We numerically compute the solutions of (3.1) on domain \([-10,10]^2\backslash E_K\), which is used to approximate the unbounded domain, \(E_K^C\), with the Dirichlet boundary condition \(\textbf{P} = \textbf{P}^*\) in (2.7) on x or \(y= \pm 10\). In Figures 10 and 11, we numerically compute the limiting profiles of \(\gamma \) for (4.2) on \([-10,10]^2\backslash E_K\). These limiting profiles of \(\gamma \) define the limiting profiles, \(\textbf{P}\) via the relation \((P_{11},P_{12}) = \frac{B}{2C}(\cos (2\gamma ),\sin (2\gamma ))\). In Figures 111213, we use these limiting profiles as initial conditions to compute solution branches for different values of \(\lambda ^2\); we do not track bifurcations and may have missed solution branches but simply focus on energy comparisons between the numerically computed solution branches. We numerically compute the solutions of the Euler–Lagrange equations (2.8), subject to the Dirichlet boundary conditions (2.5), which are necessarily critical points of (2.3). We take the solution of (3.1) ((4.2)) on \([-10,10]^2\backslash E_K\) as the initial conditions for finite but small (large) \(\lambda ^2\). We perform an increasing \(\lambda ^2\) sweep for the unique branch (in the \(\lambda \rightarrow 0\) limit) and decreasing \(\lambda ^2\) sweep for the distinct diagonal- and rotated-branches. For small \(\lambda \), a small step-size is chosen and a larger step-size is chosen for large \(\lambda \), and this is somewhat ad hoc. Except for the St branch in Fig. 13, we stop tracking the branch if the smallest eigenvalue approaches zero, i.e. near a bifurcation point. For all the numerical simulations in this manuscript including the calculation of the limiting profile on a unit disc in Figures 36, and 8 , we use the popular open-source computing software FEniCS for solving partial differential equations, using the finite element method Wells et al. (2012), which allows us to solve the weak form of the Euler–Lagrange equations or the Laplace equations, in the first order Lagrange element function space by using Newton’s method, with tolerance \(10^{-13}\) and mesh size less than 1/256. The convergence may be highly sensitive to the choice of initial condition. The energy of a solution is calculated according to the integral in (2.3). We check the stability of a solution by checking the sign of the smallest eigenvalue of the Hessian. All computed solution branches are stable in the sense that the corresponding second variation of (2.3)

$$\begin{aligned} \partial ^2F_{\lambda }[\eta ] = \int _{[-10,10]^2\backslash E_4}|\nabla \eta |^2 + \frac{\lambda ^2}{4}\left( |\textbf{P}|^2-\frac{B^2}{2C^2}\right) |\eta |^2+\frac{\lambda ^2}{2}\left( \textbf{P}\cdot \eta \right) ^2 dx \end{aligned}$$
(4.5)
Fig. 13
figure 13

Left : plots for St with \(\lambda ^2 = 5\), D, R1 and R2 with \(\lambda ^2 = 100\), the solutions of Euler–Lagrange equations in (2.8) on \([-10,10]^2\backslash E_4\) with \(\gamma ^* = \pi /4+n\pi \) (the configuration of St is zoomed in). Right: the plot of energy in (2.3) for the states on the left versus \(\lambda ^2\), where the solid lines represent stable solution branches and dashed lines correspond to unstable branches. We plot a zoom near the end of R, and the branch R doesn’t connect with unstable St branch

is strictly positive according to our numerical computations. The smallest real eigenvalue of the Hessian is computed by the SLEPc eigenvalue solver with the algorithm krylov-schur, and the tolerance is \(10^{-15}\). [https://fenicsproject.org/olddocs/dolfin/1.3.0/python/programmers-reference/cpp/la/SLEPcEigenSolver.html].

In Figure 12, we study stable solution branches of the Euler-Lagrange equations (2.8) on \(E^C_4\), with \(\gamma ^* = n\pi \) and the prescribed boundary conditions on \(\partial E_4\), for different values of \(\lambda ^2\). For small \(\lambda ^2\), there is a unique solution of the system (2.8) on the domain \([-10,10]^2\backslash E_4\), subject to the boundary conditions (2.5) and constraint \(\gamma ^* = n\pi \), labelled as the D1 state, with two bend vertices along the diagonal on \(x = 0\). This solution branch is computed using the limiting profile in Sect. 3, i.e. the solution of (3.1) on \(E^C_4\) and then using continuation methods for large \(\lambda \). This solution branch remains stable for all \(\lambda ^2>0\). When \(\lambda ^2\) is large enough, we observe the stable D2 with two bend vertices along the diagonal on \(y=0\), and the stable R solution. The D2 represents two energetically degenerate diagonal states which are related by reflection about \(x = 0\) or \(y = 0\) (see the first line of Fig. 10(a)) and the R represent four rotated states (see the second line of Fig. 10(a)). The energy of D1 solution is always lower than the D2 solution and the rotated solutions have higher energies than the diagonal solutions. We do not investigate the connections between solution branches in this manuscript further.

The case of \(\gamma ^* = \pi /4+n\pi \) is different (see Fig. 13). The unique \(\textbf{P}\)-solution for small \(\lambda \) is the St-branch with two line defects near opposite square edges. This St-branch exists for all \(\lambda ^2>0\). As \(\lambda ^2\) increases, St loses stability and bifurcates into the stable D, diagonal solution with two diagonally opposite bend vertices. When \(\lambda ^2>8.2\), the stable R1 solution branch exists. The R1 solution has two bend vertices connected by the edge of square, and this edge is perpendicular to \((\cos \gamma ^*, \sin \gamma ^*)\). There is a stable R2 branch, where the R2 has two bend vertices along the square edge, which is parallel to \((\cos \gamma ^*, \sin \gamma ^*)\). This solution branch only exists for \(\lambda ^2\) large enough. The energy of the diagonal solution is always lower than R1, which in turn is energetically preferable to the R2 solution branches. The families of solutions, labelled by D, R1 and R2, have two members each. The two states in D or R1 or R2 class are related by reflection about \(y = -x\) (see Fig. 10(b)). The diagram with \(\gamma ^* = \pi /4+n\pi \) is similar to the diagram on \(E_4\) reported in Han et al. (2020), with the St-branch on \(E_4^C\) being analogous to the WORS on \(E_4\), and the D-branch on \(E_4^C\) being analogous to the D on \(E_4\), the exception being that the R1 and R2 solutions are energetically degenerate on \(E_4\) but not on \(E_4^C\) due to the symmetry-breaking condition at infinity. Similar energy computations are performed for the Para, Meta, Ortho solution branches in Figure 11 on \(E_6^C\), for illustrative purposes.

5 Conclusion

This paper focuses on illustrative studies of nematic equilibria outside 2D regular polygonal holes, \(E_K\) with K edges, focusing on the effects of the shape and size of the polygon in a reduced LdG framework. The equilibria are modelled in terms of a reduced order parameter, \(\textbf{P}\) with two degrees of freedom, that are minimisers of a rLdG free energy (2.3), subject to homeotropic boundary conditions on the hole boundary and a fixed far-field condition, captured by \(\gamma ^*\). We focus on two distinguished limits: the \(\lambda \rightarrow 0\) limit relevant for small holes, and the \(\lambda \rightarrow \infty \) limit relevant for large holes. We study the limiting problem for \(\lambda =0\) comprehensively, using methods from complex analysis and Ginzburg–Landau theory. The limiting profile is unique for \(\lambda =0\), has an isolated interior point defect for \(K=3\), two interior point defects for \(K>4\) and for \(K=4\), the limiting profile has either line defects or no interior defects depending on \(\gamma ^*\). In the \(\lambda \rightarrow \infty \) limit, we have multiple stable equilibria, i.e. at least \(\left( {\begin{array}{c}K\\ 2\end{array}}\right) \) stable equilibria on \(E_K^C\), where K is the number of edges of the polygon hole. This is certainly an underestimate, based on two assumptions with regards to minimal topology. Compared to the interior problem studied in Han et al. (2020), where the authors study minimisers of the rLdG energy, (2.3), on \(E_K\) with tangent boundary conditions, we find that there are more possibilities for defect sets and multistability for the exterior problem. For example, the far-field condition, \(\gamma ^*\) tunes the location of defects for limiting profiles and for \(K=4\), \(\gamma ^*\) dictates the dimensionality of the defect set for limiting profiles. Further, \(\gamma ^*\) strongly enhances multistability for the exterior problem compared to the interior counterpart in Han et al. (2020). For example, in the \(\lambda \rightarrow \infty \) limit, there are only two competing equilibria on \(E_4\) but there can be up to four different classes of competing equilibria on \(E_4^C\), depending on the choice of \(\gamma ^*\). In this respect, our work provides good information about how K and \(\gamma ^*\) can tune the existence, locations and dimensionality of defects on \(E_K^C\), and how this can be used to tune multistability. Our work does not uncover solution landscapes for this exterior problem and does not give information about the connectivity of solution branches, i.e. there is no systematic bifurcation analysis. This is much needed and could be tackled using high-index optimisation shrinking dimer methods, as in Han et al. (2021). Such a study would be crucial for understanding pathways between distinct equilibria, the energy barriers and understanding the dynamics of such exterior problems on energy landscapes. In fact, a complete study of solution landscapes can give strong insight into the differences between the interior and exterior problem. Other natural generalisations include addition of the effects of elastic anisotropy (see Han et al. (2021)), other types of boundary conditions on \(\partial E_K\) and studies of arrays of polygonal holes and creation of linked defect structures. Our work provides a foundation for further detailed, comprehensive studies and our overarching goal is to propose universal theoretical frameworks for solution landscapes of confined partially ordered systems, with multiple order parameters.