We start from Onsager's model [START_REF] Onsager | The effects of shape on the interaction of colloidal particles[END_REF] adapted for rigid cylindrical discs with diameter D and thickness L exposed to some effective orientational potential U e imparted by the surface anchoring and elastic deformation each disc experiences within the molecular liquid crystal (LC). The Euler-Lagrange equation for the orientation distribution f c (ω c • n m ) of the disc normal vectors ω c around the molecular director n m reads [START_REF] Kayser | Bifurcation in onsager's model of the isotropic-nematic transition[END_REF] :

ln[f c (ω c • n m )] = λ EL -2ρ dω c B 2 (ω c , ω c )f c (ω c • n m ) -U e (ω c • n m ) (1)
with λ EL a multiplier ensuring normalisation of f c on the unit sphere and ρ = N/V the disc number density defined as the ratio of the number of discs and the volume within which they are located. Without loss of generality, all (free) energies are implicitly normalised in terms of the thermal energy k B T with k B Boltzmann's constant and T temperature. Even though Onsager's theory is known to be quantitatively inaccurate for discs [START_REF] Onsager | The effects of shape on the interaction of colloidal particles[END_REF] due to the neglect of virial coefficients beyond the second one B 2 , it should provide useful qualitative insights into the disorder-order transitions of discotic particles within a molecular host. When the disc interactions are steeply repulsive, the second-virial coefficient B 2 equals half the excluded volume v ex between the particles (Supplementary Fig. S1a and Extended Data Fig. 1k). For cylinders the second-virial coefficient takes the following form [START_REF] Onsager | The effects of shape on the interaction of colloidal particles[END_REF] :

B 2 (ω c , ω c ) = 1 2 v ex (ω c , ω c ) = π 4 D 3 + L 2 D |ω c × ω c | + LD 2 π 4 + E(|ω c × ω c |) + π 4 |ω c • ω c | (2)
where E(x) is the complete elliptic integral of the second kind. For thin discs we consider in this study (L/D 1) it suffices to retain only the dominant first term in the excluded volume. It is customary to define a dimensionless disc concentration c = ρD 3 π 2 /16 so that the Euler-Lagrange equation can be recast as follows: We parameterise the orientational phase space of each disc in terms of a polar angle 0 < θ < π, with cos θ = ω c • n m (-1 < cos θ < 1) and an azimuthal angle 0 < ϕ < 2π describing the particle orientation with respect to the molecular director n m (Supplementary Fig. S1a and Extended Data Fig. 1k).

ln[f c (ω c • n m )] =λ EL - 8 π c dω c |ω c × ω c |f c (ω c • n m ) -U e (ω c • n m ) (3)

B. Disc surface anchoring potential

The effective potential imparted by the host medium dictates that the disc normals preferentially adopt a certain equilibrium angle θ e with respect to the molecular director. We will now demonstrate that there is an orientational potential of mean force which reads up to quartic order in the deviation from the equilibrium anchoring angle:

U e (θ) ∼ ε 2 (θ -θ e ) 2 + ε 3 (θ -θ e ) 3 + ε 4 (θ -θ e ) 4 (4)
The amplitudes ε n can be established from the free energy of the surface anchoring and elastic deformation incurred by a single disc embedded in a molecular LC. Following the analysis of Refs. 53,54 we define the surface anchoring energy in terms of a Mexican-hat type potential:

W (θ) = A 2 cos 2 θ + W 2 4 cos 4 θ (5)
where A combines a molecular-field and dielectric response due to the presence of electrostatic double layers surrounding each disc [START_REF] Barbero | Selective ions adsorption and nonlocal anchoring energy in nematic liquid crystals[END_REF] :

A = W 1 -2 0 ∆ E 2 D ξ D (6)
The surface anchoring potential Eq. ( 5) depends on the anchoring energy densities W 1 and W 2 , the dielectric anisotropy ∆ > 0 of the molecular LC, the typical electric field strength E D within the double layer, and the Debye screening length ξ D . The molecular contribution is governed by the interaction between the surface-coating of the discs and the molecular LC and favours planar anchoring (θ = π/2) whereas the electrostatic contribution promotes homeotropic anchoring (θ = 0). The optimum anchoring angle, defined as θ e , follows from minimising the surface anchoring energy and reveals three anchoring scenarios:

θ e =              π 2 , A W 2 ≥ 0 (planar) arccos |A| W 2 , A W 2 < 0 (conical) 0, A W 2 < -1 (homeotropic) (7)
Inserting cos 2 θ e = |A|/W 2 into Eq. ( 5) and some algebraic manipulation allow us to recast the surface anchoring free energy per disc as follows [START_REF] Liu | Biaxial ferromagnetic liquid crystal colloids[END_REF] :

F s (θ) N ∼ πD 2 |A| 8 cos 2 θ e (cos 2 θ -cos 2 θ e ) 2 (8)
The elastic free energy associated with finite-strength conically degenerate surface anchoring reads [START_REF] Brochard | Theory of magnetic suspensions in liquid crystals[END_REF] :

F e (θ) N ∼ 2KD(θ -θ e -∆θ) 2 (9)
where K is the one-constant elastic modulus of the molecular LC and ∆θ is the angle to which the director deviates away from the equilibrium angle θ e at the disc surface. Following Liu et al. [START_REF] Liu | Biaxial ferromagnetic liquid crystal colloids[END_REF] , we substitute θ = θ e + ∆θ into the surface free energy Eq. ( 8) and minimise the total free energy F tot = F s + F e with respect to ∆θ. This gives up to linear order:

∆θ ≈ 4K(θ -θ e ) 4K + πDW 2 cos 2 θ e sin 2 θ e (10) Director distortions near the disc surface (Fig. 1f andi) are weak due to the fact that the thickness of the discotic inclusion (∼ 10 nm) is much smaller than the typical surface extrapolation length 57 K/W 2 ∼ 100 nm. The latter may be estimated by taking the experimentally measured values for the average elastic constant [START_REF] Liu | Biaxial ferromagnetic liquid crystal colloids[END_REF]58 of 5CB K = 6 pN and polar anchoring coefficient [58][59][60] W 2 = 10 -5 J/m 2 . Reinserting the above expression back into the total free energy F tot and expanding up to quartic order in the deviation from the optimal anchoring angle θ -θ e leads to the potential of mean force Eq. ( 4). The amplitudes ε n featuring in Eq. ( 4) depend on the anchoring angle itself and on the balance between the two principal energy scales: the bulk elastic one K = KD and the surface anchoring energy per particle W2 = W 2 D 2 . In explicit form they read:

ε 2 ∼ π 2 W2 w 0 sin 2 (2θ e ) ε 3 ∼ -8π W2 w 3 0 sin(4θ e ) ε 4 ∼ 16π 3 W2 w 4 0 (7 cos(4θ e ) -1) (11)
The prefactor w 0 = (4 + π 4 W2 K sin 2 (2θ e )) -1 depends itself upon the equilibrium anchoring angle θ e but remains finite at all values. While both ε 2 and ε 3 vanish in the limit of homeotropic anchoring (θ e → 0), or when approaching planar anchoring (θ e → π/2), the quartic contribution remains finite ε 4 ∼ π W2 /8 in either limit. Given that both the elastic and surface anchoring energy are considerably large ( K ∼ O(10 3 ) and W2 ∼ O(10 4 )) we infer that the amplitudes strongly exceed the thermal energy (ε n 1) and that the disc orientations will generally be strongly confined around θ e (Supplementary Fig. S1 and Extended Data Fig. 1l). Even though both K and W 2 depend in principle on temperature, their variation will be very slim across the narrow temperature range probed in experiment and hardly impact the anchoring amplitudes. Close inspection of the expressions Eq. ( 11) leads us to infer the following principal features (Supplementary Fig. S2):

• As ε 4 grows very large near 0 and π/2, the distribution should become narrower upon approaching these limits. This is corroborated by the significant narrowing of the experimental distributions upon approaching these angles (Fig. 4g and Extended Data Fig. 1j).

• The cubic term ε 3 is nonzero and the distribution is asymmetric around θ e which reflects symmetry-breaking monoclinic order within the disc subsystem (Figs. 4g,h and 5d). This will be further discussed in paragraph C-3.

C. Orientational order parameters

In this paragraph we outline our approach to describing the orientational order of the discs. In doing so, we note that while we mostly consider the nematic phases featuring in the phase diagram (Fig. 2), we also quantify orientational order in some of the positionally ordered hybrid liquid crystals we observe. We do not consider any order parameters describing long-range positional correlations related to smectic or columnar order. We further reiterate that the colloidal discs are uniaxial and their orientation is described by a single unit vector ω c (Supplementary Fig. S1a).

1. Uniaxial and biaxial disc order with respect to n m

Going back to Eq. ( 3) and expanding the kernel |ω c × ω c | up to lowest order in symmetryadapted functions using the addition theorem for spherical harmonics we find 61 :

|ω c × ω c | ∼ π 4 - 5π 32 [P 2 (t)P 2 (t ) + 1 12 P 2 2 (t)P 2 2 (t ) cos 2(ϕ -ϕ)] (12)
with t = cos θ. From this we obtain a self-consistent expression for the orientation probability of the discs: The symmetry functions are given by P 2 (ω c ) = 3 2 cos θ 2 -1 2 , and D 2 (ω c ) = sin 2 θ cos 2ϕ. The normalisation constant N ensures that dω c f c (ω c ) = 1. At infinite dilution, anchoring may occur at an angle θ e or π -θ e with equal probability (reflecting the basic apolarity of the hybrid system), but particle crowding at elevated disc concentration will force the system to select only one colloidal aligning direction which we constrain at the interval 0 ≤ θ e ≤ π/2.

f c (ω c ) =N exp 5 4 cSP 2 (ω c ) + 15 16 c∆D 2 (ω c ) - 4 n=2 ε n (θ -θ e ) n (13)
We define nematic order parameters compatible with those measured in experiment:

S = dω c f c (ω c )P 2 (ω c ) ∆ = dω c f c (ω c )D 2 (ω c ) (14)
where S measures uniaxial and ∆ describes orthorhombic biaxial nematic order 62-64 of the discs defined along the molecular director n m . The above expressions coincide with those in Eq. 1 of the main text. For completeness, we note that both order parameters naturally feature in a nematic ordering tensor 58 diagonalised along the molecular frame n m :

Q =S 3 2 n m ⊗ n m - I 2 + ∆[(n m × n ⊥ ) ⊗ (n m × n ⊥ ) -n ⊥ ⊗ n ⊥ ] (15)
with ⊗ denoting a dyadic product. This tensor describes orientational order of the colloids measured within the orthonormal reference frame spanned by n m .

Geometric monoclinic order parameter

In order to distinguish the monoclinic phases from the orthorhombic and uniaxial ones we complement the above parameters with a monoclinic order parameter that is non-zero only for the monoclinic hybrid system. A basic measure for the degree of monoclinicity can be established from considering a superposition of two uniaxial subsystem described by the ordering tensor:

Q (m,c) =S m,c 3 2 n m,c ⊗ n m,c - I 2 (16)
where the scalar order parameters S m and S c quantify the degree of uniaxial nematic order of each component (colloids or molecules) along their respective directors. Note that S c = S, except for the uniaxial hybrid system where both directors co-align and S c = S.

Let us now consider an oblique, non-orthogonal frame {n c , n m , n ⊥ } with n ⊥ = (n m × n c )/ sin θ ne denoting the 2-fold rotation axis of the monoclinic tripod representing the principal directions of alignment each component (Supplementary Fig. S3e). The monoclinic symmetry of the hybrid molecular-colloidal system is reflected by a geometric superposition of the uniaxially ordered subsystems with obliquely oriented principal ordering directions. More specifically, one may define the anti-symmetric component T of the product of both uniaxial ordering tensors 65 :

T = 1 2 Q (m) Q (c) -Q (c) Q (m) = M 2 (n m ⊗ n c -n c ⊗ n m ) (17)
The last expression features a pseudovector with a scalar monoclinic order parameter M given by 65 :

M ∼ S c S m sin(2θ ne ) (18)
with θ ne the angle between the colloidal director n c and the molecular one n m (Supplementary Fig. S3) 89 . The order parameter naturally vanishes in the limit of uniaxial order (θ ne = 0) as well as for the orthorhombic nematic (θ ne = π/2). Trivially, it also disappears when either or both components relax toward an isotropic fluid (S m and/or S c → 0).

In reality, the subtle orientational coupling between the discs and the surrounding molecular liquid crystal imprints weakly orthorhombic and monoclinic features onto the disc orientational probability as well as onto the molecular host (green zones in Extended Data Fig. 2 and Supplementary Fig. S3) that will modulate the prefactor in Eq. ( 18). A full account of these effects through a rigorous tensorial description 66 that encodes the point-group symmetry of each subsystem clearly goes beyond the scope of the present paper 65 . For the sake of simplicity, we assert that the order in Eq. ( 18) is given by a product of uniaxial components S c S m ∼ O(1) that we may set to unity without loss of generality. We then arrive at the following geometric monoclinic order parameter:

M = sin(2θ ne ) (19)
which solely depends on the oblique angle θ ne between the individual directors (Eq. 1 in the main text and Fig. 4e).

Skewness of orientational distributions

The intrinsic skewness of the orientational distribution of the colloidal discs around their mean anchoring angle θ e can be established by eliminating the azimuthal dependency of f c in Eq. ( 13) to obtain the distribution of polar angles. Integrating Eq. ( 13) over ϕ and ignoring the normalisation constant we find:

f c (θ) ϕ ∼ exp 5 4 cSP 2 (cos θ) -4 n=2 ε n (θ -θ e ) n I 0 (c∆ sin 2 θ) (20)
with I 0 a modified Bessel function. The intrinsic "obliqueness" or skewness of f c around its peak value is reflected by the moment coefficient of skewness associated with the orientation distribution of the colloids 67 :

Sk c = dω c f c (ω c )(θ -θ e ) 3 dω c f c (ω c )(θ -θ e ) 2 3/2 (21)
The skewness vanishes in the uniaxial nematic phase (θ e = 0) as well as in the orthorhombic nematic (θ e = π/2), but remains finite in the conical anchoring regime (0 < θ e < π/2) where the cubic term ε 3 in the anchoring potential Eq. ( 4) is nonzero. The skewness of the angular distributions can be gleaned from Fig. 4g,h).

The nonzero skewness of the disc orientational probaility reinforces the notion that monoclinic order is not only present on the hybrid level due to a superposition of two uniaxially ordered subsystem characterised by an oblique mutual director angle [Eq. ( 17)], and that it is also imparted at the disc level through the n m -distorted isosurface zones surrounding the disc surface (Supplementary Fig. S3) as well as by the intrinsic asymmetry of the disc polar angular probability around its mean value.

D. Temperature dependence of the equilibrium anchoring angle

We may recast Eq. ( 6) in a more convenient form that enables us to describe the temperaturedependence of the anchoring angle θ e and amplitudes ε n . We assume that the dielectric contribution to surface anchoring remains mostly insensitive to weak temperature variations. The molecular surface anchoring strengths, however, are known to decrease monotonically with temperature 68,69 . We postulate a simple exponential decay for the surface anchoring densities:

W n = W 0 n exp[-a(T -T r )], n = 1, 2 (22)
with W 0 n the off-set values at room temperature T r (say 298K) and a > 0 an unknown decay rate identical to both parameters. The surface anchoring ratio then attains the following temperaturedependence:

A W 2 = W 0 1 W 0 2 - A D W 0 2 exp[a(T -T r )] (23)
with A D = 2 0 ∆εE 2 D ξ D the dielectric contribution. Since W 0 1 /W 0 2 and A D are both positive, we find that A/W 2 is positive at low temperatures (planar anchoring, θ e = π/2), but becomes negative at higher temperature leading to conical anchoring. The critical temperature where this happens is given by:

T c -T r = a -1 ln W 0 1 A D ( 24
)
The anchoring angle scales with temperature in the following way:

θ e =      π 2 , T ≤ T c arccos W 0 1 W 0 2 1 -e a(T -Tc) , T > T c (25)
where we use the convention 90 • < θ e < 180 • so that cos θ e < 0. If we identify T 0 > T c with the temperature where θ e = 180 • (cos θ e = -1) we find:

T 0 -T c = a -1 ln 1 + W 0 2 W 0 1 (26)
Using this result to eliminate a we find after some rearranging:

cos θ e =        0, T ≤ 1 W 0 1 W 0 2 1 -1 + W 0 2 W 0 1 T -1 T0 -1 , T > 1 (27)
in terms of a renormalised temperature scale:

T = T T c ( 28
)
The variation of the conical anchoring angle then only depends on a single parameter, namely the ratio of the molecular anchoring strengths W r = W 0 1 /W 0 2 which should be independent of the disc diameter. We take T 0 = 33.5 ± 0.1 • C and T c = 29.5 ± 0.1 • C. A fit to the experimental data yields W r = 0.06±0.03 and, taking Eq. ( 26), a = (0.7±0.1)K -1 . In order to gauge the scaling of θ e close to the homeotropic-to-conical anchoring transition temperature T 0 we write T -1 T0 -1 = -( T0 -T ) T0 -1 + 1. Retaining the leading order contribution for small T0 -T we obtain:

θ e (T ) ∼ 180 • π ζ(T 0 -T ) 1 2 (29) with ζ 2 = 2a (1 + W r ).
According to Landau theory, this confirms the presence of a continuous anchoring transition at the colloidal surface in line with previous reports on anchoring transitions at confining surfaces [68][69][70] .

E. Transition from isotropic to uniaxial nematic phase

At temperatures where the molecular LC is isotropic the discs undergo a disorder-order transition to a conventional uniaxial (U) nematic phase as the interactions are unaffected by surface anchoring or elastic forces imparted by the medium. However, the discs carry surface charges and electrostatic repulsion plays a role in determining the critical concentration for the transition.

A simple recipe to account for the near-field electrostatic interactions is to assume an effective disc shape dictated by the extent of electrostatic screening which is kept more or less constant in the experimental situation. Following van der Beek 71 we introduce an effective disc diameter and thickness:

L ≈ L + 2ξ D D ≈ D + 2ξ D ( 30
)
which leads to an effective aspect ratio L/ D being larger than the bare ratio L/D ∼ 0.01; taking a typical Debye screening length ξ D ∼ 100 nm one finds L/ D ∼ 0.1. Keeping only the leading order term in the excluded volume we find that the Helmholtz free energy per particle of the fluid takes a simple form (ignoring irrelevant constant factors that are equal in both states):

F N ∼ ln c + ln f c fc + cp( φ) 4 π |ω c × ω c | fc (31)
where p(φ) = (1-3 4 φ)/(1-φ) 2 is a Parsons-Lee correction factor 72-74 that approximately accounts for higher-order steric disc interactions through an effective colloid packing fraction φ = 4 π L D c. We use a simple Gaussian representation 75 f c (θ) ∼ exp(-α 2 θ 2 ) (including its mirror form f c (π -θ)) to describe the fluctuations of the polar angle about the nematic director and keep only the leading order asymptotic contributions of the angular averages • fc valid for α → ∞. For the isotropic fluid, we simply have f c = (4π) -1 and we find that Onsager-Parsons theory predicts the transition to happen at c I ≈ 2.21 and c U ≈ 2.94 (Fig. 2e). A sample Mathematica code can be found in Appendix A. We remark that larger effective aspect ratios could be reached in experiment where the Debye screening length may go up to about 500 nm.

F. Transition from uniaxial to monoclinic and orthorhombic nematic phases Below the (uniaxial) nematic-isotropic transition temperature of the molecular host, surface anchoring sets in and the anchoring amplitude ε will be much larger than the thermal energy. The conically degenerate surface anchoring conditions the discs experience then facilitate a crossover from a uniaxial nematic (U) to a monoclinic (M) or orthorhombic (O) nematic fluid. We reiterate that the distinction between the three states lies principally in the angle θ ne between the colloidal and molecular directors which is zero (θ ne = 0 : U), oblique (0 < θ ne < π/2 : M), or right (θ ne = π/2 : O). We may then envisage two types of orientationally ordered fluids. The first is a uniaxial one where the azimuthal angles ϕ are randomly distributed while fluctuations of the polar angle around the equilibrium anchoring angle θ e are assumed to be asymptotically small:

f (U ) c (ω c ) ∼ δ(θ -θ e ) sin θ e 1 2π (32)
with δ(x) a Dirac delta distribution around the mean angle θ e , and angular phase space parame-

terised via θ ∈ [0, π] and ϕ ∈ [-π, π].
At larger disc concentration the discs tend to align and form a low-symmetry monoclinic or orthorhombic nematic fluid. In order to describe these phases, we assume that the azimuthal angle is described by a Gaussian distribution around ϕ = 0 while the polar angles are sharply distributed along the colloidal director indicated by the angle θ ne . Then one may deduce from Eq. ( 13) that for the monoclinic/orthorhombic (MO) nematic fluid:

f (M O) c (ω c ) ∼ δ(θ -θ ne ) sin θ ne β(sin θ ne ) 2 2π 1 2 exp - β 2 (ϕ sin θ ne ) 2 (33)
Here, β -1 is proportional to the strength of the azimuthal angular fluctuations. Because these fluctuations are assumed small, it is required that β 1 so that the integration domain may safely be replaced by ϕ ∈ -∞, ∞ . Clearly, in the asymptotic limit of strongly aligned discs the colloidal director angle is expected to be very close to θ e . Expanding the second Legendre polynomial in the exponent of Eq. ( 20) up to quadratic order in θ -θ e and ignoring the non-Gaussian terms we formulate the following asymptotic limit:

δ(θ -θ ne ) = lim ε2 →∞ ε2 π 1 2 e -ε 2 (θ-θe) 2 (34)
where ε2 = ε 2 + 15 8 cS cos(2θ e ) combines the single-particle surface anchoring contribution with a collective term arising from disc-pair correlations. The deviation of the colloidal director from the anchoring direction is strictly nonzero for conical anchoring angles (0 < θ e < π/2), but vanishes for asymptotically strong anchoring:

(θ e -θ ne ) = lim ε2 →∞ cS sin θ e cos θ e 2ε 2 (35)
For practical purposes, we may simply consider θ ne identical to θ e for the current system where disc orientations are chiefly governed by single-particle surface anchoring rather than by collective, excluded-volume effects. The Helmholtz free energy reads as before:

F N ∼ ln c + ln f c fc + cp( φ) 4 π |ω c × ω c | fc + δF el (36)
Here, δF el is a fluctuation term quantifying the impact of the nemato-elasticity of the molecular LC onto the ordering of the discs. Using the Gaussian representation Eq. ( 34) for the orientational distribution with θ ne = θ e we obtain for the orientational entropy in the uniaxial state:

ln f (U ) c f (U ) c ∼ lim ε 2 →∞ ln ε 1 2 2 π -1 2 - 1 2 ε 1 4 2 -ln 2π (37)
For the monoclinic/orthorhombic state we find:

ln f (M O) c f (M O) c ∼ ln f (U ) c f (U ) c + 1 2 ln β(sin θ ne ) 2 2π - 1 2 (38)
The first term attains the same asymptotic limit as the one for the uniaxial state and will therefore be irrelevant for the phase transition. The excluded volume entropy for the uniaxial fluid is proportional to a double orientation averaging of |ω c × ω c | which, upon taking the limit ε 2 → ∞ amounts to fixing the polar angle at θ = θ e so that:

|ω c × ω c | f (U ) c ∼ π -π d∆ϕ 2π 1 -[(cos θ e ) 2 + (sin θ e ) 2 cos ∆ϕ] 2 (39)
where the relative azimuthal angle ∆ϕ = ϕ -ϕ is a random variable. The integral can be solved analytically:

|ω c × ω c | f (U ) c = 2 π | sin θ e | (1 + arctanh(sin θ e ) cos θ e cot θ e ) (40)
For the monoclinic/orthorhombic state we exploit the fact that ∆ϕ 1 and we find that the typical strength of the azimuthal angular fluctuations exhibits the following scaling:

|ω c × ω c | f (M O) c ∼ | sin θ ne | β(sin θ ne ) 2 2π ∞ -∞ dϕ ∞ -∞ dϕ e -1 2 β(sin θne) 2 (ϕ 2 +ϕ 2 ) |∆ϕ| ∼ 2 √ πβ (41)
Minimisation of the free energy with respect to β yields a simple quadratic scaling of β with disc concentration:

β ∼ 64 π 3 (cp) 2 (42)
showing that the azimuthal angular fluctuations of the disc normals with respect to n m become progressively weaker at increased disc concentration. The dimensionless pressure Π and chemical potential µ for the uniaxial phase follow from standard thermodynamic derivatives of the free energy. Using Eq. ( 40) we find:

Π (U ) ∼ c + c 2 Π p 4 π |ω c × ω c | f (U ) c + cδF (U ) el µ (U ) ∼ ln c -ln 2π + cµ p 8 π |ω c × ω c | f (U ) c + 2δF (U ) el (43)
where Π p = p + φ dp dφ and µ p = p + 1 2 φ dp dφ emerge from the Lee-Parsons rescaling. Noting that

4 π |ω c × ω c | f (M O) c
∼ 1/cp one finds the following asymptotic expressions for the monoclinic/orthorhombic nematic phases:

Π (M O) ∼ c + c p Π p + cδF (M O) el µ (M O) ∼ ln c + 1 2 ln β(sin θ e ) 2 2π - 1 2 + 2 p µ p + 2δF (M O) el (44)
Results for the original second-virial approximation are easily recovered by setting all scaling factors p, Π p and µ p to unity. The pressure then interpolates between the ideal gas pressure (Π ∼ c)

and that of a conventional uniaxial nematic (Π ∼ 3c) at infinite particle alignment in the limit of vanishing packing fraction 75 . We remark that the fixed-θ e constraint used in the thermodynamic analysis likely leads to an overestimation of the critical concentration for the uniaxial-monoclinic (orthorhombic) nematic transition (sample Mathematica code in Appendix B). The agreement with the experimental diagram could therefore be further improved by including i) orientational fluctuations about θ e [cf. Eq. ( 13)] and ii) electrostatics-driven finite-thickness corrections to the excluded volume [cf. Eq. ( 2)]. In view of the technical burden these incur, both topics go beyond the scope of the present tractable theory.

G. Nematic order parameters for asymptotically strong surface anchoring

The focus on the asymptotic regime where disc orientations are driven principally by surfaceanchoring provides us the opportunity to extract analytical estimates for the relevant nematic order parameters that we may compare against experimentally measured values. The nematic order parameters for the uniaxial nematic fluid are simply given by:

S (U ) ∼ P 2 (cos θ e ) ∆ (U ) = 0 M (U ) = 0 (45)
which only features θ e , reflecting the fact that uniaxial order is primarily determined by the optimal surface anchoring orientation of a single disc. For a nematic fluid with monoclinic and orthorhombic symmetries, disc-disc correlations are of crucial importance and we recall θ ne as the angle between n m and the colloidal director n c . We then find the following asymptotic expressions for the order parameters in the monoclinic and orthorhombic nematic states:

S (M O) ∼ P 2 (cos θ ne ) ∆ (M O) ∼ (sin θ ne ) 2 - π 3 16(cp) 2 M (M ) = sin(2θ ne ) M (O) = 0 (
S c ∼ 1 -O(β -1 ) ∆ c ∼ O(β -1 ) (47)
which demonstrates that the degree of orthorhombic biaxial order ∆ c in the colloid frame can be attributed to the azimuthal angular fluctuations of strength β -1 (Figs. 3d and4h). Numerical estimates for S c , ∆ c as well as for the colloidal director n c can be obtained from the colloidal ordering tensor diagonalised with respect to n c [cf. Eq. ( 15)]:

3 2 ω c ⊗ ω c - I 2 f (M O) c = S c 3 2 n c ⊗ n c - I 2 + ∆ c [(n c × n ⊥ ) ⊗ (n c × n ⊥ ) -n ⊥ ⊗ n ⊥ ] (48)
using asymptotic expressions for the tensor components ω c ⊗ ω c . Taking a typical concentration c = 3 just above the uniaxial-orthorhombic nematic transition in Fig. 2e which roughly corresponds to the disc concentration ρ = 0.31 µm -3 for the experimental system explored in Fig. 4e we find a small biaxiality parameter ∆ c ≈ 0.03 and strong uniaxial order (S c ≈ 0.98), independent of θ e . These values, which are also consistent with the experimental measurements, provide an a posteriori justification of the monoclinic order parameter in Eq. ( 18) which presupposes the disc and molecular orientational order within their respective director frames to be principally uniaxial.

We also find that the colloidal director angle θ ne (with respect to n m ) remains very close to the equilibrium anchoring angle θ e with |θ ne -θ e | < 1 • throughout the relevant temperature range.

Finally, the intrinsic skewness of the disc orientation distribution around its peak value can be gleaned from Eq. ( 13) retaining only the surface anchoring contributions:

f (M O) c (θ) ∼ exp - 4 n=2 ε n (θ -θ e ) n (49)
with the amplitudes fully specified in Eq. ( 11). In agreement with experiments, it is straightforward to establish a basic symmetry breaking f (θ e -θ) = f (θ -θ e ) in the angular fluctuations in the monoclinic nematic for which ε 3 = 0 (Fig. 4g,h and Extended Data Fig. 1l).

H. Impact of elasticity-mediated interactions

In order to estimate the importance of elastic interactions between the discs we begin by considering the far-field quadrupolar potential for spherical particles [76][77][78] :

U el (r) ∼ KR 0 8π r R 0 -5 w(ϑ r ) (50)
with

w(ϑ r ) = 9 -90 cos 2 ϑ r + 105 cos 4 ϑ r (51)
with r the centre-of-mass distance between the colloids, cos ϑ r = r • n m , R 0 = D/2 the typical colloid radius and K the (one-constant) Frank-Oseen elastic constant of the molecular nematic LC.

Even though the expression above does not capture the complicated orientation-dependence of the elastic interaction between freely rotating discs (Supplementary Movies 4 and 5), which requires identifying a multi-dimensional potential U el (r, ϑ r , ω c , ω c ), the simple form above encapsulates the basic quadrupolar nature of director-deformation-mediated interactions between colloidal inclusions in general. Strictly, in the isotropic and nematic fluid phase where cos ϑ r is randomly distributed, the above potential reduces to zero. Note that the potential becomes attractive at tilted pair configurations (with the optimum angle being ϑ r ≈ 50 • ) while co-axial and, to a lesser extent, co-planar configurations are strongly disfavoured.

In a fluid state certain pair configurations are energetically more favourable than others, so that fluctuations around the isotropic zero-point average are non-vanishing. Let us consider a Landau free energy defined as a Boltzmann average over the quadrupolar potential:

δF el ∼ ρ 2 drU el (r) e -U el (r) dre -U el (r) (52)
Linearising for weak amplitudes ||U el || 1 we find a simple expression proportional to the square of the quadrupolar potential:

δF el ∼ - ρ 8π dr[U el (r)] 2 (53)
In order to address the effect of the nemato-elastic medium on the correlations between disc orientations, which is the main driving force behind the uniaxial-monoclinic and uniaxial-orthorhombic nematic transitions, we must find a way to generalise the elastic potential Eq. ( 50) toward discotic particles. We consider an infinitely thin circular disc, and define the position vector for a random point on the surface as follows: r s = r + D 2 t(e 1 sin ξ + e 2 cos ξ) with 0 < t < 1 and 0 < ξ < 2π parameterising the disc surface (Supplementary Fig. S4 and Extended Data Fig. 10a) and e i representing unit vectors orthogonal to each other and to the disc normal ω c . We represent the elastic interaction between discs as a superposition of point-to-point segments interactions at distance ∆r s = |r s -r s | and write:

U el (r; ω c , ω c ) ∼ dtdξ dt dξ U el (∆r s (r; ω c , ω c )) (54)
The four-fold integral is analytically intractable, but the above expression may be simplified by invoking a far-field approximation and subsequently expanding the integrand up to quadratic order in D/|r| 1. The leading order shape-dependent contribution to the fluctuation term Eq. ( 53)

can then be computed numerically both for the uniaxial (U) state, where the relative azimuthal angle ∆ϕ is randomly distributed, and for the monoclinic and orthorhombic (MO) nematic phases where it is assumed zero. The numerical results are conveniently expressed as follows:

δF (U ) el ∼ 0 δF (M O) el ∼ u 0 c[a 0 + a 1 cos(2θ e ) + a 2 cos(4θ e )] (55)
The fitted parameters (a 0 , a 1 , a 2 ) ≈ (-0.59, 0.18, 0.41) ensure that the minimum of the bracketed term is fixed at -1 and yields zero when the discs are parallel to the nematic director (θ e = 0).

The amplitude u 0 > 0 is related to the elastic constant of the molecular LC as well as the disc dimension in some complicated way. It shall be considered as a small but adjustable parameter in the actual computations (sample Mathematica code in Appendix B). We infer that the elastic interactions tend to favour azimuthally aligned disc orientations over random ones. In addition, the degree at which monoclinic orientational order is stabilised depends sensitively on the anchoring angle θ e with the maximum effect reached at intermediate angles.

I. Onset of smectic and columnar order

We now propose a simple route toward probing the stability of the monoclinic nematic fluid with respect to density fluctuations signifying smectic or columnar order. Within density functional theory 79,80 , the inverse of the static structure factor of a perfectly aligned nematic fluid is given by:

S -1 (q) = 1 -ρp( φ)ĉ 2 (q) (56)
where ĉ2 (q) is the Fourier transform of the direct pair correlation function of the discs in the fluid phase and p( φ) a Lee-Parsons factor in terms of an effective disc packing fraction φ. For simplicity we assume all discs to be aligned along the colloidal director n c = (0, -sin θ ne , cos θ ne )

with θ ne coinciding with θ e . Adopting a simple van der Waals approximation, we express ĉ2 as a combination of a near field hard-core potential plus a far-field elastic contribution (cf. Eq. ( 50)):

ĉ2 (q) = -v ex (q) - r / ∈vex dre -iq•r U el (r) (57)
In order to perform the spatial integral, it is expedient to adopt a particle-based cylindrical frame r = (r ⊥ , α, r ) corresponding to the colloidal director that is tilted with respect to the molecular one (Supplementary Fig. S4). The wave vector describing a transverse modulation of the disc density reads:

q = (q T , 0, 0) (58)
The corresponding symmetry of the incipient phase could be either smectic or columnar. Here, we shall treat both possibilities on the same footing since the simple theory presented here is incapable of accounting for the subtle free energy differences between the two. The wavenumber q T defines the typical spacing between the columns or the smectic layers. Divergence of the static structure factor at a particular wave number signals an instability of the positionally uniform nematic fluid with respect to transverse density modulations:

S -1 (q T ) = 0 (59)
Let us now proceed towards calculating ĉ2 (q T ) explicitly. First, the volume exterior to the cylindrical hard core is parameterised as follows:

r / ∈vex dr ≈ 2 ∞ 0 dr ⊥ r ⊥ ∞ L dr 2π 0 dα (60)
in terms of a lower cut-off L reflecting the distance of closest approach between two discs. The Fourier transform of the cylindrical hard core of the discs reads:

vex (q T ) ≈ 2πLD 2 F(q T ) q 2 T ( 61
)
where F(q T ) = q T sin q T + cos q T -1 and q T implicitly normalised in terms of the disc diameter D. The inverse structure factor can then be written as:

S -1 (q) = 1 + ρp( φ) vex (q T ) + LD 2 Ξ el Ûel (q T ) (62)
with Ξ el ∼ KL 1 the elastic energy (renormalised in units of the thermal energy) quantifying the relative strength of (elastic) attractions with respect to steric repulsion. We stress that in view of the lack of long-ranged electrostatic interactions in our model Ξ el should be interpreted as a purely qualitative variational parameter. Consequently, the values for Ξ el displayed in the caption of Extended Data Fig. 10 do not relate to the experimentally measured average elastic constant K and disc thickness L. In order to evaluate the quadrupolar potential we first substitute cos ϑ r = -r • n c in Eq. ( 51). After some tedious algebra, the FT of the quadrupolar potential Eq. ( 50) can be rendered in closed form:

Ûel (q) ≈ π 2 12
x 2 w(θ e ) -3q 2 I 1 (q) -q3 I 2 (q) + q3 L -4 (q) + 9q 2 L -3 (q) + 12qL -2 (q) (63)

where q = q T x and I n (x) and L n (x) denote modified Bessel and Struve functions. Expanding the above expression for small disc aspect ratio x = L/D we find up to leading order:

Ûel (q T ) w(θ e ) ∼ π 3 q 2 T x 4 + O(x 5 ) (64)
The inverse structure factor then attains a particularly simple form, namely:

S -1 (q T ) ∼ 1 + φp( φ) 8 F(q T ) q 2 T + π 3 Ξ el x 4 w(θ e )q 2 T ( 65
)
The range of concentrations at which the structure factor diverges for a given temperaturecontrolled anchoring angle is indicated in the phase diagram in Fig. 2e and illustrated in Extended Data Fig. 10b. A sample Mathematica code can be found in Appendix C. The predictions from our simplified approach qualitatively agree with the experimental observation in the sense that stable smectic order only emerges across a very narrow temperature range. The present model also predicts the smectic layer spacing σ sm = 2π/q T at which S(q T ) diverges which is found to be consistent with experimental observations (Extended Data Fig. 10c).

II. DETAILS OF COMPUTER SIMULATIONS

A. Landau-de Gennes modelling of distortions in molecular order around discs

Numerical simulations are used to identify the local configurations and defect distributions of the molecular LC around the colloidal disc surface. The simulations are based on minimising the Landau-de Gennes free energy F LC of a continuum representation of a liquid crystal, given by:

F LC = d 3 rf bulk + d 2 rf surf ( 66
)
where the bulk energy density is integrated over the three-dimensional volume occupied by the LC whereas the surface contribution is integrated over the two-dimensional surface of the colloidal disc. In the Q-tensor representation, the bulk energy density f bulk is commonly expressed as [81][82][83] :

f bulk = A LdG 2 Q (m) ij Q (m) ji + B LdG 3 Q (m) ij Q (m) jk Q (m) ki + C LdG 4 (Q (m) ij Q (m) ji ) 2 + L 1 2 ∂Q (m) ij ∂x k 2 + L 2 2 ∂Q (m) ij ∂x j ∂Q (m) ik ∂x k + L 3 2 ∂Q (m) ij ∂x k ∂Q (m) ik ∂x j + L 6 2 Q (m) ij ∂Q (m) kl ∂x i ∂Q (m) kl ∂x j (67)
Here, the summation over all indices is assumed, and the molecular tensor order parameter is given by Q (m) = S m (3n m ⊗ n m -I)/2 in line with the definition Eq. ( 16). The first three terms in the bulk energy density are thermotropic parts describing the nematic-isotropic transition of the liquid crystal, with coefficients A LdG , B LdG , and C LdG being nematic material parameters. The equilibrium scalar order parameter is found by minimisation of the thermotropic energy:

S eq m = -B LdG + B LdG 2 -24A LdG C LdG 6C LdG ( 68
)
The last four terms with spatial derivatives in Eq. ( 67) represent the elastic energy due to intermolecular interactions and reflect the energy penalty associated with distorting the molecular director. The elasticity parameter L i , i = 1, 2, 3, 6 can be related to the Frank-Oseen elastic constants of the molecular nematic host like 5CB:

L 1 = 2 27(S eq m ) 2 (K 33 -K 11 + 3K 22 ) L 2 = 4 9(S eq m ) 2 (K 11 -K 24 ) L 3 = 4 9(S eq m ) 2 (K 24 -K 22 ) L 6 = 4 27(S eq m ) 3 (K 33 -K 11 ) (69)
The elastic energy is often simplified by applying a one-constant approximation: K = K 11 = K 22 = K 33 = K 24 , leaving only L 1 nonzero, though here we use the experimentally measured elastic constants, as detailed below.

The surface anchoring energy density f surf describes the effect of the colloidal surfaces on the local alignment of the liquid crystal molecules and the molecular director. For conically degenerate surface anchoring, f surf is independent of the azimuthal angle of the principal molecular axis with respect to the surface normal. In this case, the anchoring energy can be described by 84 :

f surf = W P ik Q(m) kl P lj - 3 2 S eq m cos 2 θ e P ij 2 (70)
with W the anchoring strength, P = v ⊗ v the surface projection tensor, v the surface normal director, and

Q(m) = Q (m) + 1 2 S eq m I.
To see the connection to the anchoring potential defined before, we can recast Eq. ( 70) by assuming there is no spatial variation of the scalar order parameter S m = S eq m and get:

f surf = 9 4 W (S eq m ) 2 cos 2 θ -cos 2 θ e 2 (71)
Therefore, compared with Eq. ( 5):

W = W 2 9(S eq m ) 2
(72)

The simulation adopts a finite-difference method with equidistant grid spacing 6.43 nm in all dimensions. The typical computational volume is a box with roughly 400 × 400 × 80 grid points in the colloidal frame, although larger or smaller boxes are used depending on the structures of interest. Soft boundary conditions (W = 0) are applied for all box boundaries and the liquid crystal has a uniform initial configuration (n m = (0, 0, 1)) prior to relaxation using Euler's method. More detailed information can be found in the code in Appendix D. The disc dimensions and parameters are chosen to match the experimental values. For all simulations the following parameters are used 82 : A LdG = -1.72 × 10 5 J/m 3 , B LdG = -2.12 × 10 6 J/m 3 , C LdG = 1.73 × 10 6 J/m 3 , K 11 = 6.4 pN, K 22 = K 24 = 3 pN, K 33 = 10 pN. These values of the Frank elastic constants correspond to those reported for the used 5CB host [START_REF] Liu | Biaxial ferromagnetic liquid crystal colloids[END_REF][START_REF] Mundoor | Electrostatically controlled surface boundary conditions in nematic liquid crystals and colloids[END_REF]84,85 . From Eq. ( 69) one finds: S eq m = 0.533, L 1 = 3.29 pN, L 2 = 5.32 pN, L 3 = 0, L 6 = 3.52 pN.

B. Numerical visualisation and simulation of polarising optical images

The analysis of the numerical results includes plotting isosurfaces of S m /S eq m , which shows the region of reduction in the scalar nematic order parameter. If not stated differently, the isosurface is plotted in black with S m /S eq m = 0.95. Colloidal particles with non-zero surface anchoring strength induce local deformations of the LC field, which is characterised by the contours of distortion.

Green regions (like the ones shown in Fig. 1j) indicate LC fields having a deviation of more than 2 degrees from the uniform background. Blue and orange contours (like the ones shown in Fig. 1g,h) are calculated similarly after projecting the LC field onto the plane parallel to the disc, and the two colours indicate the opposite direction of tilt, in a manner analogous to polarising optical micrographs (Fig. 1e). We also visualise the LC field by field streamlines, whose tangent indicate the local field direction. The distorted and uniform field lines around the discs are coloured according to orientations that are consistent with colours in polarising optical micrographs obtained with a phase retardation plate (Fig. 1e,f).

By means of the ParaView visualisation software 88 (open-source freeware obtained from Kitware, Inc.), we have presented various view perspectives and animations that show the symmetry and spatial extent of such n m -and S m -perturbations around individual discs within the colloidal LC in various phases (Fig. 1f-m, Extended Data Fig. 2 and Supplementary Videos 4 and 5). When modelling columnar colloidal self-assembly within the 5CB host, discs were placed at positions extracted from the confocal 3D imaging experiments while superimposing structures of molecular perturbations due to individual discs to obtain initial conditions. Upon further minimisation of the total free energy, we found overlapping of the coronae of the molecular director distortions due to individual discs within the neighbouring columns (Fig. 5h and Extended Data Fig. 9d-f), which shows that the monoclinic columnar colloidal structures in our work emerge from the competition between repulsive electrostatic and highly anisotropic elastic interactions both within the columns and in-between columns forming the oblique lattice. The presence of weak elastic distortions of n m around individual discs also modifies the overall colloidal interactions and influences collective behaviour in the nematic and smectic phases, as we discuss above.

To compare our numerical findings with experimental results for individual discs inside the molecular LC, images of polarising optical microscopy (POM) are generated by the Jones-matrix approach 86,87 . First, the local Jones matrix of each grid point is calculated based on the optical axis (local molecular director) orientation, defined by the orientation of molecules, and the ordinary/extraordinary phase retardation, defined by the optical anisotropy of the LC. Each pixel of the simulated image is then obtained by successive multiplication of Jones matrices across the direction of the light path. To approximate experimental images, the whole visible spectrum is represented by the combination of the red, green and blue colour wavelengths (640, 540, and 430 nm) with different light source intensities matching that of the experimental source. After separate calculation for each wavelength, the three single-coloured images are combined to give the final computer-simulated image. The Matlab code for POM simulation is given in Appendix D at the end of Supplementary Information. An optical anisotropy ∆n = 0.2 of 5CB is used in this modelling of POM micrographs, such as the ones shown in the insets of Fig. 1d,e.

III. PARAMETER SPECIFICATIONS FOR THE EXPERIMENTAL PHASE DIAGRAM

The phase diagram presented in Fig. 2a provides an overview of different phases and coexistence regions observed in experiment, with the coordinate axes representing temperature, disc number density and electrostatic charge of the colloids. Supplementary Tables S1 andS2 provide detailed information on the specific parameter ranges within which different phases and coexistence regions have been observed. The data in each Table corresponds to either one of the planar projections (denoting an effective particle charge of +20e and +80e, respectively) depicted in Fig. 2a. The relative error of measuring the disc number density based on 3D imaging and counting particles is less than ±5% , whereas the error of measuring temperature is ±0.1 o C. The Table S2. Temperature ranges (in ºC) of existence of different phases at Z*e=+20e

Table S2: Same lay-out as Supplementary Table S1, but for a disc charge of +20e

Appendix: Numerical codes A. Mathematica code for isotropic-uniaxial coexistence (Fig. 2e) a c t u a l = c e l l ( 1 , 3 ) ; % t h e i n d e x t o t a k e o u t t h e a c t u a l p a r t

Remove
f o r i = 1 : 3 % x y z a c t u a l { i } = (2-PBC ( i ) ) : ( n ( i ) -1+PBC ( i ) ) ; end i s L C r e a l = isLC ; isLC = z e r o s ( n ) ; isLC ( a c t u a l { 1 } , a c t u a l { 2 } , a c t u a l { 3 } ) = i s L C r e a l ;
% now t h i n k i t a s a l l PBC r = r e p m a t ( z e r o s ( n ) , 1 , 1 , 1 , 3 ) ; n = s i z e ( isLC ) ;

f o r i = 1 : 3 % l e t i == d i r e c t i o n o f l o o k i n g f o r i n d e x 2 = 1 : n ( pb ( i -1 ,3) ) f o r i n d e x 3 = 1 : n ( pb ( i + 1 , 3 ) ) nmax = n ( i ) + 4 ; % PBC : d u p l i c a t e 1 ˜4 t o t h e end i n d = { pb ( 1 : nmax , n ( i ) ) , i n d e x 2 , i n d e x 3 } ; % i n d e x 1 == d i r e c t i o n o f l o o k i n g i n d e x 1 = f i n d (1-isLC ( i n d { i } , i n d { pb ( i -1 ,3) } , i n d { pb ( i + 1 , 3 ) } ) , 1 ) ; % f i n d f i
r ( i n d { i } , i n d { pb ( i -1 ,3) } , i n d { pb ( i + 1 , 3 ) } , i ) = -4; end w h i l e i n d e x 1 <= nmax i n d = { pb ( ( i n d e x 1 + 1 ) : nmax , n ( i ) ) , i n d e x 2 , i n d e x 3 } ; i n d e x 1 = i n d e x 1 + f i n d ( isLC ( i n d { i } , i n d { pb ( i -1 ,3) } , i n d { pb ( i + 1 , 3 ) } ) , 1 ) ; % f i n d LC s u r f i n d = { pb ( ( i n d e x 1 + 1 ) : nmax , n ( i ) ) , i n d e x 2 , i n d e x 3 } ; t h = f i n d (1-isLC ( i n d { i } , i n d { pb ( i -1 ,3) } , i n d { pb ( i + 1 , 3 ) } ) , 1 ) ; % f i n d d i s t a n c e t o n e x t p a r t i c l e s u r f == t h i c k n e s s o f LC i f t h i n d = { pb ( [ i n d e x 1 i n d e x 1 + t h -1] , n ( i ) ) , i n d e x 2 ,
i n d e x 3 } ; % LC [ p o s i t i v e n e g a t i v e ] s u r f

r ( i n d { i } , i n d { pb ( i -1 ,3) } , i n d { pb ( i + 1 , 3 ) } , i ) =
[ min ( t h , 4 ) -min ( t h , 4 ) ] ;

i n d e x 1 = i n d e x 1 + t h ; % p a r t i c l e s u r f e l s e

i n d = { pb ( i n d e x 1 , n ( i ) ) , i n d e x 2 , i n d e x 3 } ; % LC p o s i t i v e s u r f o n l y r ( i n d { i } , i n d { pb ( i -1 ,3) } , i n d { pb ( i + 1 , 3 ) } , i

  Figure S1: a, Overview of the main vectors and angles parameterising the disc orientation. The excluded volume of two thin repulsive discs with L D at random orientations ω c and ω c is proportional to |ω c × ω c |. b, The surface anchoring potential U e causes the polar angle θ between the disc normal vector ω c and the molecular director n m to adopt a preferred value θ e . Notably, U e [given by Eq. (4)] is asymmetric about the equilibrium anchoring angle θ e .

Figure S2 :

 S2 Figure S2: Surface anchoring amplitudes in units of the thermal energy k B T [cf. Eq. (11)] versus equilibrium anchoring angle θ e (in radians) for discs with diameter D = 2 µm, polar surface anchoring coefficient W 2 = 10 -5 J/m 2 and average host elastic modulus K = 6 pN.

Figure S3 :

 S3 Figure S3: Schematic of the different symmetries of the hybrid nematic phase. a-c, the reference frame of a single disc at θ e = 90 • , 75 • , and 0 • , respectively. When 0 • < θ e < 90 • , the orthogonal disc frame {ω c , e 1 , e 2 } is defined by e 2 = ω c × n m and e 1 = e 2 × ω c . The green isosurface delineates the zone where the n m -field distortion from the far-field background is more than 2 • . d-f, the nematic phases at θ e = 90 • , 75 • , and 0 • , respectively. The monoclinicity of the hybrid system is demonstrated by a single mirror symmetry plane spanned by n c and n m when 0 • < θ ne < 90 • (M = 0).

  46) which additionally depend on the disc concentration c. Results are presented in the Fig. 4e, along with the corresponding experimental data. Defining uniaxial (S c ) and orthorhombic biaxial (∆ c ) order of the colloidal discs measured along the colloidal director n c we find the following scaling expressions:

Figure

  Figure S4: Schematic of the principal vectors and angles used for the estimation of the elastic interaction between a pair of thin discs. The surface of each disc can be parameterised within a particle-based frame {ω c , e 1 , e 2 }.

  * (4/Pi) * ldeff ; phin = cn * (4/Pi) * ldeff ; ffi = (1-(3/4) * phii)/(1-phii ) ˆ2; ffn = (1-(3/4) * phin)/(1-phin)ˆ2; dffi = (-5+3 phii)/(4 (-1+phii)ˆ3) ; dffn = (-5+3 phin)/(4 (-1+phin)ˆ3); pii = ci+ciˆ2 * ( ffi + dffi * phii ) ; mui = Log[ci]+2 * ci * ( ffi +(1/2) * phii * dffi ) ; pnn = cn+2 * cn * (1+ffnˆ-1 * dffn * phin); ( * equality of pressure and chemical potential * ) mun=Log[cn] + 2 * Log[cn * ffn] + Log[4/Pi] -1 + 4 * (1+(1/2) * phin * ffnˆ-1 * dffn);sv={ci,cn}/.FindRoot[{pii == pnn,mui == mun},{ci,cis},{cn,cns}]; Mathematica code for uniaxial-monoclinic/orthorhombic coexistence (Fig.2e) Remove[All] ( * starting values for coexistence concentrations at tilt angle tm = Pi/2 * ) cus = 2.65; cbs = 5.67; ( * relative strength of elasticity -mediated interactions * ) = Pi/2, tm > 0.1, tm = tmstep; phiu = (4/Pi) * ldeff * cu; phib = (4/Pi) * ldeff * cb; ffu = (1 -(3/4) * phiu)/(1 -phiu)ˆ2; ffb = (1 -(3/4) * phib)/(1 -phib)ˆ2; dffu = (-5 + 3 phiu)/(4 (-1 + phiu)ˆ3) ; dffb = (-5 + 3 phib)/(4 (-1 + phib)ˆ3) ; bb = (64/Piˆ3) * (cb * ffb) ˆ2; rhou = (4/Pi) * (2/Pi) * (Sin[tm]) * (1 + ArcTanh[Sin[tm]] * Cos[tm] * Cot[tm]); rhon = (4/Pi) * (1/( cb * ffb) ) ; pressu = ( ffu + dffu * phiu); pressb = (1 + ffb ˆ-1 * dffb * phib); muu = ffu + (1/2) * phiu * dffu; mub = 1 + (1/2) * phib * ffbˆ-1 * dffb; flucu = 0; flucb = cb * u0 * (-0.59 + 0.18 * Cos[2 * tm] + 0.41 * Cos[4 * tm]); ( * equal chemical potential and pressure * )sv = {cu,cb}/.FindRoot[{cu+cuˆ2 * rhou * pressu-flucu * cu==cb+cb * pressb-flucb * cb,Log[cu]+2 * cu * muu * rhou-Log[2 * Pi]-2 * flucu==Log[cb]+(1/2) * Log[bb * Sin[tm]ˆ2/(2 * Pi)]-(1 /2)+2 * mub-2 * flucb},{cu, cus},{cb,cbs}]; ( * concentration U and M/O phase: * ) cus = sv [[1]]; cbs = sv [[2]]; ffbb = (1-(3/4) * (4/Pi) * ldeff * cbs)/(1-(4/Pi) * ldeff * cbs)ˆ2; bbs = (64/Piˆ3) * sv [[2]]ˆ2 * ffbb ˆ2; ] Exit [] C. Mathematica code for the smectic stability analysis (Fig. 2e and Extended Data Fig. 10) = Pi/2, tm > 0.005, tm = tmstep; ww = 9 -90 * Cos[tm]ˆ2 + 105 * Cos[tm]ˆ4; dsf = D[(8 * (-1 + Cos[qq] + qq Sin[qq])/qqˆ2) + xiel * qqˆ2 * ww,qq]; qqs = qq/.FindRoot[dsf == 0,{qq,qqstart}]; qqstart = qqs; phis = (-1/(8 * (-1 + Cos[qq] + qq Sin[qq])/qqˆ2) + xiel * qqˆ2 * ww)/.{qq->qqs}; If [phis < 0,phis = 0]; phisparsons = (2 * (1+2 * phis-Sqrt[1 + phis]))/(3 + 4 * phis); e p m a t ( r e s h a p e ( [ 0 , s i n ( t h e t a ) , c o s ( t h e t a ) ] , 1 , 1 , 1 , 3 ) , nx , ny , nz , 1 ) ; % u n i f o r m z n = n + 0 . 0 0 5 * ( rand ( s i z e ( n ) ) -0.5) ; % p e r t u r b %% LC p a r a m e t e r s % t h e r m o t r o p i c % u n i t : J /mˆ3 == kg /ma s t i c % u n i t : N == J /m == kg * m/ s % d o v e r p r a t i o , number o f i d e a l p i t c h i n box ( zd i r e c t i o n ) % c o n i c a l l y d e g e n e r a t e s u r f a c e a n c h o r i n g Wc = 5E-4; % a n c h o r i n g s t r e n g t h % u n i t : J /mˆ2 t h e t a = 0 * p i / 1 8 0 ; % p o l a r a n g l e b e t w e e n s u r f a c e n o r m a l and p r e f e r e d a n g l e % u n i t : r a d i a n %% d e r i v e d c o n s t a n t s U = 0 ; r = r o l e ( PBC , isLC ) ; [ b u l k d , p s u r f d , n s u r f d , b u l k , s u r f a c e ] = r e g i o n s ( r , isLC ) ; [ n , v ] = u p d a t e n v ( n , v , r ) ; Seq = (-B+ s q r t (Bˆ2 -24 * A * C) ) / 6 / C ; % e q u i l i b r i u m s c a l a r o r d e r p a r a m e t e r m a x i t e = 1E5 ; % max number o f i t e r a t i o n s maxtime = 2 4 ; % u n i t : h o u r s %% i t e r a t i o n % i n i t i a l i z e e n e r g y l i s t = z e r o s ( 1 , m a x i t e ) ; % f r e e e n e r g y d e n s i t y , u n i t : J /mˆ3 Q ( : , : , : , 6 ) = -Q ( : , : , : , 1 )-Q ( : , : , : , 4 ) ; % c h a n g e t o 6 v a r i a b l e s % f i x e d s t e p s i z e f o r c o u n t 2 = 1 : m a x i t e % r e l a x a t i o n w i t h f i x e d s t e p s i z e f o r i = 1 : c h e c k i t e [ Qx , Qy , Qz , Qxx , Qyy , Qzz , Qxy , Qyz , Qzx ] = u p d a t e d e r i v a t i v e (Q, b u l k d , p s u r f d , n s u r f d , s p a c i n g ) ; fQ = dfdQ (Q, Qx , Qy , Qz , Qxx , Qyy , Qzz , Qxy , Qyz , Qzx , Ux , Uy , Uz , v , u , A , B , C , Wc , t h e t a , L1 , L2 , L3 , L4 , L6 , e0 , eamol , b u l k , s u r f a c e ) ; Q = Qs t e p s i z e * fQ . * ( u p d a t e s u r f + ( ˜u p d a t e s u r f * b u l k ) ) ; Q ( : , : , : , [ 1 4 6 ] ) = Q ( : , : , : , [ 1 4 6 ] )-mean (Q ( : , : , : , [ 1 4 6 ] ) , 4 ) ; % t r a c e l e s s n e s s end % c a l c u l a t e e n e r g y f = F r e e E n e r g y (Q, Qx , Qy , Qz , Ux , Uy , Uz , v , A , B , C , Wc , t h e t a , L1 , L2 , L3 , L4 , L6 , e0 , eamol , eavg , b u l k , s u r f a c e ) ; % u n i t : L1 * du e n e r g y l i s t ( c o u n t 2 ) = mean ( f , ' a l l ' ) * L1 / du ˆ2 ; % u n i t : J /mˆ3 % c h e c k t i m e i f t o c ( t 0 ) / 3 6 0 0 > maxtime % f p r i n t f ( ' l o o p c o u n t %d : maxtime r e a c h e d \n ' , c o u n t 2 ) ; break end % c o n v e r g e c o n d i t i o n i f max ( fQ , [ ] , ' a l l ' ) < 1E-9 % f p r i n t f ( ' l o o p c o u n t %d : e n e r g y c o n v e r g e \n ' , c o u n t 2 ) ; c t i o n p e r i o d i c b o u n d a r y = pb ( i n d e x , nmax ) p e r i o d i c b o u n d a r y = mod ( i n d e x -1 ,nmax ) + 1 ; end f u n c t i o n r = r o l e ( PBC , isLC ) % s t r u c t u r e % r o l e : +/-1 , 2 , 3 , 4 i s s u r f , 5 i s b u l k , 0 i s p a r t i c l e % 1 , 2 , 3 , ( > = ) 4 i s LC t h i c k n e s s b e t w e e n two w a l l s , % d e t e r m i n e s t h e o r d e r o f f i n i t e d i f f e r e n c e % +/-i s t h e d i r e c t i o n o f s u r f % make a s h e l l i f n o t PBC n = s i z e ( isLC ) +2 * (1 -PBC ) ;

  r s t p a r t i c l e i f i n d e x 1 >= 5 % f i r s t n e g a t i v e b o u n d a r y i n d = { i n d e x 1 -1 , i n d e x 2 , i n d e x 3 } ;

  isLC . * ( r ==0) ; % b u l k % t a k e a c t u a l p a r t r = r ( a c t u a l { 1 } , a c t u a l { 2 } , a c t u a l { 3 } , : ) ; end f u n c t i o n [ b u l k d , p s u r f d , n s u r f d , b u l k , s u r f a c e , bulkV , s u r f A ] = r e g i o n s ( r , isLC ) % s t r u c t u r e s p e c i f i c a t i o n u s e d i n f u n c t i o n s % f o r d and dd : l o g i c a l nx-ny-nz -3(-4) % s u r f d ( x , y , z , d i r e c t i o n , t h i c k n e s s ) b u l k d = ( r == 5 ) ; % u n i t l e s s [ nx , ny , nz , ˜] = s i z e ( r ) ; p s u r f d = f a l s e ( nx , ny , nz , 3 , 4 ) ; p s u r f d ( : , : , : , : , 1 ) = ( r == 1 ) ; p s u r f d ( : , : , : , : , 2 ) = ( r == 2 ) ; p s u r f d ( : , : , : , : , 3 ) = ( r == 3 ) ; p s u r f d ( : , : , : , : , 4 ) = ( r == 4 ) ; n s u r f d = f a l s e ( nx , ny , nz , 3 , 4 ) ; n s u r f d ( : , : , : , : , 1 ) = ( r == -1) ; n s u r f d ( : , : , : , : , 2 ) = ( r == -2) ; n s u r f d ( : , : , : , : , 3 ) = ( r == -3) ; n s u r f d ( : , : , : , : , 4 ) = ( r == -4) ; % f o r dFdQ : l o g i c a l nx-ny-nz b u l k = ( min ( abs ( r ) , [ ] , 4 ) == 5 ) ; % b u l k i n a l l d i r e c t i o n s % u n i t : du ˆ2 / L1 s u r f a c e = isLC & ( ˜b u l k ) ; % s u r f a c e w i t h t h i c k n e s s o f 1 % u n i t : du / L1 % f o r F r e e E n e r g y : d o u b l e nx-ny-nz bulkV = isLC ; % u n i t : du ˆ3 s u r f A = isLC . * s q r t ( sum ( r < 5 , 4 ) ) ; % u n i t : du ˆ2 end f u n c t i o n [ n , v ] = u p d a t e n v ( n0 , v0 , r ) % n o r m a l i z e and f i x v a t box w a l l s n = n o r m a l i z e ( n0 , 4 , ' norm ' ) ; % n o r m a l i z e n ( ( i s n a n ( n ) ) ) = 0 ; v = v0 ; s u r f a c e = ( min ( abs ( r ) , [ ] , 4 ) ˜= 0 ) & ( min ( r , [ ] , 4 ) ˜= 5 ) ; % u n i t : s p a c i n g / L1 v = v . * s u r f a c e ; % k e e p o n l y b o u n d a r y v = v -2 * ( sum ( v . * ( r <5) . * r , 4 ) < 0 ) . * v ; % c h e c k v o r i e n t a t i o n v = n o r m a l i z e ( v , 4 , ' norm ' ) ; v ( ( i s n a n ( v ) ) ) = 0 ; end f u n c t i o n S = s c a l a r O P ( n , r , b u l k d , p s u r f d , n s u r f d , Seq ) nb = sum ( ( r ==5) + ( abs ( r ) >1) , 4 ) ; % number o f n e i g h b o r s nb ( nb ==0) = 1 ; % ˜isLC a v g c o s 2 = ( sum ( n . * n ( [ end 1 : end -1 ] , : , : , : ) , 4 ) . ˆ2 . * ( b u l k d ( : , : , : , 1 ) | any ( n s u r f d ( : , : , : , 1 , [ 2 3 4 ] ) , 5 ) ) + . . . sum ( n . * n ( [ 2 : end 1 ] , : , : , : ) , 4 ) . ˆ2 . * ( b u l k d ( : , : , : , 1 ) | any ( p s u r f d ( : , : , : , 1 , [ 2 3 4 ] ) , 5 ) ) + . . . sum ( n . * n ( : , [ end 1 : end -1 ] , : , : ) , 4 ) . ˆ2 . * ( b u l k d ( : , : , : , 2 ) | any ( n s u r f d ( : , : , : , 2 , [ 2 3 4 ] ) , 5 ) ) + . . . sum ( n . * n ( : , [ 2 : end 1 ] , : , : ) , 4 ) . ˆ2 . * ( b u l k d ( : , : , : , 2 ) | any ( p s u r f d ( : , : , : , 2 , [ 2 3 4 ] ) , 5 ) ) + . . . sum ( n . * n ( : , : , [ end 1 : end -1 ] , : ) , 4 ) . ˆ2 . * ( b u l k d ( : , : , : , 3 ) | any ( n s u r f d ( : , : , : , 3 , [ 2 3 4 ] ) , 5 ) ) + . . . sum ( n . * n ( : , : , [ 2 : end 1 ] , : ) , 4 ) . ˆ2 . * ( b u l k d ( : , : , : , 3 ) | any ( p s u r f d ( : , : , : , 3 , [ 2 3 4 ] ) , 5 ) ) ) . / nb ; S = 0 . 5 * ( 3 * a v g c o s 2 -1) . * Seq ; end f u n c t i o n Q 5 = n2Q 5 ( n , S ) [ nx , ny , nz , ˜] = s i z e ( n ) ; Q 5 = z e r o s ( nx , ny , nz , 5 ) ; Q 5 ( : , : , : , 1 ) = S / 2 . * ( 3 * n ( : , : , : , 1 ) . ˆ2 -1 ) ; Q 5 ( : , : , : , 2 ) = S / 2 . * ( 3 * n ( : , : , : , 1 ) . * n ( : , : , : , 2 ) ) ; Q 5 ( : , : , : , 3 ) = S / 2 . * ( 3 * n ( : , : , : , 1 ) . * n ( : , : , : , 3 ) ) ; Q 5 ( : , : , : , 4 ) = S / 2 . * ( 3 * n ( : , : , : , 2 ) . ˆ2 -1 ) ; Q 5 ( : , : , : , 5 ) = S / 2 . * ( 3 * n ( : , : , : , 2 ) . * n ( : , : , : , 3 ) ) ; %Q ( : , : , : , 6 ) = S / 2 . * ( 3 * n ( : , : , : , 3 ) . ˆ2 -1 ) ; end f u n c t i o n Qx = dx (Q, b u l k d , p s u r f d , n s u r f d , s p a c i n g ) % d e r i v a t i v e ( r e g u l a r mesh ) ( p s u r f d ( : , : , : , 3 , 3 ) | p s u r f d ( : , : , : , 3 , 4 ) ) . * ( -1 . 5 * Q+2 * Q ( : , : , [ 2 : nz 1 ] , : ) -0.5 * Q ( : , : , [ 3 : nz 1 : 2 ] , : ) ) + . . . n s u r f d ( : , : , : , 3 , 2 ) . * ( -Q ( : , : , [ nz 1 : ( nz -1) ] , : ) +Q) + . . . ( n s u r f d ( : , : , : , 3 , 3 ) | n s u r f d ( : , : , : , 3 , 4 ) ) . * ( 0 . 5 * Q ( : , : , [ nz -1 nz 1 : ( nz -2) ] , : ) -2 * Q ( : , : , [ nz 1 : ( nz -1) ] , : ) + 1 . 5 * Q) ; Qz = Qz / s p a c i n g ; end f u n c t i o n Qxx = ddx (Q, b u l k d , p s u r f d , n s u r f d , s p a c i n g ) % s e c o n d d e r i v a t i v e ( r e g u l a r mesh ) [ nx , ˜, ˜, ˜] = s i z e (Q) ; Qxx = b u l k d ( : , : , : , 1 ) . * ( Q ( [ nx 1 : ( nx -1) ] , : , : , : ) -2 * Q+Q ( [ 2 : nx 1 ] , : , : , : ) ) + . . . p s u r f d ( : , : , : , 1 , 4 ) . * ( 2 * Q-5 * Q ( [ 2 : nx 1 ] , : , : , : ) +4 * Q ( [ 3 : nx 1 : 2 ] , : , : , : )-Q ( [ 4 : nx 1 : 3 ] , : , : , : ) ) + . . . n s u r f d ( : , : , : , 1 , 4 ) . * ( -Q ( [ ( nx -2) : nx 1 : ( nx -3) ] , : , : , : ) +4 * Q ( [ nx -1 nx 1 : ( nx -2) ] , : , : , : ) -5 * Q ( [ nx 1 : ( nx -1) ] , : , : , : ) +2 * Q) ; Qxx = Qxx + p s u r f d ( : , : , : , 1 , 3 ) . * Qxx ( [ 2 : nx 1 ] , : , : , : ) + n s u r f d ( : , : , : , 1 , 3 ) . * Qxx ( [ nx 1 : ( nx -1) ] , : , : , : ) ; % FDCoef ( 2 , [ 0 1 2 ] ) == FDCoef (2 ,[ -1 0 1 ] ) Qxx = Qxx / s p a c i n g ˆ2 ; end f u n c t i o n Qyy = ddy (Q, b u l k d , p s u r f d , n s u r f d , s p a c i n g ) [ ˜, ny , ˜, ˜] = s i z e (Q) ; Qyy = b u l k d ( : , : , : , 2 ) . * ( Q ( : , [ ny 1 : ( ny -1) ] , : , : ) -2 * Q+Q ( : r f d ( : , : , : , 2 , 4 ) . * ( -Q ( : , [ ( ny -2) : ny 1 : ( ny -3) ] , : , : ) +4 * Q ( : , [ ny-1 ny 1 : ( ny -2) ] , : , : ) -5 * Q ( : , [ ny 1 : ( ny -1) ] , : , : ) +2 * Q) ; . ˆ2 ) ) /2-L6 . * Q1 . * Qxx1 -2. * L6 . * Q2 . * Qxy1-L6 . * Qx2 . * Qy1-L6 . * Qx1 . * Qy2 -L6 . * Qy1 . * Qy4-L6 . * Q4 . * Qyy1 -2. * L6 . * Q5 . * Qyz1-L6 . * Qx3 . * Qz1-L6 . * Qy5 . * Qz1+L4 . * ( Qy3-Qz2 )-L6 . * Qx1 . * Qz3-L6 . * Qy1 . * Qz5-L6 . * Qz1 . * Qz6 -2. * L6 . * Q3 . * Qzx1 -(L2+L3 ) . * ( Qxx1+Qxy2+Qzx3 )-L6 . * Q6 . * Qzz1-L1 . * ( Qxx1+Qyy1+Qzz1 ) -(Ux . ˆ2 . * e0 . * eamol ) / 3 ; d f d Q b u l k ( : , : , : , 2 ) =A. * Q2+C . * Q1 . ˆ2 . * Q2+2 * C . * Q2 . ˆ3 + 2 * C . * Q2 . * Q3 . ˆ2 + C . * Q2 . * Q4 . ˆ2 + 2 * C . * Q2 . * Q5 . ˆ2 + B . * ( Q1 . * Q2+Q2 . * Q4+Q3 . * Q5 ) +C . * Q2 . * Q6 .ˆ2 -L6 . * Qx1 . * Qx2-L6 . * Q1 . * Qxx2-2 * L6 . * Q2 . * Qxy2 + 0 . 5 * L6 . * Qx1 . * Qy1-L6 . * Qx2 . * Qy2+L6 . * Qx3 . * Qy3 + 0 . 5 * L6 . * Qx4 . * Qy4-L6 . * Qy2 . * Qy4+L6 . * Qx5 . * Qy5 + 0 . 5 * L6 . * Qx6 . * Qy6-L6 . * Q4 . * Qyy2-2 * L6 . * Q5 . * Qyz2-L6 . * Qx3 . * Qz2-L6 . * Qy5 . * Qz2-L6 . * Qx2 . * Qz3+L4 . * ( -0 . 5 * Qx3 + 0 . 5 * Qy5 + 0 . 5 * Qz1 -0.5 * Qz4 )-L6 . * Qy2 . * Qz5-L6 . * Qz2 . * Qz6-2 * L6 . * Q3 . * Qzx2 + ( L2+L3 ) . * ( -0 . 5 * Qxx2 -0.5 * Qxy1 -0.5 * Qxy4 -0.5 * Qyy2 -0.5 * Qyz3 -0.5 * Qzx5 ) +L1 . * ( -Qxx2-Qyy2-Qzz2 )-L6 . * Q6 . * Qzz2 -1 / 3 . * Ux . * Uy . * e0 . * eamol ; d f d Q b u l k ( : , : , : , 3 ) =A. * Q3+C . * Q1 . ˆ2 . * Q3+2 * C . * Q2 . ˆ2 . * Q3+2 * C . * Q3 . ˆ3 + C . * Q3 . * Q4 . ˆ2 + 2 * C . * Q3 . * Q5 . ˆ2 + C . * Q3 . * Q6 . ˆ2 + B . * ( Q1 . * Q3+Q2 . * Q5+Q3 . * Q6 )-L6 . * Qx1 . * Qx3-L6 . * Q1 . * Qxx3-2 * L6 . * Q2 . * Qxy3-L6 . * Qx3 . * Qy2-L6 . * Qx2 . * Qy3-L6 . * Qy3 . * Qy4-L6 . * Q4 . * Qyy3-2 * L6 . * Q5 . * Qyz3 + 0 . 5 * L6 . * Qx1 . * Qz1+L6 . * Qx2 . * Qz2-L6 . * Qx3 . * Qz3-L6 . * Qy5 . * Qz3 + 0 . 5 * L6 . * Qx4 . * Qz4+ L4 . * ( 0 . 5 * Qx2 -0.5 * Qy1 + 0 . 5 * Qy6 -0.5 * Qz5 ) +L6 . * Qx5 . * Qz5-L6 . * Qy3 . * Qz5 + 0 . 5 * L6 . * Qx6 . * Qz6-L6 . * Qz3 . * Qz6-2 * L6 . * Q3 . * Qzx3+L1 . * ( -Qxx3-Qyy3-Qzz3 ) + ( L2+L3 ) . * ( -0 . 5 * Qxx3 -0.5 * Qxy5 -0.5 * Qyz2 -0.5 * Qzx1 -0.5 * Qzx6 -0.5 * Qzz3 )-L6 . * Q6 . * Qzz3 -1 / 3 . * Ux . * Uz . * e0 . * eamol ; d f d Q b u l k ( : , : , : , 4 ) =A. * Q4+B . * ( Q2 . ˆ2 + Q4 . ˆ2 + Q5 . ˆ2 ) +C . * Q4 . * ( 2 . * Q2 . ˆ2 + 2 . * Q3 . ˆ2 + Q4 . ˆ2 + 2 . * Q5 . ˆ2 + Q6 . ˆ2 )-L6 . * Qx1 . * Qx4+Q1 . * ( C . * Q1 . * Q4 -L6 . * Qxx4 ) -(L2+L3 ) . * ( Qxy2+Qyy4+Qyz5 ) +L4 . * ( -Qx5+Qz2 )-L1 . * ( Qxx4+ Qyy4+Qzz4 ) -(L6 . * ( 4 . * Q2 . * Qxy4-Qy1 . ˆ2 -2 . * Qy2 . ˆ2 -2 . * Qy3 . ˆ2 -2 . * Qy5 .ˆ2 -Qy6 . ˆ2 + 2 . * Qx4 . * ( Qy2+Qz3 ) + 2 . * Qy5 . * Qz4+Qy4 . * ( 2 . * Qx2+Qy4 + 2 . * Qz5 ) + 2 . * ( Q4 . * Qyy4 + 2 . * Q5 . * Qyz4+Qz4 . * ( Qx3+Qz6 ) + 2 . * Q3 . * Qzx4+Q6 . *Qzz4) ) ) /2 -(Uy . ˆ2 . * e0 . * eamol ) / 3 ; d f d Q b u l k ( : , : , : , 5 ) =A. * Q5+B . * ( Q2 . * Q3+Q4 . * Q5+Q5 . * Q6 ) +C . * Q5 . * ( 2 * Q2 . ˆ2 + 2 * Q3 . ˆ2 + Q4 . ˆ2 + 2 * Q5 . ˆ2 + Q6 . ˆ2 )-L6 . * Qx1 . * Qx5+Q1 . * ( C . * Q1 . * Q5-L6 . * Qxx5 ) +L4 . * ( 0 . 5 * Qx4 -0.5 * Qx6 -0.5 * Qy2 + 0 . 5 * Qz3 ) +L1 . * ( -Qxx5-Qyy5-Qzz5 ) + ( L2+L3 ) . * ( -0 . 5 * Qxy3 -0.5 * Qyy5 -0.5 * Qyz4 -0.5 * Qyz6 -0.5 * Qzx2 -0.5 * Qzz5 ) +L6 . * ( -2 * Q2 . * Qxy5-Qx5 . * Qy2-Qx2 . * Qy5-Qy4 . * Qy5-Q4 . * Qyy5-2 * Q5 . * Qyz5 + 0 . 5 * Qy1 . * Qz1+Qy2 . * Qz2-Qx5 . * Qz3+Qy3 . * Qz3 + 0 . 5 * Qy4 . * Qz4-Qx3 . * Qz5-Qy5 . * Qz5 + 0 . 5 * Qy6 . * Qz6-Qz5 . * Qz6-2 * Q3 . * Qzx5-Q6 . * Qzz5 ) -1/3. * Uy . * Uz . * e0 . * eamol ; d f d Q b u l k ( : , : , : , 6 ) =A. *

( 2 .

 2 * t h e t a ) ) / 2 ; d f d Q s u r f ( : , : , : , 2 ) = -0.25. * L4 . * Q3 . * u1-L1 . * Qx2 . * u1 -0.5 * L2 . * Qx2 . * u1 -0.5 * L3 . * Qx2 . * u1-L6 . * Q1 . * Qx2 . * u1 -0.5 * L3 . * Qy1 . * u1-L6 . * Q2 . * Qy2 . * u1 -0.5 * L2 . * Qy4 . * u1-L6 . * Q3 . * Qz2 . * u1 -0.5 * L2 . * Qz5 . * u1 + 0 . 2 5 . * L4 . * Q5 . * u2 -0.5 * L2 . * Qx1 . * u2-L6 . * Q2 . * Qx2 . * u2 -0.5 * L3 . * Qx4 . * u2-L1 . * Qy2 . * u2 -0.5 * L2 . * Qy2 . * u2 -0.5 * L3 . * Qy2 . * u2-L6 . * Q4 . * Qy2 . * u2-L6 . * Q5 . * Qz2 . * u2 -0.5 * L2 . * Qz3 . * u2 + 0 . 2 5 . * L4 . * Q1 . * u3 -0 . 2 5 . * L4 . * Q4 . * u3-L6 . * Q3 . * Qx2 . * u3 -0.5 * L3 . * Qx5 . * u3-L6 . * Q5 . * Qy2 . * u3 -0.5 * L3 . * Qy3 . * u3-L1 . * Qz2 . * u3-L6 . * Q6 . * Qz2 . * u3+v1 . * v2 . * ( -0 . 5 * Seq +2 * Q1 . * v1 . ˆ2 + 4 * Q2 . * v1 . * v2 +2 * Q4 . * v2 . ˆ2 + 4 * Q3 . * v1 . * v3 +4 * Q5 . * v2 . * v3 +2 * Q6 . * v3 . ˆ2 ) . * Wc -1 . 5 . * Seq . * v1 . * v2 . * Wc. * c o s ( 2 . * t h e t a ) ; d f d Q s u r f ( : , : , : , 3 ) = 0 . 2 5 . * L4 . * Q2 . * u1-L1 . * Qx3 . * u1 -0.5 * L2 . * Qx3 . * u1 -0.5 * L3 . * Qx3 . * u1-L6 . * Q1 . * Qx3 . * u1-L6 . * Q2 . * Qy3 . * u1 -0.5 * L2 . * Qy5 . * u1 -0.5 * L3 . * Qz1 . * u1-L6 . * Q3 . * Qz3 . * u1 -0.5 * L2 . * Qz6 . * u1 -0 . 2 5 . * L4 . * Q1 . * u2 + 0 . 2 5 . * L4 . * Q6 . * u2-L6 . * Q2 . * Qx3 . * u2 -0.5 * L3 . * Qx5 . * u2-L1 . * Qy3 . * u2-L6 . * Q4 . * Qy3 . * u2 -0.5 * L3 . * Qz2 . * u2-L6 . * Q5 . * Qz3 . * u2 -0 . 2 5 . * L4 . * Q5 . * u3 -0.5 * L2 . * Qx1 . * u3-L6 . * Q3 . * Qx3 . * u3 -0.5 * L3 . * Qx6 . * u3 -0.5 * L2 . * Qy2 . * u3-L6 . * Q5 . * Qy3 . * u3-L1 . * Qz3 . * u3 -0.5 * L2 . * Qz3 . * u3 -0.5 * L3 . * Qz3 . * u3-L6 . * Q6 . * Qz3 . * u3+v1 . * v3 . * ( -0 . 5 * Seq +2 * Q1 . * v1 . ˆ2 + 4 * Q2 . * v1 . * v2 +2 * Q4 . * v2 . ˆ2 + 4 * Q3 . * v1 . * v3 +4 * Q5 . * v2 . * v3 +2 * Q6 . * v3 . ˆ2 ) . * Wc-1 . 5 . * Seq . * v1 . * v3 . * Wc. * c o s ( 2 . * t h e t a ) ; d f d Q s u r f ( : , : , : , 4 ) = ( L4 . * ( -( Q5 . * u1 ) +Q2 . * u3 ) -2. * ( L3 . * Qy2 . * u1+L3 . * Qy4 . * u2+L2 . * ( Qx2+Qy4 ) . * u2+L2 . * Qz5 . * u2+L3 . * Qy5 . * u3+L1 . * ( Qx4 . * u1 +Qy4 . * u2+Qz4 . * u3 ) +L6 . * ( Q1 . * Qx4 . * u1+Q2 . * Qy4 . * u1+Q3 . * Qz4 . * u1+Q2 . * Qx4 . * u2+Q4 . * Qy4 . * u2+Q5 . * Qz4 . * u2+Q3 . * Qx4 . * u3+Q5 . * Qy4 . * u3+Q6 . * Qz4 . * u3 ) ) +v2 . ˆ2 . * ( -Seq + 4 . * ( Q1 . * v1 . ˆ2 + 2 . * Q2 . * v1 . * v2+Q4 . * v2 . ˆ2 + 2 . * Q3 . * v1 . * v3 + 2 . * Q5 . * v2 . * v3+Q6 . * v3 . ˆ2 ) ) . * Wc-3. * Seq . * v2 . ˆ2 . * Wc. * c o s ( 2 . * t h e t a ) ) / 2 ; d f d Q s u r f ( : , : , : , 5 ) = 0 . 2 5 . * L4 . * Q4 . * u1 -0 . 2 5 . * L4 . * Q6 . * u1-L1 . * Qx5 . * u1-L6 . * Q1 . * Qx5 . * u1 -0.5 * L3 . * Qy3 . * u1-L6 . * Q2 . * Qy5 . * u1 -0.5 * L3 . * Qz2 . * u1-L6 . * Q3 . * Qz5 . * u1 -0 . 2 5 . * L4 . * Q2 . * u2 -0.5 * L2 . * Qx3 . * u2-L6 . * Q2 . * Qx5 . * u2-L1 . * Qy5 . * u2 -0.5 * L2 . * Qy5 . * u2 -0.5 * L3 . * Qy5 . * u2-L6 . * Q4 . * Qy5 . * u2 -0.5 * L3 . * Qz4 . * u2-L6 . * Q5 . * Qz5 . * u2 -0.5 * L2 . * Qz6 . * u2 + 0 . 2 5 . * L4 . * Q3 . * u3 -0.5 * L2 . * Qx2 . * u3-L6 . * Q3 . * Qx5 . * u3 -0.5 * L2 . * Qy4 . * u3-L6 . * Q5 . * Qy5 . * u3 -0.5 * L3 . * Qy6 . * u3-L1 . * Qz5 . * u3 -0.5 * L2 . * Qz5 . * u3 -0.5 * L3 . * Qz5 . * u3-L6 . * Q6 . * Qz5 . * u3+v2 . * v3 . * ( -0 . 5 *

  i t : L1 / du ( e n e r g y / a r e a ) f s u r f = (Wc . * ( Seq -4 . * (Q1 . * v1 . ˆ2 + 2 . * Q2 . * v1 . * v2+Q4 . * v2 . ˆ2 + 2 . * Q3 . * v1 . * v3 + 2 . * Q5 . * v2 . * v3+Q6 . * v3 . ˆ2 ) + 3 . * Seq . * c o s ( 2 . * t h e t a ) ) . ˆ2 ) / 1 6 ; % u n i t : L1 / du ˆ2 f = b u l k . * ( f t h e r m o + f g r a d + f f i e l d ) + s u r f a c e . * f s u r f ; end f u n c t i o n [ n , S ] = Q2n (Q, isLC ) % c o n v e r t Q-t e n s o r t o d i r e c t o r and s c a l a r o r d e r p a r a m e t e r [ nx , ny , nz , ˜] = s i z e (Q) ; n = z e r o s ( nx , ny , nz , 3 ) ; S = z e r o s ( nx , ny , nz ) ; f o r x = 1 : nx f o r y = 1 : nyf o r z = 1 : nz i f isLC ( x , y , z ) Q t e n s o r = [Q( x , y , z , 1 ) Q( x , y , z , 2 ) Q( x , y , z , 3 ) ; . . . Q( x , y , z , 2 ) Q( x , y , z , 4 ) Q( x ,y , z , 5 ) ; . . . Q( x , y , z , 3 ) Q( x , y , z , 5 ) Q( x , y , z , 6 ) ] ; [V, D] = e i g ( Q t e n s o r ) ; [ s , i n d ] = max ( d i a g (Dne . * o n e s ( nx , ny , nz ) ; n o = no . * o n e s ( nx , ny , nz ) ; %% POM s i m u l a t i o n p h i = a t a n 2 ( n ( : , : , : , 2 ) , n ( : , : , : , 1 ) ) + a n g l e p s ; t h e t a = a c o s ( n ( : , : , : , 3 ) ) ; i n t e n s i t y = z e r o s ( ny , nx , l e n g t h ( RGBwavelength ) ) ; % POM f o r c o l o r = 1 : l e n g t h ( RGBwavelength ) E f i e l d = r e p m a t ( r e s h a p e ( [ 1 ; 0 ] , [ 1 , 1 , 2 ] ) , nx , ny , 1 ) ; % l i g h t b e f o r e e n t e r i n g s a m p l e do = 2 * p i / RGBwavelength ( c o l o r ) * s p a c i n g ( 3 ) * n o ; de = 2 * p i / RGBwavelength ( c o l o r ) * s p a c i n g ( 3 ) * n e . * n o . / ( n o . ˆ2 . * s i n ( t h e t a ) . ˆ2 + n e . ˆ2 . * c o s ( t h e t a ) . ˆ2 ) . ˆ0 . 5 ; % p r o p o g a t i o n o f l i g h t f o r z = 1 : nz Ex=exp ( 1 i . * do ( : , : , z ) ) . * s i n ( p h i ( : , : , z ) ) . * ( -( E f i e l d ( : , : , 2 ) . * c o s ( p h i ( : , : , z ) ) ) + E f i e l d ( : , : , 1 ) . * s i n ( p h i ( : , : , z ) ) ) +exp ( 1 i . * de ( : , : , z ) ) . * c o s ( p h i ( : , : , z ) ) . * ( E f i e l d ( : , : , 1 ) . * c o s ( p h i ( : , : , z ) ) + E f i e l d ( : , : , 2 ) . * s i n ( p h i ( : , : , z ) ) ) ; Ey=exp ( 1 i . * do ( : , : , z ) ) . * c o s ( p h i ( : , : , z ) ) . * ( E f i e l d ( : , : , 2 ) . * c o s ( p h i ( : , : , z ) )-E f i e l d ( : , : , 1 ) . * s i n ( p h i ( : , : , z ) ) ) +exp ( 1 i . * de ( : , : , z ) ) . * s i n ( p h i ( : , : , z ) ) . * ( E f i e l d ( : , : , 1 ) . * c o s ( p h i ( : , : , z ) ) + E f i e l d ( : , : , 2 ) . * s i n ( p h i ( : , : , z ) ) ) ; E f i e l d = c a t ( 3 , Ex , Ey ) ; end i f w a v e p l a t e t h i c k = 7 * 5 . 8 8 8 9 e -05; de wp = 2 * p i / RGBwavelength ( c o l o r ) * t h i c k * 1 . 5 5 3 3 8 ; do wp = 2 * p i / RGBwavelength ( c o l o r ) * t h i c k * 1 . 5 4 4 2 5 ; Ex=exp ( 1 i . * do wp ) . * ( 0 . 5 . * E f i e l d ( : , : , 1 ) -0 . 5 . * E f i e l d ( : , : , 2 ) ) + exp ( 1 i . * de wp ) . * ( 0 . 5 . * E f i e l d ( : , : , 1 ) + 0 . 5 . * E f i e l d ( : , : , 2 ) ) ; Ey=exp ( 1 i . * do wp ) . * ( -0 . 5 . * E f i e l d ( : , : , 1 ) + 0 . 5 . * E f i e l d ( : , : , 2 ) ) + exp ( 1 i . * de wp ) . * ( 0 . 5 . * E f i e l d ( : , : , 1 ) + 0 . 5 . * E f i e l d ( : , : , 2 ) ) ; E f i e l d = c a t ( 3 , Ex , Ey ) ; end % a n a l y z e r i n t e n s i t y ( : , : , c o l o r ) = ( abs ( E f i e l d ( : , : , 1 ) * c o s ( a n g l e a p ) + E f i e l d ( : , : , 2 ) * s i n ( a n g l e a p ) ) . ˆ2 ) . ' ; % y by x i n t e n s i t y ( : , : , c o l o r ) = i m g a u s s f i l t ( i n t e n s i t y ( : , : , c o l o r ) , 9 * r a t e ( 1 ) ) * R G B s e n s i t i v i t y ( c o l o r ) ; end % p l o t i m a g e s c ( [ 1 , nx * dx ] , [ 1 , ny * dy ] , i n t e n s i t y ) a x i s e q u a l t i g h t x l a b e l ( ' x ' ) y l a b e l ( ' y ' ) s e t ( gca , ' YDir ' , ' n o r m a l ' )

Table S1 .

 S1 Temperature ranges (in ºC) of existence of different phases at Z*e=+80e

	Number density ρ (mm -3 )	0.14	0.16	0.18	0.20	0.22	0.24	0.26	0.28	0.30-0.42	0.42-0.48	0.48
	Volume fraction (%)	0.42	0.48	0.54	0.60	0.66	0.72	0.78	0.84	0.90-1.27	1.27-1.45	1.45
	Uniaxial nematic	24.0< T<34.7	24.0<T<31.5 & 33.4<T<34.7	24.0<T<30.5 & 33.5<T<34.7	24.0<T<30.2 & 33.5<T<34.7	33.3<T<34.6	33.5<T<34.6	33.5<T<34.6	33.5<T<34.6	33.5<T<34.6	33.5<T<34.5	33.5<T<34.5
	Orthorhombic nematic									24<T<30	24<T<30	24<T<30
	Monoclinic nematic					32.0<T<33.3	30.5<T<33.5	30.2<T<33.5	30.0<T<31.0 & 32.4<T<33.5	30.0<T<31.0 & 32.5<T<33.5	30.0<T<31.0 & 32.5<T<33.5	30.0<T<31.0 & 32.5<T<33.5
	Monoclinic smectic								31.0<T<32.4	31.0<T<32.5	31.0<T<32.5	31.0<T<32.5
	Isotropic	T>34.7	T>34.7	T>34.7	T>34.7	T>34.6	T>34.6	T>34.6	T>34.6	T>34.6
	Colloidal nematic											T>34.5
	Uniaxial-monoclinic nematic Co-ex		31.5<T<33.4	30.5<T<33.5	30.2<T<33.5	30.0<T<33.3	30.0<T<30.5	30.0<T<30.2
	Uniaxial-
	orthorhombic					24.0<T<30.0	24.0<T<30.0	24.0<T<30.0	24.0<T<30.0
	nematic Co-ex
	Isotropic-colloidal nematic Co-ex										T>34.5

Table S1 :

 S1 Overview of the temperature ranges for the observed hybrid colloidal-molecular phases and their co-existences corresponding to an effective colloid charge of +80e. Number densities and volume fractions are specified in the top rows.

	Number density ρ (mm -3 )	<0.01	0.01-0.14	0.16	0.18	0.20	0.22	0.24	0.26	0.28	0.30	0.32	0.34	0.36	0.38	0.40	0.42	0.44	0.46	0.48
	Volume fraction (%)	<0.03	0.03-0.42	0.48	0.54	0.60	0.66	0.72	0.78	0.84	0.90	0.97	1.03	1.09	1.15	1.21	1.27	1.33	1.39	1.45
	Uniaxial nematic with individual discs	24.0<T< 34.7	34.0<T< 34.7	34.0<T< 34.7	34.0<T< 34.7	34.0<T< 34.7
	Uniaxial nematic with colloidal columns of discs		24.0<T< 34.0	32.2<T< 34.0	33.2<T< 34.0	33.8<T< 34.0
	Orthorhombic nematic with individual discs											34.0<T< 34.6	34.0<T< 34.6	34.0<T< 34.6	34.0<T< 34.5	34.0<T< 34.5	34.0<T< 34.5	34.0<T< 34.5	34.0<T< 34.5
	Orthorhombic					24.0<T<	24.0<T<	24.0<T<	24.0<T<	24.0<T<	24.0<T<	31.5<T<	32.0<T<	32.4<T<	32.6<T<	32.8<T<	33.0<T<	33.2<T<	33.4<T<
	columnar nematic					31.5	33.4	33.8	34.0	34.0	34.0	34.0	34.0	34.0	34.0	34.0	34.0	34.0	34.0
	Columnar											24.0<T< 31.5	24.0<T< 32.0	24.0<T< 32.4	24.0<T< 32.6	24.0<T< 32.8	24.0<T< 33.0	24.0<T< 33.2	24.0<T< 33.4	24.0<T< 34.0
	Isotropic	T>34.7 T>34.7 T>34.7 T>34.7 T>34.7 T>34.6 T>34.6 T>34.6 T>34.6 T>34.6 T>34.6 T>34.6 T>34.6 T>34.5 T>34.5 T>34.5
	Colloidal nematic																			T>34.5
	Uniaxial -orthorhombic columnar Co-ex			24.0<T< 32.2	24.0<T< 33.2	24.0<T< 33.8
	Uniaxial-orthorhombic nematic Co-ex						34.0<T< 34.6	34.0<T< 34.6	34.0<T< 34.6	34.0<T< 34.6	34.0<T< 34.6
	Isotropic-colloidal nematic Co-ex																	T>34.5 T>34.5

volume fraction data are obtained from the measured number density while taking the average geometric dimensions of the used discs.

  Qx1 . ˆ2 + L6 . * Q1 . * Qx1 . ˆ2 + L3 . * Qx2 . ˆ2 + 2 . * L6 . * Q1 . * Qx2 . ˆ2 + L3 . * Qx3 . ˆ2 + 2 . * L6 . * Q1 . * Qx3 . ˆ2 + L6 . * Q1 . * Qx4 . ˆ2 + 2 . * L6 . * Q1 . * Qx5 . ˆ2 + L6 . * Q1 . * Qx6 . ˆ2 + 2 . * L6 . * Q2 . * Qx1 . * Qy1 + 2 . * L3 . * Qx2 . * Qy1+L6 . * Q4 . * Qy1 . ˆ2 + 4 . * L6 . * Q2 . * Qx2 . * Qy2 + 2 . * L3 . * Qx4 . * Qy2+L3 . * Qy2 . ˆ2 + 2 . * L6 . * Q4 . * Qy2 .ˆ2 -L4 . * ( Q2 . * Qx3 + (Q4-Q6 ) . * Qx5+Q3. * ( -Qx2+Qy1 ) +Q5. * ( -Qx4+Qx6+ Qy2 ) ) +L4 . * Q1 . * Qy3-L4 . * Q6 . * Qy3 + 4 . * L6 . * Q2 . * Qx3 . * Qy3 + 2 . * L3 . * Qx5 . * Qy3 + 2 . * L6 . * Q4 . * Qy3 . ˆ2 + 2 . * L6 . * Q2 . * Qx4 . * Qy4+L3 . * Qy4 . ˆ2 + L6 . * Q4 . * Qy4 . ˆ2 + L4 . * Q2 . * Qy5 + 4 . * L6 . * Q2 . * Qx5 . * Qy5+L3 . * Qy5 . ˆ2 + 2 . * L6 . * Q4 . * Qy5 . ˆ2 + L4 . * Q3 . * Qy6 + 2 . * L6 . * Q2 . * Qx6 . * Qy6+L6 . * Q4 . * Qy6 . ˆ2 + L4 . * Q2 . * Qz1 + 2 . * L6 . * Q3 . * Qx1 . * Qz1 + 2 . * L3 . * Qx3 . * Qz1 + 2 . * L6 . * Q5 . * Qy1 . * Qz1+L6 . * Q6 . * Qz1 .ˆ2 -L4 . * Q1 . * Qz2+L4 . * Q4 . * Qz2 + 4 . * L6 . * Q3 . * Qx2 . * Qz2 + 2 . * L3 . * Qx5 . * Qz2 + 4 . * L6 . * Q5 . * Qy2 . * Qz2 + 2 . * L3 . * Qy3 . * Qz2 + 2 . * L6 . * Q6 . * Qz2 . ˆ2 + L4 . * Q5 . * Qz3 + 4 . * L6 . * Q3 . * Qx3 . * Qz3 + 2 . * L3 . * Qx6 . * Qz3 + 4 . * L6 . * Q5 . * Qy3 . * Qz3+L3 . * Qz3 . ˆ2 + 2 . * L6 . * Q6 . * Qz3 .ˆ2 -L4 . * Q2 . * Qz4 + 2 . * L6 . * Q3 . * Qx4 . * Qz4 + 2 . * L6 . * Q5 . * Qy4 . * Qz4 + 2 . * L3 . * Qy5 . * Qz4+L6 . * Q6 . * Qz4 .ˆ2 -L4 . * Q3 . * Qz5 + 4 . * L6 . * Q3 . * Qx5 . * Qz5 + 4 . * L6 . * Q5 . * Qy5 . * Qz5 + 2 . * L3 . * Qy6 . * Qz5+L3 . * Qz5 . ˆ2 + 2 . * L6 . * Q6 . * Qz5 . ˆ2 + 2 . * L6 . * ( Q3 . * Qx6+Q5 . * Qy6 ) . *

	Seq +2 * Q1 . * v1 . ˆ2 + 4 * Q2 f g r a d = ( L3 . * Qz6 + ( L3+L6 . * Q6 ) . * Qz6 . ˆ2 + L1 . * ( Qx1 . ˆ2 + 2 . * Qx2 . ˆ2 + 2 . * Qx3 . ˆ2 + Qx4 . ˆ2 + 2 . * Qx5 . ˆ2 + Qx6 . ˆ2 + Qy1 . ˆ2 + 2 . * Qy2 . ˆ2 + 2 . * Qy3 . ˆ2 + Qy4 . ˆ2 + 2 . * Qy5 . ˆ2 + Qy6 . ˆ2 + Qz1 . ˆ2 + 2 . * Qz2 . ˆ2 + 2 . * Qz3 . ˆ2 + Qz4 . ˆ2 + 2 . * Qz5 . ˆ2 + Qz6 . ˆ2 ) +L2 . * ( ( Qx1+Qy2+Qz3 ) . ˆ2 + ( Qx2+Qy4+Qz5 ) . ˆ2 + ( Qx3+Qy5+Qz6 ) . ˆ2 ) ) / 2 ; . + 3 . * Q4 . * Q5 . ˆ2 + 3 . * ( Q3 . ˆ2 + Q5 . ˆ2 ) . * Q6+Q6 . ˆ3 ) ) / 1 2 ; f f i e l d = ( e0 .

* v1 . * v2 +2 * Q4 . * v2 . ˆ2 + 4 * Q3 . * v1 . * v3 +4 * Q5 . * v2 . * v3 +2 * Q6 . * v3 . ˆ2 ) . * * ( -2 . * ( Q1 . * Ux . ˆ2 + 2 . * Q2 . * Ux . * Uy+Q4 . * Uy . ˆ2 + 2 . * Q3 . * Ux . * Uz + 2 . * Q5 . * Uy . * Uz+Q6 . * Uz . ˆ2 ) . * eamol -3 . * (Ux . ˆ2 + Uy . ˆ2 + Uz . ˆ2 ) . *

Supplementary Information