key: cord-0600309-rqno6zhg
authors: Kumari, Sunita; Ye, Fangfu; Podgornik, Rudolf
title: Ordering of adsorbed rigid rods mediated by the Boussinesq interaction on a soft substrate
date: 2020-07-19
journal: nan
DOI: nan
sha: 53705d1b20f170cbcf47c28b582c13f6f3a6b4fd
doc_id: 600309
cord_uid: rqno6zhg

Orientational ordering driven by mechanical distortion of soft substrates plays a major role in material transformation processes such as elastocapillarity and surface anchoring. We present a theoretical model of the orientational response of anisotropic rods deposited onto a surface of a soft, elastic substrate of finite thickness. We show that anisotropic rods exhibit a continuous {sl isotropic-nematic phase transition}, driven by rod-rod interactions mediated by the deformation of the underlying elastic substrate, and quantified by the Boussinesq solution adapted to the case of surface deposited rods. From the microscopic rod-rod interactions we derive the appropriate Maier-Saupe mean-field description, which includes the Boussinesq elastic free-energy contribution due to the substrate elasticity, derive the conditions for the existence of a continuous orientational ordering transition and discuss the implication of results in the soft (bio) systems context.

Living organisms interact sensitively with their environments in particular with numerous surfaces and substrates permeating their natural milieu. The surface compliance triggers varied responses that play key roles in regulating their behaviors and fates. For example, the novel coronavirus (SARS-CoV2) deposited on various surfaces can proliferate for varied amounts of time, depending on the nature of the substrate [1] . Spinal cord neurons extend more primary dendrites and shorter axons on stiffer gels [2] and astroglia are less adherent to the softer hydrogels compared to hard gels [3] . Rod-like fd virus particles form three distinct interaction regimes (linear chains, collapsed globules and chain-like aggregation) on lipid membranes depending on the adhesion strength [4] . Interestingly, stationary cells plated on an elastic substrate which is cyclically stretched re-orientate away from the stretching direction [5, 6] , and locomotive cells on a strained elastic substrate re-orient in the strain direction [7] .

As design and technology entails progressively more sophisticated manipulation of materials, the role of surfaces is becoming more prominent. In this context, a soft elastic substrate is of particular interest because of its flexibility and deformability. Its elastic deformation energy can mediate long-range elastic interactions between substrate-embedded particles [8] [9] [10] [11] that often gives rise to different structural motifs, such as an analog of elastic "Cheerios effect" in which spherical [12] and cylindrical [13] particles are suspended in an elastic gel. The long-range elastic interactions between particles deposited on a soft gel come about from the local surface depression of the substrate, created either by their weight or by adsorption energy, and leading to long-range attractive interaction [12, 13] . Similar interactions are found for liquid droplets on elastomeric interfaces that give rise to an inverted Cheerios effect [14] . The effect of the surface elasticity has been studied also in the context of wetting transition and adsorption phenomena [15] [16] [17] . Elastic substrates introduce complex liquid film dynamics with surface topography as the substrate shape changes with adsorption and in turn, this change in shape alters its wetting properties, consequently modifying the adsorption properties of the vicinal liquid [15] . Combined effects of surface elasticity and wetting have been demonstrated also in the context of ultrahydrophobic surfaces [18] [19] [20] and fast condensation on soft surfaces [21] . Elastic interactions in conjunction with defects in liquid crystals control the coordinates of embedded cylinders and tune the directional interactions [22] so that the trapped cylinders orient along the director axis on thin nematic films and pairs of cylinders form chains by capillarity, which forces them to align along the director by an elastic torque [22] .

The properties of membranes containing anisotropic inclusions were examined by considering the membrane as a continuous medium [23] . In the linear response approximation the membranes were shown to be weakly perturbed by inclusions of particles and to mediate long-range elastic interactions between the inclusions [8] . Using an analytical method, a wrapping -unwrapping transition of an elastic filament around cylinders has been studied in detail and a change of entropy due to wrapping of the filament around the adsorbing cylinders as they move closer together was identified as an additional source of interactions between them [11] . Recent computer simulation work shows that the rigid rod-like particles on an elastic substrate are predicted to show membrane-driven self-organization controlled by the tension and curvature of the underlying membrane [24] . An effective field theory was derived for the analytical computation of the membrane-mediated interactions. The formalism effectively reduces the rigid particles to points and computes the interaction free energy as an asymptotic expansion in inverse separations in a systematic way [25] . Within the framework of the free energy formulation, the intricate shapes of hydrogel menisci due to the indentation of rigid particle has also been analyzed and shown to result from a competition between surface tension, elasticity and hydrostatic pressure inside the gel, which drive the mutual interaction between particles [26] . In a very recent work, a mean-field theory (MFT) based on the minimization of a free energy functional was implemented to explain the wetting transition of elastic substrates [17] . The substrate was described in terms of the linear theory of elasticity, and the fluid contribution was based on the van der Waals theory. A long-range attractive interaction between fluid particles was caused by substrate elasticity, which moreover induced a decrease of the critical wetting temperature as compared to the case of wetting of an inert substrate, or precludes critical wetting of the elastic substrate altogether [17] .

However, despite significant progress in recent years on soft surfaces, the study of orientation of particles -especially in contact with deformable substrates -is comparatively scant. Moreover, apart from the nature of the substrate, the shape of the anisotropic inclusions implying perturbed orientational degrees of freedom plays and equally important role. In what follows, we will analyze the possibility of a substrate mediated, ordered or nematic-like phase of surface deposited rods. To our knowledge, no work has been performed yet to discuss an isotropic-nematic phase transition of rods deposited onto an elastic substrate. Studying such systems can provide insights into the influence of confined substrate deposition on physical properties of soft matter. We consider a system of anisotropic rods deposited on the surface of an elastic substrate of finite thickness. The substrate deforms under the loading of rods. We apply the Boussinesq theory [27] to calculate the deformation and consequent elastic displacements within the substrate. The Boussinesq theory describes the displacement equilibrium inside an elastic half-space in the presence of a surface inclusion, generating point-like forces normal to the surface [28] . We will generalize this theory to account for a finite thickness of the elastic medium as well as the rod-like shapes of the inclusions. As each rod generates a deformation field around it, the elastic coupling between different rods creates mutual interactions and hence long-range directional orientation forces. Including these interactions in a statistical field theory yields a continuous orientational phase transition of rods interacting via an orientational Boussinesq pair potential.

The paper is organized as follows: To demonstrate the range of the theory and for an intuitive illustration of the methodology, we first consider a single particle interaction with the elastic substrate in Sec. II. The phenomenology of the elastic deformation is described briefly and we derive explicitly the Green functions associated with an incompressible linear elastic solid by solving the pertaining Navier equation with the Galerkin Ansatz. In the next Sec. (III), the pairwise interactions between two rod-like particles, obtained by analogy with electrostatics [29] , are considered together with the model system of many-particle generalization. Afterwards, in Sec. (IV) we address the topic of mutual interactions of deposited rods, stemming from the elastic distortion induced in the substrate. We subsequently formulate an analytical theory, based on a field theoretical representation of the canonical partition function, that yields a MFT (saddle-point) solution for the orientational order parameter. Finally, we delineate some conclusions in the last Sec. (V) and in the Appendices present further computational details.

We start with the the problem of elastic deformations and stresses arising from a single, rigid point-like particle deposited on the surface of an elastic substrate of finite thickness, see Fig. (1). The later condition is not considered in the standard Boussinesq problem, which is formulated for a semi-infinite substrate [27] . The Boussinesq theory is based on elastic displacement potential function for a point-like normal force, acting at the surface, which satisfies the Navier equation and which can be used to find the displacement components [28] . A related problem of a tangential force acting at the surface is considered in the Cerruti problem [30] .

FIG. 1. The Boussinesq problem for a finite thickness elastic substrate. (a) The deformed surface profile (schematic) of a substrate for a point-like source F at r(x, y) acting at the upper z = 0 plane of a substrate. (b) Contour plot of the Boussinesq iso-uz (uz(r, z)) surfaces from Eq. (D10) for σ = 0.5. We plot uz(r, z), z and r in the units of (ρs(1 + σ)P )/2πE. Throughout this paper, we consider an isotropic and homogeneous elastic substrate of a finite thickness h. We assume a point-like rigid particle centered at (x, y, z = 0) = (r, z = 0), with r = (x, y) at the upper free surface of the substrate, chosen to be at the z = 0 plane. The z axis points downward (negative) so that points within the substrate have z < 0. The point-like force F = (0, 0, f ) imposed by this particle creates a stress in the substrate and as a result the substrate is deformed. The exact nature of this applied force is immaterial at this point. It can be either externally imposed, can be the weight of the inclusion, or can be a result of the adhesion energy between the particle and the surface, see Sec. (V) for discussion.

The elastic stresses are obviously distributed with cylindrical symmetry and we formulate the equations for the displacement field at any arbitrary point in the substrate due to the applied force. The surface density of the force in the direction normal to the surface is given by f = P δ 2 (r − r ′ ), where δ 2 (r) is the 2D Dirac delta function so that in extenso F = (0, 0, F z (r, z)) = (0, 0, P δ(z)δ 2 (r − r ′ )),

with P its magnitude. In the linear regime, the elastic displacement fields must satisfy the Navier equation [31] ,

where ρ s is the density of the substrate material, σ is the Poisson ratio related to the compressibility of the substrate and E is its elastic modulus. For isotropic materials, the Poisson ratio must satisfy −1 ≤ σ ≤ .5. In general, Eq.

(2) is difficult to solve because of the coupling of various components of the displacement vector. However, it can be solved straightforwardly by introducing the Galerkin vector representation method, which transforms the second order Navier differential equation into a fourth-order Galerkin equation (for detailed derivation see Appendix (B)) by introducing g [32] ,

in terms of which the Navier equation (Eq. (2)) assumes the form of a vector inhomogeneous biharmonic differential equation

Because of the cylindrical symmetry of the problem this equation can be simplified further by introducing the Love potential function [28] as,

in terms of which we remain with a scalar inhomogeneous biharmonic equation for the Love potential,

The connection between the displacement vector u, the stress tensor p ij and the Love potential can be derived straightforwardly (for detailed derivation see Appendix (B)). Since the external force is located at the upper surface of the substrate, it can be subsummed into a boundary condition, leaving a homogeneous biharmonic equation (Eq. 

with r = |r|, in terms of which the homogeneous Love equation is reduced to

and has a general solution: where the coefficients can be obtained from the boundary conditions Eqs. (7a) and (7b). The displacement field then follows as, see Eqs. (B2a) and (B2b),

The corresponding contours of the u z (r, z) profiles, obtained from Eqs. (11b), are displayed in Fig. (2) . Clearly the effect of the boundary condition at the lower z = −h surface of the substrate, is to "compress" the contours in the z-direction. The extent of the deformation due to the surface source is thus cut down mostly to the vicinity of the top layer. Building on the obvious analogy with electrostatics -a valid analogy since both theories are linear -but taking into account that the fields and sources are vectorial in the Boussinesq problem, one can obtain the resulting deformation vector for two point-like vector sources P 1 , P 2 acting at a fixed separation |r 1 − r 2 |, with

by simply superposing the solutions for the two single sources. This leads to the following contour plot of u z and u x,y in the (z, x) (or (z, y)) plane, see Fig. (3). Solving the Navier equation Eq. (2) with a Green's function dyadic G ik (r − r ′ ), for details see Appendix (A), and taking into account that the external force acts only in the z direction, one derives the elastic interactions free energy, in a completely analogous way as in the standard electrostatics, in the form

where dS is the surface integral and G zz (|r − r ′ |, z = z ′ = 0) is the zz component of the dyadic Green function for u z , so that we obtain from Eq. (A2) that the z-component of the total elastic deformation vector for the two point sources Eq. (12) is

where the indices 1, 2 stand for either of the two point source in Eq. (12) . G zz (|r 1 − r 2 |, z 1 = z 2 = 0) entering the interaction free energy, Eq. (13), can be obtained straightforwardly from the Love potential as . The combined deformation for the two sources u(r, z), mostly in the region between the sources, leads to a more negative total energy and thus to an effective elastic attraction between the sources. We plot uz(r, z), ux,y(r, z), z and r in the units of (ρs(1 + σ)P )/2πE.

where Z(r, z = 0) is the solution of Eq. (6) for a point-like surface force.

In the elastic interaction free energy, Eq. (13), we excluded the self-interaction terms that do not depend on the configuration of the sources. We refer to this interaction mediated by a finite slab of elastic material the Boussinesq interaction between two point-like sources. In what follows we will generalize this to the Boussinesq interaction between two rigid rods.

The dependence of the Boussinesq interaction on the separation between the two point sources is straightforwardly analyzed in terms of the function H σ (h/r), where r = |r 1 − r 2 |, connected to G zz (r, z = 0), see Eq. (E3), through

The dependence of H σ (x), where x = h/r, is shown in Fig. (4) . Obviously H σ (x) is a very slowly convergent function to its limit lim x−→∞ H σ (x) = 1. It is also evident that the vertical strain decreases with increase in the thickness of the substrate. Furthermore, the vertical strain also decreases with the increase in the distance from the axis of the load in the lateral direction at any depth. From Eq. (16) it follows that the Boussinesq interaction between two point surface sources is Coulombic for small separations but because of the finite thickness of the substrate it is cut off or screened for r ≥ h. This long-range nature of the interaction makes it distinct from, e.g., the excluded volume interaction which is short range. Also, notice the difference between 2D electrostatic interaction that depends logarithmically on separation, and the Boussinesq interaction which is also effectively 2D but varies as the inverse for small seprations.

We now turn to rigid rods of length L adsorbed onto the upper surface of the elastic substrate of finite thickness. An illustration of the instantaneous position and geometry of two rods is shown in Fig. (5) . The rods are identical in all aspects except for their orientation, and have by assumption a vanishing thickness. p(r) is now the linear density of the normal force acting along the rod and pointing downwards.

Each rod, i = 1, 2, is described by three parameters: the director, n, that defines the preferred orientation of the rod, assumed to be located in the plane perpendicular to the normal of the substrate surface, the center of mass located at r i , while s i is the contour parameter describing the position of the points along the rods in the frame of the rod itself, see the Fig. (5) , as For

runs over all the points along the i − th rod. The normal linear force density is then specified as

Each two points along the two rods are now assumed to interact via the Boussinesq interaction potential W 12 (|R 1 (n 1 , s 1 ) − R 2 (n 2 , s 2 )|). Following the analysis of Deutsch and Goldenfeld for the case of charged rods with screened Coulomb interaction, see Ref.

[29], we use the principle of superposition for the Boussinesq interaction to calculate the total interaction free energy between the two rods by summing over the point-wise interactions along their lengths. This leads to the total interaction free energy as

where G zz (r, z = 0) is again the Green's function as introduced in the previous section, see Eq. (15) . The above expression is a complicated integral over the contour parameters of the two rods. One could in principle carry on the analysis with this form of the interaction energy, but we choose a much more transparent development, based on the assumption that the separation between the rods is much larger then their length. This approximation is consistent with the fact that we will ignore the orientational Onsager interaction and concentrate solely on the long-range interactions.

In this limit we can expand the bracket term of Eq. (19) into a Taylor series up to the second order in the separation between the rods, yielding s1 s2

as all the odd powers of s 1,2 in the Taylor expansion integrate out to zero. After evaluating all the derivatives of the Love potential explicitly, what we remain with is an expansion of the interaction free energy in the form

where we introduced the following shorthand r = |r 1 − r 2 | andq 12 = (r 1 − r 2 )/|r 1 − r 2 |. As the above interaction free energy is quadratic in the director n, it describes the long-range nematic interaction. Explicitly, by introducing

The interpretation of the three different terms in the expansion is straighforward: in Eq. (22a) V (0) is the scalar interaction that depends only on the positions of the two rods but not on their orientations; in Eq. (22b) V (1) ij is the scalar-vector coupling interaction proportional to the square of the scalar product of the unit separation vector and the normal of the particle; in Eq. (22c) V (2) ijkl is the tensor interaction proportional to the product of the scalar products of the unit separation vector and the directors of the particles. The three terms together specify the elastic interaction free energy up to the order (L/r) 4 as the product Lp i can be introduced as the total "elastic charge" of the rod.

Since it will become clear that only the tensor interaction potential is relevant for the orientational ordering transition, taking into account Eq. (16) we deduce from Eq. (22c)

where x = h/r and H σ (x) has been defined in Eq. (E3). From Fig. (6) it is clear that F σ (x) has a behavior that exhibits a monotonic and an oscillating component for a range of small values of x, the relative importance of the two dependent on σ. The oscillatory component can actually lead to negative values of F σ (x). The orientational component of the interaction is thus non-monotonic for large lateral separation between the rods. In addition, for large values of the separation between the rods, r ≫ h, the interaction is effectively cut off.

Analyzing the different components of F σ (x) it is clear that apart from the first term they mostly cancel each other for large x, and that lim x→∞ F σ (x) = 90 lim x→∞ H σ (x) = 90. In any case, the strength of the tensorial Boussinesq interaction decays very steeply with the separation and is in addition cut off for large separations between the rods by the effects of the finite thickness of the substrate. 

In this section, we confine ourselves to a system of N indistinguishable rigid rods, differing only with respect to their orientations and positions, adsorbed onto an elastic substrate of finite thickness. The rods are oriented randomly on the surface in the X-Y plane at z = 0, see Fig. (7) . We assume that the interaction potential between rods is pairwise and is given by the Boussinesq interaction in the form valid for large separation between the rods, Eq. (21). The total elastic interaction energy between rods is then represented by

where all the interactions depend solely on V (i) (r µ − r ν ) ≡ V (i) (|r µ − r ν |). The Greek indices are used for particle identification and Latin indices are for Cartesian components of the vectors. The interaction potential is obviously complicated and we will introduce some additional approximations to simplify the analysis of the thermal properties of this system.

Instead of writing the configurational energy of the system in terms of particle coordinates and orientations, we will introduce two collective order parameters: the scalar particle density order parameter, ρ(r), and the traceless tensor nematic order parameter density, Q ij (r), where r = (x, y) is the position vector pertaining to the surface of the substrate. The two collective order parameters are defined as follows,

where the traceless quadrapolar nematic tensor order parameter for a two dimensional (2D) system is defined standardly asQ

FIG. 7. A schematic diagram of rods on an elastic substrate. Each rod is associated with a deformation field around it, leading to mutual interactions between them. Rods with their centers located at ri = (xi, yi, ) and with orientations in the (x, y) plane (represented by red arrows).

One notes the missing factor 3 2 in front of the first term in the above equation: this is due to the fact that the embedding space is 2D (surface of the substrate) and not 3D.

From its definition it is clear that the quadrupolar nematic tensor order parameter satisfies the following identitieŝ

Q µ ijQ µ kl = n µ i n µ j n µ k n µ l − 1 2 δ ij n µ k n µ l − 1 2 δ kl n µ i n µ j + 1 4 δ kl δ ij .

With the two collective order parameters we can now write the pairwise interaction energy of the system, Eq. 24, in a succinct form

Substituting the definitions of the collective order parameters, Eqs. (25), we are clearly back to the coordinate form Eq. (24). The first term in the above equation is the self-energy of the system defined as

which formally diverges, but we need to remember that we left out the short-range repulsions that would regularize it. The self-energy term corresponding to the Boussinesq interactions is completely analogous to the self-energy of the Coulomb system, and its only consequence is to renormalize the chemical potential [33] . The other terms in Eq. (29) are obtained as where V (0) , V (1) and V (2) have been defined in Eqs. (22a), (22b) and (22c), and depend only on the absolute value of the argument. In Eq. (29) the various area integrals contain an orientational part and radial part:

where ϕ is the polar angle in cylindrical coordinates and dS = 2πrdr is the area element. Clearly the polar angle integrals in Eq. (29) can be evaluated explicitly as the orientational dependence of the interaction free energy is also explicit. This can be done by noticing that in 2D

Furthermore, we assume that the bulk state of the system is homogeneous, so that the order parameter densities do not depend on the coordinate. Thus we finally remain with a much simplified form for the interaction energy

where S is the area of the elastic substrate and

where a (i) are the various orders of the second virial coefficient: a (0) is the positional virial coefficient, a (1) is the positional-orientational coupling virial coefficient and a (2) is the orientational virial coefficient. It is the latter that will play a fundamental role for the orientational ordering transition in what follows. Taking into account that the nematic order parameter density is traceless, we finally remain with

This has the form consistent with the Lebwohl-Lasher interaction energy [34] and will much simplify later calculations. The orientational phase transition depends only on the second term of the free energy and thus we only need to calculate a (2) which equals

where r = |r − r ′ |. In Appendix (F) we derive the final form of the dimensionless second virial coefficient, based on Eq. (F3),ã

and F σ (x), was defined in Eq. (23) . The orientational second virial coefficient, which will be the determining factor in our analysis of the orientational ordering transition below, then assumes the form, Eq. (F3),

The non-monotonic structure of the second virial coefficient, see Fig. (8) , is a direct consequence of the effective orientational interaction F σ (x), but is inconsquential in the macroscopic limit of h/a ≫ 1. At this point it is illuminating to compare the present theory with the 2D problem of thin rods interacting only through excluded volume interactions. The Onsager free energy can be conveniently expressed in the form Eq. (39) [35] , where the 2D virial coefficient scales as a (2) ≃ k B T L 2 [36] . Clearly the virial coefficient pertaining to the Boussinesq interactions is much larger for the same length of rods, and ignoring the Onsager interactions is thus entirely valid in this case.

What we have accomplished in this section is to write the interaction energy due to Boussinesq elastic deformation of the substrate in the form containing only the collective order parameters. This leads naturally to the introduction of different second virial coefficients connected with the Boussinesq interaction, a (0) and a (2) . From here we then need to proceed to the evaluation of the appropriate partition function of the canonical system.

To determine the nature of the equilibrium state of a system consisting of N rods, interacting via the Boussinesq elastic interactions, we must evaluate the corresponding canonical partition function defined as

where β = (k B T ) −1 is the inverse thermal energy and the product is over all N rods and the integration is over all configurations. We have shown that the interaction energy can be written in terms of collective order parameters ρ(r) and Q ij (r), defined in Eq. (25) , in the form

We now introduce two decompositions of unity

with the shorthandρ

and cast the functional delta function with its integral representation

where we designate φ(r) as the scalar molecular field (note that in the case of Coulomb fluids this would be the fluctuating electrostatic potential, see Ref. [33] ) and φ ij (r) as the nematic tensor molecular field. We can then rewrite the partition function as

This formidable looking expression can be reworked a bit to yield

where we introduced the partition function of the system in an scalar and tensor external field as

as the trace is diagonal and the single particle partition function in the external fields is theñ

From hereZ

with the effective field action defined as

This is the final form of the canonical partition function for a system of N rods interacting via the Boussinesq interaction. We note here that the derivation of the field representation of the partition function follows closely the related developments in the field representation of the Coulomb fluid partition function [33] , the main difference being that the orientational part of the Boussinesq interactions between rods brings fourth not only the scalar collective order parameter but also the tensorial collective order parameter.

Since the field action is highly non-linear there is no direct way to evaluate Eq. (53). Further approximations are therefore necessary. Just as in the case of the Coulomb fluid, where the saddle point approximation leads to the well known Poisson-Boltzmann theory, the saddle point approximation of Eq. (52) will naturally yield an equivalent of the Maier-Saupe theory. In the next section we will see how.

We first invoke the approximation of a homogeneous "bulk", meaning that the collective order parameters do not depend on the position along the surface. By invoking Eq. (39) and Eq. (53) this yields the field action as

where, to cut down on the clutter, we expressed the virial coefficients in terms of the thermal energy. In addition we can write (see Eq. (51))Z

where the orientational average has been defined before, see Eq. (34) . The saddle-point approximation is now defined as the stationary point of the field action, leading to the following two scalar equations and two tensorial equations

From the first pair of Eq. (56a), we see that the molecular field φ at the saddle point is purely imaginary, φ −→ iφ, just as in the case of the saddle-point approximation for Coulomb fluids [33] . The second identity then implies

where c R is the surface density of the rods (notably different from ρ s , the density of the elastic matrix). Let's concentrate on Eq. (56b) then. The first equation implies that at the saddle point φ ij −→ iφ ij , so that it is also purely imaginary. The second equation, after cancelling the imaginary unit, then becomes

where we introduced the orientational average

In Eq. (58) we recognize the Maier-Saupe approximation that describes the equilibrium properties of the system. We now assume that the saddle-point solution of the molecular field φ ij , from the first identity of Eq. (56b), has the form

where S is proportional to the orientational order parameter, which quantifies how well-aligned the molecules are, andn is the director of the ordered phase. Therefore

Multiplying this with (n inj − 1 2 δ ij ) on both sides, taking into account thatn · n = cos ϕ and contracting over the tensorial indices leads us to

where we took into account that the orientational second virial coefficient is negative. Except for the dimensionality, 2D vs. 3D, this equation is equivalent to the Kleinert formulation of the Maier-Saupe theory [37] . As it describes a lyotropic transition, the strength of nematic interaction quantified by the Boussinesq second virial coefficient is multiplied by the density of the rods, c R , i.e. the larger the local density of the rods, the smaller needs to be the orientational interaction driving the isotropic-nematic transition. Notably, the above solution is just a minimum of an effective free energy

where we took into account the integral representation of the modified Bessel function of the first kind I 0 (z). The solution of equation Eq. (62), see Fig. (9) , exhibits a continuous ordering transition at a critical value of (a (2) ρ) * ≃ 32.04, characterized by a non-vanishing value of S. Interestingly enough, the only formal difference between the 2D and the 3D cases is that instead of I 0 (S/2) in Eq. (63), one ends up with Exp(− 1

is the Dawson function [37] . However the solution of the Maier-Saupe equation then leads to a first order transition, with a finite jump in S.

In this article, we studied theoretically the orientational ordering of anisotropic rigid rods deposited onto the surface of a soft substrate of finite thickness. The substrate is assumed to be a homogeneous and isotropic elastic medium. The profile of a deformed substrate and the corresponding elastic displacements due to adsorbed rods are calculated by the linear Boussinesq theory for a finite thickness substrate. In this model the rods are represented by a concentrated force normal to the substrate surface, uniformly distributed over their length with a constant linear density. The longrange, Coulomb-like effective interactions between two rods, mediated by the elastic displacement of the substrate, are then pairwise additive on the linear level. At the macroscopic length scale we formulate the partition function of a canonical ensemble of rods interacting via these substrate elastic deformation mediated Boussinesq interactions and derive a mean-field type theory, yielding a continuous isotropic-nematic transition.

The model of a constant, normal force linear density along the adsorbed thin rods, similar to the model of rodlike polyelectrolyte rods with uniform linear charge density [29] , is a coarse grained and simplified version of the physical reality in the sense that, microscopically, there is an adhesion energy distributed over the area of the finite diameter rod in contact with the substrate. The local substrate surface deformation is then very similar to the case of adhesion of colloids to a soft membrane [10, [38] [39] [40] or wrapping of colloids along an elastic filament [11] and could be computed in the same way. In scaling terms the product of the normal force and the surface depression of the substrate surface should be proportional to the product of the adhesion (wrapping) surface energy density times the contact surface area. While in principle the calculation of effective rod-rod interactions with explicit adhesion to the surface is possible, it would make the analysis of a canonical ensemble of such rods unfeasible.

We first restrict ourselves to interactions between two point-like surface forces acting at r 1 and r 2 , respectively, and discuss the influence of substrate thickness on the corresponding strain fields in Sec. (II). The dependence of the vertical strain u z (r, z) on the substrate thickness is plotted in Fig. (2) , where different curves correspond to different values of the thickness, h(= 12, 6, 3). We set the Poisson ratio σ = 0.5 in all three cases. The surface pressure on the top surface of the substrate, z = 0, is set by the normal force density distribution of two point like surface forces, while deformation vectors at the bottom surface are constrained to vanish at z + h = 0. We found that the vertical strain decreases with increasing substrate thickness, while for a fixed thickness, it decreases with the radial distance, r. Importantly, it can be seen that as the substrate thickness decreases, the iso-deformation contours shrink towards the z = 0 surface. For example, for the coordinate z = -1.5 and r = 0 the u z (r, z) increases by a factor of 0.85P/(2πE) for h = 3 and 1.68P/(2πE) for h = 12, which implies that the iso-deformation contour contraction for h = 3 is significantly larger (about 50%) than that for h = 12. Fig. (3) depicts the iso-deformation contours of u z (r, z) and u x,y (r, z) due to the two point-like sources at dimensionless separation |r 1 − r 2 | = 4 for σ = 0.01. Because on the Boussinesq level the problem is linear, the deformation fields around each point-like source force are obtained form a superposition principle and in turn result in an effective elastic attraction between the sources. The situation is similar to electrostatics, the differences stemming from the fact that electrostatics is a scalar theory while elasticity is a tensor (vector) theory, thus the change in sign in the interaction.

It is obvious from Fig. (3) that the thinner the substrate, the higher the value of the strain. The separation dependence of the Boussinesq interaction between two point sources as a function of the separation between them is investigated in Sec. (II). Fig. (4) shows the dependence of the interaction scaling function H σ (h/r) on h/r for four different values of Poisson ratio σ =(0.5, 0.01, -0.5, -1.0). Clearly, H σ (h/r) saturates asymptotically to 1 as h/r → ∞ and approaches 0 for r >> h. This yields an effective Boussinesq interaction decaying as the inverse power of the separation just as in 3D electrostatics. For smaller value of h/r, dramatic changes are observed, with H σ (h/r) exhibiting a cut-off dependent on the thickness of the substrate. This seems to describe a "screening" effect of the finite thickness of the substrate, meaning that if the two point surface sources are close enough the interaction between them is Coulombic but is then effectively cut off for r h.

In Sec. (III), we discuss the substrate deformation caused by two adsorbed rods as well as the effective Boussinesq elastic interaction between them, as a function of their separation and orientation. Fig. (6) , shows the variation of the orientational nematic interaction in terms of the scaling function F σ (h/r). We notice that for small separation between rods (r << h), F σ (h/r) asymptotically approaches a constant value. In general, F σ (h/r) is composed of a monotonic and a superimposed oscillatory component that yields a negative interaction energy if the rods are away from each other (large value of r). It is clear that the rod-rod interaction decay sharply with the separation between the rods, and in fact for small enough h/r, F σ (h/r) limits to 0. It is worth mentioning that the material properties of the substrate have also a significant impact on the interaction between surface inclusions mostly through the Poisson ratio. As shown in Fig. (6) , the non-monotonic interaction component of F σ (h/r) is more pronounced for σ = 0.5 as compared to σ = 0.

In Sec. (IV), we extend our model to a canonical N rod system. The interaction potential between rods is assumed to be pairwise and is given by the Boussinesq interaction, leading to elastic interaction energy between rods given by Eq. (24) . We investigate the critical behavior of the orientational order parameter and the free energy through a solution of the Maier-Saupe equation that we derive analytically and exactly. For different values of the product of the orientational virial coefficient and the macroscopic surface density of the rods, a (2) c R in Eq. (63), as illustrated in Fig. (9(a) ), we conclude that for higher values of a (2) c R , there is a continuous lyotropic transition from an isotropic to a nematic phase. This happens at the critical value a (2) c R = 32.04, at which the orientational order parameter decreases sharply but continuously to zero, implying a second order transition. In order to have a nematic surface phase, according to Eq. F3, the elastic parameters of the system need to satisfy

where F s = ρ s P = ρ s (pL) is the total normal force per rod that can be interpreted as an "elastic charge" of the rod. While the value of this force is a phenomenological parameter of our model, the above limit is fulfilled even if it is very small. In fact taking reasonable values for the Young modulus of a gel, length of the rod on the order of 100 nm, and a on the order of nm, we are led to

which is easily fulfilled even for small surface concentrations of rods. We can also investigate the ratio between the Onsager excluded volume interaction in 2D and the orientational virial coefficient [36, 41, 42] . This ratio scales as

being much larger then 1 for reasonable values of the parameters. This also implies that the long-range Boussinesq elastic interactions are more important for surface ordering on soft substrates then the standard (Onsager) steric effect.

At higher values of a (2) c R (dense system) the stable phase is thus a nematic phase and the transition to this nematic phase is lyotropic at constant temperature. This second order lyotropic transition is purely an effect of a 2D nature of the problem. The same analysis in 3D would yield a first order transition, see [37] . These results are in line with previous studies of continuous istropic-nematic transition in 2D [43] [44] [45] [46] [47] [48] . A continuous transition was predicted by the 2D version of the Maier-Saupe theory for cylinders interacting via a weak anisotropic pair potential [43] and for particles in strong confinement [45, 46] . In addition, a new type of 2D nematic, a tetratic phase was hypothesized, by using an interaction of the Maier-Saupe type with terms representing both twofold and fourfold anisotropic interactions [47] . It was found that the isotropic-tetratic transition is always second-order and independent of the symmetry breaking parameter κ (the ratio between fourfold and twofold interactions). The tetratic-nematic transition is second-order for small κ, but becomes first-order at the tricritical point [47] . In a recent work, a first-order isotropic-nematic transition is observed for intermediate aspect ratio and a small anisotropic interaction parameter, while a continuous isotropic-tetratic transition is found for small aspect ratio and small anisotropic interaction parameter, followed by a continuous tetratic-nematic transition at higher densities [48] .

In conclusion, we would like to emphasise that the deformations and the mechanical response of elastic substrates significantly alters the directionality of the adsorbed rod-rod interactions and strongly influences their phase behavior. Throughout this work, we aimed to understand the governing parameters that tune the phase behavior of adsorbed rods on soft elastic surfaces. These results are therefore likely to have a significant impact on the design of smart surfaces or diagnostic tools for the detection of pathogenic molecules such as those based on surface binding of rodshaped viruses and bacteria. Our study can be used as a model system to understand their biological counterparts. For example, responses of cells to surface elastic forces, or surfaces of hydrogels are just a few to name. In these cases a dynamic generalization of our model, possibly along the lines of the recent work by Bar-Haim and Diamant [49] , would yield a rich enough model to be able to describe that kind of phenomenology. we obtain from Eq. (A4) the interaction free energy between these two point like sources in the form W = − 1 2 ρ 2 s G zz (r 1 − r 2 , z = z ′ = 0) P 1 P 2 ,

where we ignored the two formally infinite self-energy terms that do not modify the interaction between the sources. This is exactly Eq. (13) in the main text. The whole derivation very much follows the standard approach in electrostatics.

expressed in terms of the Fourier-Bessel components of the stress tensor. The first two identities of Eq. (C2) are satisfied by

while the last one implies and Eq. (C4),

At the lower surface, z = −h, the boundary condition is simply that all the components of the deformation vector vanish at this plane. is non-vanishing and an explicit solution can therefore be obtained.

Proceedings of the National Academy of Sciences of the United States of America

Proceedings of National Academy of sciences, United State of America

Application des potentielsá lëtude de léquilibre et du mouvement des solidesélastiques

Theory of Elasticity

Collective Classical and Quantum Fields

We start with minimization of the elastic free energy in the presence of external forces, see Ref. [31] for details, that implies the Navier equation, Eq. (2), aswhich, since it is linear, can be solved by the Ansatzwhere the elastic Green's dyadic has to satisfyInserting the Green's dyadic into the elastic free energy we remain withIf the external force is given by two point like sources acting at the upper surface, z = 0, in the z direction, then with Eq. (12)Appendix B: Elastic displacement and stress tensorIn this appendix, we provide explicit expressions for the elastic displacement and stress tensor with the Galerkin-Love Ansatz [32] , The Galerkin vector and the Love solution are introduced as followsFrom here, the elastic displacement components are expressed asandFrom standard definitions the stress tensor the off diagonal components p xy , p xz , p yz then follow aswhile the diagonal components readwhere the arrows indicate the two quantities have the same form except for the substitution indicated by the arrow.

Consistent with the cylindrical symmetry of the problem and isotropy in the (x, y) plane, one can introduce the Fourier-Bessel transform Z(Q = |Q|, z), see Eq. (8) , so that the homogeneous Love equation, Eq. (9), is solved bywhere the integration constants A, B, C and D are determined from the boundary conditions, which we consider next.At the upper surface, z = 0, the boundary condition Eqs. (7b) and (7a) specify the components of the stress tensor. As the force is applied only in the z− direction the boundary condition at that surface reads,Appendix D: Solution of the boundary conditions and the Boussinesq limit The values of coefficients A, B, C and D in Eq. (C12) are then found to be as follows,The Love potential can then be evaluated from Eq. (C1) asThis is the solution of the Boussinesq problem for a finite thickness elastic matrix with a point force acting perpendicular to the upper surface. Note that in the limit h −→ ∞ only the coefficients A and B remain finite, so that the original Boussinesq solution is recoveredExpressing the displacement field with the Love potential, Eqs. (B2b), we getand analogously for u x,y . Inserting the solution Eq. (D5) we are left withand again analogously for u x,y with Eq. (B2b). From here we get back the solution u z (r, z) to the original Boussinesq problem of a point force at the origin with h −→ ∞,wherefrom it followsand analogously for u (x,y) (r, z)where we wrote the solution with the components of the Green's dyadic, see Eq. (A2). It is clear form this solution that the elastic deformation diverges at the origin for r = 0, z = 0. This is of course an artifact, as the Boussinesq calculation works only for linear elasticity valid for small deformations, while the deformation close to the application point of the external force can be large. This spurious divergence can be eliminated by assuming a finite region of the application of the force and truncating the solution within this region. In Fig. (2) we display the contours of the z component of the deformation vector, u z (r, z). Clearly the effect of the boundary condition at the bottom surface, viz., u x,y,z (r, z = h) = 0 is the contraction of the profile in the z direction, see the contour plot in Fig. (2) . This affects also the interaction potential between the point sources or consequently the extended rods.

What is an observable is the deformation vector at the top surface, whose z component we consider next. The deformation u z of the top surface is now obtained from the Green's function G zz defined in Eq. (15) . This Green's function at the same time defines also the interaction between two point-like surface inclusions. In the Boussinesq limit of an infinitely thick substrate, Eq. (D10), we can extractwhich is of an electrostatic form. We now proceed to the general case of finite thickness of the substrate where we getFrom here we obtainwhere H σ (x) fully quantifies the dependence on the thickness of the substrate on top of the Coulombic dependence for an infinitely thick substrate, Eq. (E1). The dependence H σ (x) is shown on Fig. (4) . From this functional form it is clear that for large thickness of the substrateand we are back to the Coulombic dependence for the pure Boussinesq case, Eq. (E1). On the other hand for x ≤ 0.65 H σ (x) changes sign and approaches zero. In the case of two surface point sources at r 1 and r 2 , the elastic interaction between them is proportional to G zz (|r 1 − r 2 |, z = 0), see Eq. (13) . We thus conclude that the interaction between two point surface sources is almost Coulombic for small separations but is then effectively cut off for ρ ≥ h. It looks as if the interactions are "screened". This effect is due purely to the finite thickness of the substrate.Appendix F: Calculation of a (2) We start from Eq. (40) and Eq. (23), where we derived explicitlyThis furthermore implies thatwhereã (2) is the dimensionless orientational virial coefficient and F σ (x) has been defined in Eq. (23), while the infinite upper bound has been substituted by h/a, where a is the microscopic cutoff scale, corresponding roughly to the molecular size. This cutoff is necessary since otherwise the integral is divergent. The asymptotic limit of the orientational virial coefficient for large values of the argument is lim h a −→∞ãwhich follows from the definition of the integrand F σ (x) and its properties in the asymptotic limit. As the integrand is also a very slowly limiting function, its limit is also reached only asymptotically. The above results are used in the analysis of the orientational virial coefficient in the main text.