### Simulation method

#### Event-chain Monte-Carlo

To simulate hard disks in a suspension of hard needles the event-chain MC algorithm for many-body interactions is used^{53, 54, 62}. This is a rejection-free MC method that performs very well for dense systems^{57}, which makes it an excellent choice for needle systems in the nematic phase, for which it has been adapted in Ref.^{54}. The event-chain MC is a Markov-chain MC scheme, which fulfills maximal (rejection-free) global balance rather than detailed balance. Global balance is achieved by introducing lifting moves^{62}. For hard spheres or needles a lifting move is the transfer of MC displacement from one particle to another particle. This means in a MC move only one particle at a time is moved along a line until it contacts another object. For the application to hard needles, one of the two endpoints are moved and, upon collision with another needle, the remaining MC move distance is transferred to one of the endpoints of this needle. To which of the endpoints it is transferred depends on the collision point along the needle^{54}. In the presence of additional disks, MC displacement is also transferred to disks if a needle collides with the disk and vice versa. Collision detection is the computational bottleneck, and we use a sophisticated neighbor list system to achieve high simulation speeds also in the nematic phase (see Supplementary Information).

#### Hard needle liquid crystal

We model the molecules of a lyotropic LC as hard needles, which consist of two endpoints connected by an infinitely thin, hard line (see Fig. 1a). For efficient sampling, the distance of the two endpoints, i.e. the length *l* of the hard needle can fluctuate around its characteristic length (l_0) in order to allow for independent motion of both endpoints in the MC simulation; the needle length *l* is restricted by a hard potential (V_text {n}(l)) with (V_text {n}(l) = 0) for (l/l_0 in [0.9, 1.1]) and infinite (V_text {n}(l)) else, such that (l_0) is the effective needle length. Restricting the needle length is also necessary for our neighbor lists, which accelerate the simulation.

The interaction potential (V_text {nn}) between two hard needles is either infinite when they overlap or zero else. This results in a lyotropic transition from isotropic to nematic upon increasing the needle density (rho _text {n}). In a system with (N_text {n}) needles the density (rho _text {n} = N_text {n}/A_text {free}) is defined using the available free area (A_text {free}) (reduced by the area occupied by disks). The simulation box is a square with edge length *L* and periodic boundary conditions. In a typical system of size (A_text {free} = 900 l_0^2) we simulate (N_text {n}sim 9000) needles at a needle density (rho _text {n}’=10).

Liquid crystalline ordering is described by the standard tensorial order parameter

$$begin{aligned} Q_{alpha beta } = frac{1}{N_text {n}} sum _{i=1}^{N_text {n}} left( vec {nu }^i_alpha vec {nu }^i_beta – frac{delta _{alpha beta }}{2} right) , , end{aligned}$$

which is invariant under needle inversion; (vec {nu }^i) is the orientation of needle *i*.

The scalar order parameter ({S in [0, 1]}) is calculated by diagonalizing the matrix, where ({lambda =2S}) is the biggest eigenvalue. The corresponding eigenvector is the director (vec {n}). The scalar order parameter *S* measures the degree of alignment in the system and is zero in a perfectly isotropic phase and unity in a perfectly aligned nematic phase. Figure 5 shows the order parameter *S* as a function of density for the two-dimensional hard needle system^{54}; the system quasi-orders above ({rho _text {n}’ approx 6})^{55}. In our two-dimensional system, this transition is a Kosterlitz-Thouless transition into a quasi-ordered nematic state^{56}.

### Elastic energy, weak anchoring and elastic anisotropy

For a nematic LC the state of perfect alignment is the energetically preferred configuration. But this state is almost never reached because of thermal fluctuations, boundary conditions, or topological defects. Often effective free energies from field theories are used to describe deviations from homogeneous alignment. For a director field (vec {n}(vec {r})) in three dimensions, the energy of small distortions can be captured with the Frank-Oseen free energy (F_text {FO})^{63} with three elastic constants (K_i), which determine the energy cost of the three different distortions splay ((K_1)), twist ((K_2)) and bend ((K_3)). In a two-dimensional system, twist distortions are not possible and, therefore, ({K_2=0}) resulting in^{60}

$$begin{aligned} F^text {2D}_text {FO}&= int {mathrm{d}} A left[ frac{K_1}{2} (nabla vec {n})^2 + frac{K_3}{2} (nabla times vec {n})^2 right] approx int {mathrm{d}} A left[ frac{K_3}{2} (partial _x Phi )^2 +frac{K_1}{2} (partial _y Phi )^2 right] end{aligned}$$

(5)

with the needle orientation field (Phi (vec {r})) related to the director via (vec {n}=(cos Phi , sin Phi )). We assume (Phi =0) corresponding to (vec {n}||vec {e}_x) asymptotically at the system boundary. The last approximation is the leading order in an expansion in deviations from (Phi =0) and neglects non-linear terms in (Phi )^{61}.

The elastic constants (K_1) and (K_3) can be calculated by adapting a method described by Straley^{59} to two dimensions. Strictly speaking the pronounced logarithmic orientational fluctuations in two dimensions will give rise to a renormalization at large scales^{64}, which we ignore in our finite size system. We obtain ({K_1 / rho _text {n}’^2 k_text {B}T = 0.086}) and ({K_3 / rho _text {n}’^2 k_text {B}T = 0.63}) for the nematic phase with (rho _text {n}’=10). This result is based on the angular distribution (p_A(theta ) propto exp (A cos ^2theta )) of needle orientations (vec {nu }=(cos theta , sin theta )) with (A=10) matching our MC simulations in the nematic phase with (rho _text {n}’=10) (see Fig. 5). In the one-constant approximation, we assume (K=frac{1}{2}(K_1+ K_3)) resulting in ({K / rho _text {n}’^2 k_text {B}T = 0.358}); for this value of *K* the quadrupolar interaction (2) is shown in 4.

In order to describe colloidal disks in a nematic needle phase faithfully, we need too take into account elastic anisotropy (K_1 ne K_3) and a finite parallel anchoring with an anchoring strength *W*. In order to quantify the elastic anisotropy (K_1/K_3) and the relative anchoring strength (w=Wsigma /2K_3), we consider the needle orientation field (Phi (vec {r})) around a single disk, which is suspended in a nematic phase with (Phi =0) asymptotically at the system boundary and quantitatively compare MC simulation data with minimization of the free energy (3) containing both the anisotropic 2D Frank-Oseen elastic energy and the anchoring potential. Free energy minimization using the linearized approximation in the Frank-Oseen part (5) results in the partial differential equation

$$begin{aligned} left( partial _{tilde{x}}^2 + frac{K_1}{K_3} partial _{tilde{y}}^2 right) Phi&=0~~(tilde{r}>1),~~~~ partial _{tilde{x}} Phi + frac{K_1}{K_3} partial _{tilde{y}} Phi = frac{w}{2} sin left( 2(Phi -Phi _0)right) ~~(tilde{r}=1). end{aligned}$$

(6)

in dimensionless units (tilde{vec {r}} equiv 2vec {r}/sigma ). Analytical solutions are only available in the one-constant approximation (K_1 = K_3), where^{65}

$$begin{aligned} Phi (r, vartheta )&= – arctan left[ frac{(sigma /2r)^2 p(w) sin (2vartheta )}{1 – (sigma /2r)^2 p(w) cos (2vartheta )} right] qquad text {with} qquad p(w) = frac{2}{w}left( sqrt{1 + frac{w^2}{4}} – 1 right) end{aligned}$$

(7)

(both (Phi ) and (vartheta ) are defined with respect to the asymptotic director orientation, which determines the *x*-axis). For the full problem (6) we have to resort to numerical solutions by finite element methods (using MATHEMATICA).

In Fig. 6 we plot MC simulation data for the director orientation field (Phi (r,vartheta )) as a function of (vartheta ) (i.e., along circles) for different rescaled (r/sigma ) and for different system sizes such that (L/sigma =6) remains fixed. Rescaling of lengths with (sigma ) results in good data collapse for different (sigma ). The characteristic shoulder of (Phi (r,vartheta )) as a function of (vartheta ) can only be explained by an elastic anisotropy (K_1/K_3approx 0.1). Fits with Eq. (7) from elastically isotropic one-constant approximation ((K_1=K_3)) remain unsatisfactory. We perform fits for specific (r/sigma ) and overall fits within the whole range (r/sigma =1.5-3.0) with the numerical solution of the PDE (6) in the same geometry as the MC simulation (square box with (L/sigma =6) and director oriented along the diagonal as fixed by Dirichlet boundary conditions) using (K_1/K_3=0.1) and *w* as fit parameter. The fit results for *w* are weakly *r*-dependent (dashed blue lines) and the result of the overall fit (solid blue line) is (wapprox 0.5), see Eq. (4).

### Quadrupolar elastic interaction

If two colloidal disks are embedded, boundary conditions on the disk surface induce deformations in the director field and, thus, induce an elastic interaction given by the energy (5) of the deformation. There are some results on these elastic interactions in the one-constant approximation (({K = K_1 = K_3})). For two disks with perfect planar boundary conditions there is an approximate analytical result which is the quadrupolar long-range interaction from Eq. (2)^{43, 58}. Since this is a field theoretical result, the needles are assumed to be infinitely small relative to the colloids. There are no analytical results for the full two-constant free energy (5).

Therefore, we performed numerical simulations based on the finite element method (using MATHEMATICA), i.e., solving Eq. (6) for a system with two disks with distance *r* and with an angle (vartheta ) with respect to the director axis and using the same anisotropy (K_1/K_3=0.1) and weak anchoring strength (w=0.5) that we determined numerically. We place both disks symmetrically in a system in a square box with (L/sigma =10) and with a director axis (Phi =0) oriented parallel to one edge of the simulation box as fixed by Dirichlet conditions at the system boundaries. The total energy of the system is the sum of elastic and anchoring energy (see Eq. (3)). The interaction energy (U_text {quad}(r, vartheta )) is obtained as the difference in energy between a system containing two disks and the non-interacting system as obtained by the sum of the energies of two systems containing only one of the disks each. This results in Fig. 4b.

### Density-dependent depletion interaction

We use the results of Biben et al.^{66} and generalize them to anisotropic depletants with a rotational degree of freedom (varphi ) to get a density-dependent depletion interaction for disks in a suspension of hard needles. We consider a system of hard disks with positions ({ vec {X}_I }) and (N_text {n}) hard needles with positions ({ vec {x}_i }) and orientations ({ varphi _i }). Upper case indices refer to disks, lower case indices to needles. The energy of the system is given by

$$begin{aligned} H = sum _{I<J} V_text {dd}(vec {X}_I – vec {X}_J) + sum _{i<j} V_text {nn}(vec {x}_i – vec {x}_j, varphi _i, varphi _j) + sum _{iI} V_text {dn}(vec {x}_i – vec {X}_I, varphi _i) , . end{aligned}$$

The disk-disk interaction is given by (V_text {dd}), the needle-needle interaction by (V_text {nn}) and the disk-needle interaction by (V_text {dn}). By integrating over the needle degrees of freedom one can derive the effective part of the interaction ({mathcal V} ({ vec {X}_I })) between the disks ^{66},

$$begin{aligned} beta mathcal V ( { vec {X}_I } ) = – ln left[ int prod _i {mathrm{d}} vec {x}_i {mathrm{d}} varphi _i exp left( – beta left[ sum _{iI} V_{dn}(vec {x}_i – vec {X}_I, varphi _i) + sum _{i<j} V_{nn}(vec {x}_i – vec {x}_j, varphi _i varphi _j) right] right) right] end{aligned}$$

((beta equiv 1/ k_text {B} T)). The corresponding force (mathcal F_K ( { vec {X}_I } )) on disk *K* is given by

$$begin{aligned} mathcal F_K ( { vec {X}_I } )&= – nabla _{vec {X}_K} mathcal V ( { vec {X}_I } ) = – frac{1}{N_text {n}} sum _l int rho ^{(1)}(vec {x}_l, varphi _l | { vec {X}_I }) nabla _{vec {X}_K} V_text {dn}(vec {x}_l – vec {X}_K, varphi _l) {mathrm{d}} vec {x}_l {mathrm{d}} varphi _l. end{aligned}$$

Here, we used the single particle density of needles with angle (varphi _l) at (vec {x}_l) for fixed disk positions:

$$begin{aligned} rho ^{(1)}(vec {x}_l, varphi _l | { vec {X}_I }) = N_text {n} frac{ int prod _{i ne l} {mathrm{d}} vec {x}_i {mathrm{d}} varphi _i exp left( – beta left[ sum _{iI} V_text {dn}(vec {x}_i – vec {X}_I, varphi _i) + sum _{i<j} V_text {nn}(vec {x}_i – vec {x}_j, varphi _i, varphi _j) right] right) }{ int prod _i {mathrm{d}} vec {x}_i {mathrm{d}} varphi _i exp left( – beta left[ sum _{iI} V_text {dn}(vec {x}_i – vec {X}_I, varphi _i) + sum _{i<j} V_text {nn}(vec {x}_i – vec {x}_j, varphi _i, varphi _j) right] right) } , . end{aligned}$$

By using ({nabla _{vec {X}_K} V_text {dn}(vec {x_l}-vec {X_K}, varphi _l) = – nabla _{vec {x}_l} V_text {dn}(vec {x_l}-vec {X_K}, varphi _l)}), evaluating the sum to a factor (N_text {n}) and defining the average over the needle angles as ({langle A rangle _varphi = int {mathrm{d}} varphi A(varphi )}), we obtain

$$begin{aligned} mathcal F_{vec {r}} ( vec {0}, vec {r} )&= leftlangle int rho ^{(1)}(vec {r}’, varphi | vec {0}, vec {r}) nabla _{vec {r}’} V_text {dn}(vec {r}’ – vec {r}, varphi ) {mathrm{d}} vec {r}’ rightrangle _varphi end{aligned}$$

for the case of two disks at (vec {0}) and (vec {r}). We use the superposition approximation ({rho ^{(1)}(vec {r}’, varphi | vec {0}, vec {r}) approx rho (vec {r}’, varphi | vec {0})rho (vec {r}’-vec {r}, varphi | vec {0}) / rho _text {n}}), where ({rho (vec {r}’, varphi ) equiv rho (vec {r}’, varphi | vec {0})}) is the density distribution around a single disk. For a single disk the needles are distributed according to the direct interaction potential (V_text {dn}(vec {r}, varphi )),

$$begin{aligned} rho (vec {r}, varphi )&= rho _text {n} exp left( – beta V_text {dn}(vec {r}, varphi ) right) . end{aligned}$$

This finally leads to

$$begin{aligned} mathcal F_{vec {r}} ( vec {0}, vec {r} )&approx nabla _{vec {r}} frac{1}{rho _text {n} beta } leftlangle int rho (vec {r}’, varphi ) rho (vec {r}’-vec {r}, varphi ) {mathrm{d}} vec {r}’ rightrangle _varphi = – nabla _{vec {r}}U_text {dep}(vec {r}) , . end{aligned}$$

This effective potential is the density-dependent depletion interaction, which we further approximate by employing a factorization approximation for the angular averages, ({langle rho (vec {r}’, varphi ) rho (vec {r}’-vec {r}, varphi ) rangle _varphi approx langle rho (vec {r}’, varphi ) rangle _varphi langle rho (vec {r}’-vec {r}, varphi ) rangle _varphi }), which is valid for the isotropic phase and the ideal nematic phase. Since we investigate the effective interaction in the nematic phase this should be a good approximation. This gives our final result

$$begin{aligned} beta U_text {dep}(vec {r})&approx – frac{1}{rho _text {n}} int langle rho (vec {r}’, varphi ) rangle _varphi langle rho (vec {r}’-vec {r}, varphi ) rangle _varphi {mathrm{d}} vec {r}’ = – rho _text {n} int left( 1 – frac{rho (vec {r}’)}{rho _text {n}} right) left( 1 – frac{rho (vec {r}’-vec {r})}{rho _text {n}} right) {mathrm{d}} vec {r}’. end{aligned}$$

More intermediate steps of the derivation are given in the Supplementary Information.