Orientational order and glassy states in networks of semiflexible polymers

Motivated by the structure of networks of cross-linked cytoskeletal biopolymers, we study the orientationally ordered phases in two-dimensional networks of randomly cross-linked semiflexible polymers. We consider permanent cross-links which prescribe a finite angle and treat them as quenched disorder in a semi-microscopic replica field theory. Starting from a fluid of un-cross-linked polymers and small polymer clusters (sol) and increasing the cross-link density, a continuous gelation transition occurs. In the resulting gel, the semiflexible chains either display long range orientational order or are frozen in random directions depending on the value of the crossing angle, the crosslink concentration and the stiffness of the polymers. A crossing angle $\theta\sim 2\pi/M$ leads to long range $M$-fold orientational order, e.g.,"hexatic"or"tetratic"for $\theta=60^{\circ}$ or $90^{\circ}$, respectively. The transition is discontinuous and the critical cross-link density depends on the bending stiffness of the polymers and the cross-link geometry: the higher the stiffness and the lower $M$, the lower the critical number of cross-links. In between the sol and the long range ordered state, we always observe a gel which is a statistically isotropic amorphous solid (SIAS) with random positional and random orientational localization of the participating polymers.


I. INTRODUCTION
The cytoskeleton is a network of linked protein fibers which plays an important role in several functions of eukaryotic cells such as maintenance of morphology, mechanics and intracellular transport [1]. Cytoskeletal fibers, such as F-actin can be described as semiflexible polymers with a behavior intermediate between the two extreme cases of rigid rods and random coils. The function of the actin cytoskeleton is modulated by a large number of actin-binding proteins (ABPs) [2,3]. The organization of actin filaments into networks is regulated by ABPs which can be classified into two broad categories; cross-linkers promote binding of the filaments at finite crossing angles whereas bundlers promote formation of bundles consisting of parallel or antiparallel filaments. F-actin is a polar semiflexible polymer and some ABPs bind filaments with a specific polarity whereas others are not affected by the filament's polarity. In order to unravel the physics of these complex aggregates, in vitro solutions of actin filaments with controlled cross-linkers have been studied [4].
The stiffness of semiflexible filaments gives rise to orientational correlations and promotes the formation of structures with long range orientational order. The isotropic-nematic transition in solutions of partially flexible macromolecules has been studied theoretically using the Onsager approach [5] or the Maier-Saupe approach [6,7]. In Ref. [7], the role of the solvent is taken into account and a very rich phase transition kinetics is predicted. Aggregation and orientational or-dering of Lennard-Jones macromolecules with bending and torsional rigidity, with or without solvent, has recently been investigated in Ref. [8] using molecular dynamics simulations. Experimental investigations of the isotropic-nematic transition in lyotropic F-actin solutions have been carried out in Ref. [9] and measurements of the associated order parameter are presented in [10].
The excluded volume effect is not the only mechanism which drives an isotropic-nematic transition in solutions of semiflexible polymers. Assemblies of cytoskeletal filaments exhibit a structural polymorphism due to interactions mediated by a wide range of ABPs, as shown in Ref. [11]. Theoretical attempts to study the formation of ordered structures in this kind of systems involve a generalized Onsager approach [12][13][14], a Flory-Huggins theory [15], and a semi-microscopic replica field theory [16]. In Refs. [12][13][14][15] the filaments are modeled as perfectly rigid rods whereas in Ref. [16] they are considered to be semiflexible polymers and the thermal bending fluctuations (finite persistence length) are fully taken into account.
In vitro structural studies by Wong et al. [17] of Factin in the presence of counterions have shown the formation of sheets ("rafts") with tetratic order without any cross-linking proteins. It appears that electrostatic interactions are the main mechanism behind the effective "π/2" cross-linking of the actin filaments in this experiment. A theory for tetratic raft formation by rigid rods and reversible sliding "π/2" cross-linkers has been proposed by Borukhov and Bruinsma [13].
In Ref. [16], a three-dimensional melt of identical, fixed-contour-length semiflexible chains is considered and random permanent cross-links are introduced. The cross-links are such that they fix the relative positions of the corresponding polymer segments and constrain their orientations to be parallel or antiparallel to each other. The aim of the present study is to examine the effect of crosslinks which prescribe a finite crossing angle on a twodimensional version of the previous model. We distinguish the case of sensitive cross-links that are perceptive of the polymers' polarity and unsensitive cross-links that are not. In contrast to previous studies where the rigid rods have a finite width which allows for nematic orderingà la Onsager, our semiflexible polymers are considered one-dimensional objects and the sole cause for the emergence of orientationally ordered phases is the interplay of the finite persistence length of the polymers and the cross-link geometry.
Using the polymer stiffness and the cross-link density as control parameters, we obtain a phase diagram which involves a sol and various types of orientationally ordered gels. For appropriate values of the control parameters and crossing angle θ ∼ 2π/M ; M ∈ Z, we predict the emergence of an exotic gel with random positional localization and M -fold orientational order (e.g., hexatic for M = 6 or tetratic for M = 4). Similar phases have been predicted [14] for a different system in three space dimensions: a semi-dilute solution of charged rods with finite diameter in the presence of polyvalent ions that have the function of (non-permanent) cross-linkers and may favor various crossing angles. Besides the orientationally ordered phase, we also predict a statistically isotropic amorphous solid (SIAS) with random positional and orientational localization of its constituent polymers appearing right at the gelation transition.
This article is organized as follows. In Section II, we introduce our model and the Deam-Edwards distribution which parametrizes the quenched disorder associated with the cross-links. In Section III, the disorder-averaged free energy is presented as a functional of a coarsegrained field which plays the role of an order parameter. Then, in Section IV, we discuss different variational Ansätze that express positional and orientational localization. The symmetries imposed by the cross-linking constraints allow the emergence of specific orientationally ordered phases. The corresponding free energies are calculated variationally in the saddle-point approximation and a phase diagram is obtained. We summarize and give an outlook in Section V. Finally, details of the calculations are given in the appendices.

II. MODEL
We consider a large rectangular two-dimensional volume (area) V which contains N identical semiflexible polymers. A single polymer of total contour length L is represented by a curve in 2d space with r(s) denoting the position vector at arc-length s. Bending a polymer costs energy according to the effective free energy functional ("Hamiltonian") for the wormlike chain (WLC) model [18] (1) Here we have introduced the tangent vector t(s) = dr(s)/ds and chosen a parametrization of the curve, in accord with the local inextensibility constraint of the WLC, such that |t(s)| = 1. The position vector is recovered from t(s) as r(s) = r(0) + s 0 t(τ )dτ . Hence the conformations of a single polymer, that can be bent but not be stretched, are alternatively characterized by t(τ ), 0 ≤ τ ≤ L, and r(0). The bending stiffness is denoted by κ; it determines the persistence length L p according to κ = L p k B T /2. Throughout the rest of the paper we set k B T = 1. The WLC model can describe stiff rods, obtained in the limit L/L p → 0 and fully flexible coils, obtained in the limit L p /L → 0.
The mutual repulsion of all monomers is modeled by the excluded-volume interaction being the Fourier transformation of the positionalorientational density. t i (s) = (cos ψ i (s), sin ψ i (s)) denotes the orientation of monomer s on polymer i. The coefficients λ 2 |k|,m depend only on the absolute value of the vector k in order to preserve the rotational symmetry of the system. They are later chosen large enough to provide stability with respect to density modulations. For details, see Appendix C.
The Hamiltonian H W LC is invariant with respect to interchanging head and tail of the filaments, i.e., the energy is unchanged under reparametrizations of the contour of one polymer i by {r i (s) → r i (L − s), ∀s ∈ [0, L]}. However, in the following, we want to consider filaments with a definite polarity which could for example arise due to the helicity of F-actin. The WLC Hamiltonian is not sensitive to such a polarity, but the cross-links may or may not differentiate between the two states of the filament, as discussed below.
We now introduce M permanent (chemical) cross-links between pairs of randomly chosen monomers. A single cross-link, say between segments s and s ′ on two chains i and j, constrains the two polymer segments to be at the same position, i.e. r i (s) = r j (s ′ ), and fixes their relative orientation by a constraint expressed by the function ∆ (t i (s), t j (s ′ ), θ). The constraint of relative orientation is most easily formulated in polar coordinates where ∆ just fixes the crossing angle ψ − ψ ′ to some prescribed value θ. In the case of cross-links sensitive to polarity we simply have where the weight of the delta function is such that 2π 0 dψ 2π δ(ψ) = 1. In the unsensitive case on the other hand, changing head and tail of either filament leads to four equivalent situations that correspond to two different crossing angles (see Fig. 1). In this case the cross-link (4) thus appears with the two equally likely cross-linking angles θ and θ + π. We point out that in our model the polarity is only recognized by the cross-links and does not enter into the Hamiltonian. The delta function in (4) is the simplest way to model the angular constraint of the cross-links. It is however much more realistic to introduce an effective angular cross-linking potential that allows for fluctuations around the preferred direction. A simple model for these "soft" cross-links is given by The cross-link stiffness parameter γ controls the variance of the fluctuations of the angle and it is clear that in the limit γ → ∞ we recover the simple delta function model of "hard" cross-links. In the following we will first explore the case of hard cross-links, because it is mathematically simpler. It gives rise, however, to artifacts which disappear when considering the more physical model which favors certain angles but does not enforce them strictly. The cross-links are permanent and do not break up or rebuild. Hence we are led to study the equilibration of the thermal degrees of freedom {r i (s)} in the presence of quenched disorder represented by a given cross-link configuration C = {i e , j e ; s e , s ′ e } M e=1 which is characterized by the set of pairs of polymer segments that are involved in a cross-link. The central quantity of interest is the canonical partition function Here the thermal average < ... > H is taken with respect to the Hamiltonian H := H W LC + H ev and the partition function depends on the cross-link realization C. The free energy F is expected to be self-averaging so that we are allowed to compute the disorder averaged free energy .] denotes the "quenched" average over all cross-link configurations according to some distribution P(C). We assume that the different realizations obey a Deam-Edwards-like [19] distribution implying that polymer segments that are likely to be close to each other in the uncross-linked melt also have a high probability of getting cross-linked. Note that the probability of cross-linking is independent of the relative orientation of the two segments, so that the cross-link actually reorients the two participating segments. The parameter µ 2 controls the mean number of cross-links in the system: Given the Hamiltonian, the constraints due to crosslinking and the distribution of cross-link realizations, the specification of the model is complete and we proceed to calculate the disorder averaged free energy [F ].

III. REPLICA FREE ENERGY
The standard method to treat the quenched disorder average is the replica method, representing the disorder averaged free energy as [F ] = −[ln Z] = − lim n→0 ([Z n ] − 1)/n, i.e., we have to calculate the simultaneous disorder average of the partition sums of n non-interacting copies of our system. In the end, we extract the disorder averaged free energy [F ] from [Z n ] as the linear order coefficient of its expansion in the replica number n.
In this section, we present only the essential formulas and present the details of the calculation in Appendix B. The disorder average, [Z n ], gives rise to an effectively uniform theory, however, with a coupling of different replicas: where the average ... H n+1 is over the (n + 1)-fold repli-cated Hamiltonian H. To simplify the notation we have introduced the abbreviationsr ≡ r 0 , r 1 , . . . , r n , ψ ≡ ψ 1 , . . . , ψ n , and the shorthand notation s ≡ (1/L) L 0 ds. Furthermore, δ(r) ≡ n α=0 δ(r α ) and ∆ denotes for sensitive cross-links and (10) if they are unsensitive to polarity of the filaments. For soft cross-links we introduce the corresponding definitions.
Different polymers are decoupled via a Hubbard-Stratonovich transformation, introducing collective fields, Ω(x,φ), whose expectation values are given by The field Ω quantifies the probability to find monomer s on chain i at position x 0 in replica 0, at position x 1 in replica 1,... and at position x n in replica n and similarly to find it oriented in the direction e 1 = (cos(ϕ 1 ), sin(ϕ 1 )) in replica 1,... and oriented in the direction e n = (cos(ϕ n ), sin(ϕ n )) in replica n. Sometimes it is convenient to use its equivalent representation in Fourier space In terms of these collective fields the effective replica theory is given by The replica free energy per polymer reads with the effective single chain partition function The average ... W LC n+1 refers to the (n + 1)-fold replicated Hamiltonian of a single wormlike chain. The sum k over replicated wave vectors is restricted to the combination of the so-called "0-replica sector" (0RS) that contains only the pointk = (0, . . . , 0) and the "higherreplica sector" (HRS) where at least two wave vectors in different replicas are non-zero: k α = 0 and k β = 0 with α = β.
Note that in the above overview we left out the "1replica sector" (1RS) that consists of (n + 1)-fold replicated vectorsk α where only one entry (the αth) is nonzero, i.e.k α = (0, ..., k α , ..., 0), andm being arbitrary. The corresponding fieldsΩ α (k,m) describe regular macroscopic density fluctuations (modulated states). In this work we focus on the properties of the macroscopically translationally invariant amorphous solid state and assume that the inter-polymer interactions (2) are such that periodic density fluctuations are suppressed. See Appendices B and C for more details on that issue.

IV. VARIATIONAL APPROACH
We shall only discuss the saddle-point approximation to the field theory on the right hand side of (13) replacing the field Ω by its saddle point Ω sp which has to be calculated from δF /δΩ| Ωsp = 0. In fact, even the saddle point equation is too hard to solve, because the conformational distribution of the WLC is very complex. Similarly, it is not possible to perform a complete stability analysis of the Gaussian theory. Hence we restrict ourselves to a variational approach and in the following we are going to construct Ansätze which capture the symmetry of the different physical states which may emerge.
What behavior do we expect? Upon increasing the number of cross-links up to about one per polymer, there should be a gelation transition from a liquid to an amorphous solid state. A finite fraction of the polymers form the percolating cluster and are localized at random positions. The other polymers belong to finite clusters or remain un-cross-linked and are still free to move around in the volume. This scenario has been found to be valid for a variety of different models [20][21][22][23][24]. A similar replica field-theoretic approach has beeen used by Panyukov and Rabin to study the properties of well cross-linked macromolecular networks [25,26].
For semiflexible polymers the positional localization in the macroscopic cluster is accompanied by orientational localization. The cross-links create locally an orientational structure. Supposing that the semiflexible polymers are rather stiff, a long range orientationally ordered state can be established. Otherwise we may find an orientational glass, i.e., the directions of the polymer segments are frozen in random directions in analogy to the low temperature phase of a spin glass [27]. We call such a phase statistically isotropic amorphous solid (SIAS).
Can we expect long range orientational order for arbitrary crossing angles θ? We first consider a special case, namely that the crossing angle is such that an integer multiple of the crossing angle, M θ, equals 2π, e.g. θ = 120 • and M = 3. Such a choice of crossing angle allows for orientationally ordered states provided the chains are sufficiently stiff. In the case of sensitive cross-links we expect M -fold discrete rotational symmetry, in the case of unsensitive cross-links and odd M we expect 2M -fold symmetry because cross-links are established including both the angles θ and θ + π (see Fig. 1). Sketches of gels with long range four-fold order (θ = 90 • ) or long range three-fold or six-fold order respectively (θ = 120 • / 60 • ) are shown in Fig. 2 and Fig. 3. The case θ = k 2π M with k = 1, . . . , M − 1 leads in a similar fashion to an M-fold symmetric phase.
Let us now consider the constraints on the order parameter imposed by symmetries. To this end we note that the Hamiltonian (1), the cross-link constraints (6) and the Deam-Edwards distribution (7)  Only the fluid state is expected to have the full symmetry of the partition function Eq.(8) because here, the polymers and finite clusters of polymers are free to sample the complete volume and can take any orientation. This implies for the order parameter i.e., each replica is separately invariant under translations and rotations. We will only investigate gels with random localization of a fraction of the particles. In other words, we restrict ourselves to amorphous solids and do not consider periodic density fluctuations. It might be of interest to study one-dimensional periodic density modulations with an overall orientation of the WLCs perpendicular to the wave vector,-reminiscent of bundles. However this is not the topic of the present paper, where we stick to the incompressible limit.
If the translational and rotational symmetry is broken spontaneously, then we expect a non-trivial expectation value of the local density in a particular equilibrium state. For example, in the amorphous solid phase a finite fraction of particles should be spatially localized at preferred positions and possibly oriented along preferred directions. A simple Ansatz quantifying such a scenario is the following: The vector a is the preferred mean position of monomer s and the thermal fluctuations around the preferred position are controlled by the localization length ξ. The preferred orientation is given by ϕ 0 and η parametrizes the variance of the orientational distribution.
An amorphous solid is macroscopically translational invariant, so that the spontaneous symmetry breaking of translational invariance is local. All macroscopic observables should display translational symmetry and in particular all moments of the local densitỹ are non-zero only, if the wave vectors add up to zero (see reference [21] for a more detailed discussion).
As far as the rotational symmetry is concerned, we will consider two possibilities: a statistically isotropic state (SIAS) as well as states with true long range orientational order. In the former case the spontaneous symmetry breaking of rotational invariance is local and restored globally analogously to what happens in the case of translational symmetry. All macroscopic properties, such as the moments, are non-zero only, if the angular momentum "quantum" numbers add up to zero: For the long range ordered case the simplest Ansatz consists in assuming a local M -fold symmetry. This implies for the moments i.e. each m α has to equal an integer multiple of M .
Note that a rotation affects both r i (s) and t i (s) and that consequently, a system with e.g. M -fold symmetry should be described by an Ansatz Ω(k,φ) which is only invariant under common rotations of the k α and ϕ α . We have chosen our approach for simplicity and incorporated the symmetry with respect to the orientational and position vector separately. Consequently, the Ansatz is symmetric under individual rotations of the spatial and angular argument. An interesting extension of the present work would involve the study of Ansätze which couple positional and orientational degrees of freedom. Some preliminary work in this direction has been done in the context of random networks of covalently connected atoms or molecules [28,29].
Rephrasing our results in replica language, frozen fluctuations in a single equilibrium state correspond to a non-zero expectation value of the density, ρ i,s (k, m), within one replica. The full statistics of the local static fluctuations that we need for the characterization of the macroscopic symmetries in the gel is specified by the order parameter field Ω(k,m), encoding moments of arbitrary order. In particular macroscopic translational invariance requires n α=0 k α = 0; macroscopic rotational invariance requires n α=1 m α = 0; long range orientational order for crossing angle θ = 2π M requires m α = lM with l ∈ Z.
Let us now encode these physical expectations in a variational Ansatz: The liquid state is characterized by Ω sp (k,m) = δk ,0 δm ,0 . It becomes unstable at a critical cross-link density, µ 2 c , where a percolation transition takes place and a macroscopic cluster is formed containing a finite fraction of the polymers. To account for the fraction of localized particles Q in the percolating cluster coexisting with mobile particles (fraction 1 − Q) in finite clusters, we make the following general Ansatz for the expectation value of the order parameter: The first term describes the sol phase which is characterized by perfect spatial homogeneity and orientational isotropy. The second part describes the amorphous solid which is macroscopically translational invariant as reflected in δ 0, n α=0 k α . The details of the amorphous solid phase under consideration have to be implemented in the order parameter ω(k 2 ,φ).
Plugging the generic form of the order parameter (22) into the replica free energy (14) we get We expect the gel fraction to be small close to the gelation transition and thus expand the log-trace contribution of the free energy in Q. The terms linear in Q cancel, as they should for the expansion around Q = 0 to be justified.
A. M -fold orientationally ordered amorphous solid

Hard cross-links
Assuming that the thermal fluctuations of the polymers are to a certain degree suppressed by their stiffness and a sufficient number of cross-links has been formed, long range orientational order may be present as sketched in Fig. 2 and Fig. 3. The Ansatz for the M -fold orientationally ordered amorphous solid favors M preferred orientational axes separated by angles 2π M . A simple way to incorporate this symmetry is the following Ansatz for the replica order parameter Here I 0 denotes the modified Bessel function of the first kind; it ensures the proper normalization. Fork =0, the above order parameter is a probability distribution and thus specifies the local orientational order completely. In experiment on the other hand one has access to low order moments only. The simplest physical order parameter being sensitive to the degree of long range M -fold orientational order is given by where in the second line, we replaced the average over all monomers of the system by the disorder average of one arbitrary monomer. This expression allows to relate S M to our effective one-particle theory. . . . denotes the thermal expectation value of the system for a given instance of disorder C. In order to establish the connection between the variational parameter η and the physical order parameter S M we evaluate (26) using (24) and find that S M is in leading order proportional to η.
We plug the Ansatz (24) in (23) and perform thek summations andφ integrations. If we have chosen M according to the symmetry considerations of the previous section the constraint ∆ drops out. We are left with the free energy as a function of three variational parameters: the gel fraction Q, the spatial localization length ξ and the degree of orientational order as measured by η. To determine these parameters we minimize the free energy. The resulting equation for the gel fraction is universal and has been derived previously [21]. There are two solutions Q = 0 and Q ∼ 2 µ 4 (µ 2 − 1) which implies that the phase transition from sol to gel takes place at µ 2 c = 1. In order to determine the remaining parameters 1 ξ 2 and η we assume them to be small near the transition and do a Taylor expansion of the free energy where we leave out terms which do not depend on ξ or η because these terms are irrelevant for the variation of F . More precisely, it will turn out that 1 ξ 2 ∝ Q, being thus small close to the transition, and we will keep contributions up to order 1/ξ 2 in our expansion. As for the variational parameter η, we will find that it actually jumps from 0 to a finite value, i.e that there is a first order orientational transition. We expand all the same in η and obtain at least a qualitative picture. Details on the calculations are given in Appendix D.
Since our Ansatz is replica-symmetric it is straightforward to extract the part of F linear in the replica index n and we find Note that for M > 1 there are no terms coupling the orientational order specified by η and the positional localization specified by ξ 2 . Such terms appear in higher order but will not be considered here because we restrict ourselves to the vicinity of the gel point. The persistance length L p and the localization length ξ are both rescaled by the contour length L of the polymers.
The functions g, h,h,h and l,l,l go to zero for large argument and are for small argument approximately given by For the definitions of these functions see Appendix F.
Minimizing F for M > 1 with respect to ξ 2 yields Hence the filaments are localized as soon as a percolating cluster of cross-linked chains has formed and the gel fraction is finite. The localization length is independent of M and its scale is set by the radius of gyration R g of the filament which is ∼ L 2 for stiff chains and ∼ L for random coils. For finite polymer flexibility L Lp the incipient gel has no long range orientational order. Increasing the number of cross-links beyond µ 2 = 1 we find a first-order transition that corresponds to a minimum of the free energy at nonzero η, which is first metastable and eventually becomes the global minimum (see Fig. 4). The critical values of µ 2 for the first appearance of a metastable ordered state (µ 2 1 ), the first order phase transition (µ 2 2 ) and the disappearance of the metastable disordered state (µ 2 3 ) are shown in Fig. 5. In the limit L Lp → ∞ the phase transition coincides with the sol-gel transition, whereas for more flexible filaments higher cross-link densities are required. For rather stiff polymers we find that the transition takes place at In addition, Fig. 5 shows the dependence of η c , which is the value of the variational parameter η at the phase transition, on the polymer flexibility. There is a tendency to a lower degree of orientational localization for higher values of the polymer stiffness. This behavior is qualitatively different from the behavior predicted for the lyotropic nematic ordering of partially flexible rods by Khokhlov and Semenov in the corresponding range of flexibility [5]. However, a direct comparison of the two models is not feasible. In our case, orientational ordering is due to the cross-links and the critical density of cross-links µ 2 2 of the transition decreases for increasing polymer stiffness. There are thus two competing effects: stiffer polymers should lead to a stronger orientational localization of the polymers whereas the smaller number of cross-links should have the opposite effect. It seems that the lower number of cross-links plays the dominant role in the dependence of η c (L/L p ).
We point out that the cross-linking angle θ = 2π/M enters the free energy as a rescaling of the polymer flexibility L/L p through the parameter M 2 . This implies that the higher the value of M (the smaller the angle), the higher the polymer stiffness required for the transition at a given cross-link density. The reason for this scaling is that within our mean-field description we are dealing with one single chain in an effective medium. Mfold order has to be propagated along the chain from one cross-link to the next and the M 2 scaling reflects the properties of the WLC in the calculation of the corresponding correlator.  where a metastable ordered minimum appears, where it becomes the global minimum (phase transition) and where the minimum at zero corresponding to the disordered state becomes unstable (plot for M > 1). Note that always µ 2 1 < µ 2 2 < µ 2 3 for finite stiffness.
For a first-order transition the expansion in η 2 is not really justified and can only give qualitative results. Even worse, the expansion of F in η 2 exhibits an oscillating behavior: the orders 2, 6, 10, . . . provide a stable orientational free energy (for large η) in a region around µ 2 = 1 and L Lp = 0, but the orders 4, 8, 12, . . . are always unstable in the above region because they diverge asymptoti-cally to minus infinity. However, considering the purely orientational free energy it is obvious that for any finite κ the Gaussian part is larger than the log-trace contribution as long as the crosslink density µ 2 does not become too large. So, as long as we restrict ourselves to a region close enough to the sol-gel transition asymptotic stability is garantied and a qualitative picture can be obtained by truncating the expansion at order 6, 10 or even higher. For the polar case, i.e M = 1, there are additional terms ins the free energy (27) which couple spatial and orientational part. At first sight it might be tempting to argue that they can be neglected close to the transition because they are proportional to 1 ξ 2 ∼ Q. As it turns out this is not correct: the orientational transition for a given polymer stiffness is shifted to significantly higher values of µ 2 . But all the same, the qualitative picture of the transition is still valid.
In order to analyze the orientational transition we first calculate the stationarity equation with respect to 1/ξ 2 and plug it again into the free energy. Keeping contributions up to order η 6 we obtain again a power series in η 2 . The transition scenario that we found for M > 1 is still valid, but the transition takes place at larger values of µ 2 2 . For rather stiff polymers we find approximately µ 2 2 ∼ 1 + 0.94 L/L p (33) and the critical µ 2 is considerably larger than what we would obtain from (30) without the coupling terms. At first sight, it seems a bit curious that the case M = 1 is set apart by its coupling term. This is, however, only due to the low order expansion and the coupling terms for M > 1 have not appeared yet. Because of the higher symmetry of the orientational contributions, only higher order spatial contributions may lead to nonvanishing coupling terms.
If long range orientational order could exist only for rational values of θ/2π, then it would be inaccessible in experiment. We show in the next section that the long range ordered states discussed above are also present if we introduce crosslinks that favor given crossing angles on average only.

Soft cross-links
Which changes do we expect when using soft crosslinks (5) that do not rigidly fix the intersection angle to one particular value, but allow for fluctuations about a given mean direction? Soft cross-links surely will make it harder to establish long range order, so the first expectation is that the orientational transition will need a higher cross-link density µ 2 to take place. In particular, the behaviour in the limit of stiff rods will change qualitatively: In the case of hard cross-links the macroscopic network is a completely rigid object, even at cross-link densities just above the critical value of µ 2 = 1, whereas in the case of soft cross-links there are still the degrees of freedom of fluctuations around the preferred directions of the crosslinks left. Instead of approaching a combined spatial and orientational transition right at µ 2 = 1 upon increasing the polymers stiffness, it appears possible that for soft cross-links the long range order transition will take place not at µ 2 = 1 but at a higher cross-link densitiy because the networks needs to be stabilized. However, in the limit of stiff rods and hard cross-links we should recover the combined transition right at µ 2 = 1.
We now discuss the modifications of the free energy due to the soft cross-links: Keeping terms only up to order Q 3 as before, we find that in the case M > 1 there is still an effective decoupling of spatial and orientational transition. The spatial part does not change at all, but the purely orientational contribution becomes f or,sof t := (34) The soft cross-links give rise to additional factors IMq(γ) I0(γ) that equal 1 in the limit of hard cross-links γ → ∞. They appear linearly in the Gaussian and quadratically in the contribution of the log-trace term. Their range lying between 0 and 1 it is clear that the log-trace contribution is smaller with respect to the Gaussian so that the transition occurs at a higher cross-link concentration (as compared to the case of hard cross-links). This is confirmed by a numerical analysis of the Taylor expansion up to 10th order in η.
The effect of soft cross-links on the orientational phase transition is illustrated in Fig. 6 for M = 3: First of all, it shows that for finite γ the phase transition takes place at a finite distance from µ 2 = 1, even in the limit of stiff rods. Moreover, comparing the curves corresponding to increasing values of γ, i.e. to harder and harder orientational cross-link constraints, the curves converge towards the solid curve at the bottom that was drawn for perfectly hard cross-links as they should.
The case M = 1 involves additional coupling terms between spatial and orientational parameters that need to be calculated. But as before, they don't lead to a behavior that differs qualitatively from what we found for M > 1. In the preceding section we found that the higher M the more cross-links are needed to get into the long range ordered phase. This holds true for soft cross-links, too, as we checked numerically. The related calculations are presented in Appendices D and F.

B. Statistically isotropic amorphous solid (SIAS)
In the preceding two sections, we found that (given a suitable cross-linking angle θ) there is a phase boundary µ 2 ( L Lp ) above which long range orientational order becomes possible. For lower cross-link density, long range orientational order vanishes. But the positional localization of the polymer segments that takes place in the macroscopic cluster is always accompanied by orientational localization. The corresponding alternative to long range order are glassy states, where the average orientations of the polymer segments are frozen in random directions, so that isotropy is restored on a macroscopic level. We expect such states for WLCs with a small persistence length, such that the order induced by a cross-link cannot be sustained along the contour length up to the next cross-link. This will be particularly severe, if M is large.
Frozen orientations for polar filaments are described by the distribution (18) for a single site. However the direction of localization varies from chain to chain, so that averaging over the whole sample implies an average over the locally preferred orientation ϕ 0 assuming all directions to be equally likely. For convenience we introduce the unit vectors u α = (cos ϕ α , sin ϕ α ) and denoting the local preferential axes by e we get de 2π exp ηe · ( n α=1 u α ) = I 0 η| n α=1 u α | .
Altogether the order parameter for polar filaments in the glassy state then reads What is the physical order parameter and how does it scale with the parameter η?
Locally, we have for each localized polymer segment a polar moment t i (s) ∼ η = 0, but because of the orientational disorder it vanishes globally. We thus need an Edwards-Anderson like order parameter and choose As for the M -fold order parameter in the preceding section, we relate the variational parameter η to q and find that in lowest order We plug the Ansatz (36) for the replicated order parameter into the saddle point free energy and keep only terms which depend on ξ or η. We obtain the simplest non-trivial free energy by including terms up to second order in η 2 .
where we have introduced the shorthand notation Λ(γ) := I 2 1 (γ) I 2 0 (γ) . The functions g, h, l are given by (28) in the M-fold section. Minimizing the above free energy with respect to ξ 2 , yields qualitatively the same result as for the long range ordered state, namely The orientational part shows a behavior different from the M-fold case: The stationarity equation with respect to η 2 gives rise to The coupling terms implies a non-zero value η 2 , i.e. orientational localization, as soon as positional localization sets in. Hence glassy orientational order is enslaved to positional localization and the orientational order parameter grows continuously at the gelation transition. Varying γ we see that the softer the cross-links are, the smaller is also the variational parameter η, i.e. the degree of orientational order. Note that, in contrast to the M -fold ordered case, the limit n → 0 leads in the SIAS free energy to an extra minus sign for all the orientational contributions because of the coupling of the replicas in the corresponding order parameter (36). As a consequence, we have to maximize the free energy with respect to η instead of minimizing it as it is well known from spin glasses [27].

C. Phase diagram
The results of the previous sections can be summarized in a phase diagram, presented in Fig. 7. The control parameters are the cross-link density measured by µ 2 , the polymer flexibility measured by L/L p , and the crosslinking angle θ. Irrespective of the cross-linking angle, independent of the stiffness of the filaments and of the softness of the cross-links there is a continuous gelation transition accompanied by random local orientational ordering (SIAS phase) at the critical cross-link density µ 2 c = 1. This glassy ordering has been encountered previously for randomly linked molecules with many legs (p-Beine) [20,22]. The free energy of the SIAS is above the free energy of the sol, which however is unstable beyond µ 2 = 1 and hence, is not available in this region of the phase diagram.
What is new, is the appearance of a state with long range orientational order, if the crossing angle θ = k M 2π where k, M ∈ Z. This phase is characterized by a spontaneous breaking of the rotational symmetry. For sensitive cross-links the symmetry of the orientational order is M -fold, for unsensitive cross-links and odd M the resulting phase has 2M -fold symmetry. The free energy of the long ranged ordered state is below the free energy of the isotropic sol as well as below the free energy of the SIAS. Hence, we expect it to win as soon as it appears, even though we cannot do a complete stability analysis beyond the variational Ansatz.
We find that the appearance of long range order is pushed to higher values of µ 2 as the constraint for the crossing angle is softened. The same effect is observed for increasing flexibility of the polymers, -because it becomes more difficult to sustain the orientation of the polymers -, and increasing M .
When reading the phase diagram, we should keep in mind that the parameters µ 2 and L/L p are not thermodynamic ones like a temperature or a chemical potential because changing either of them changes the disorder ensemble, too. This means that two points in the phase diagram correspond effectively to two different systems.

V. CONCLUSIONS -DISCUSSION
In this paper, we studied the role of the cross-linking angle in the formation of orientational order in random networks of semiflexible polymers in two dimensions. We have used a variational Ansatz to map out a phase diagram. Besides a statistically isotropic amorphous solid (SIAS) we find more exotic gels with random positional order coexisting with long ranged orientational order, whose symmetry is dictated by the crossing angle. In analogy to liquid crystals with long range orientational order and thermal centre of mass motion like in a fluid, these gel phases might be termed "glassy crystals"-with long range orientational order and frozen in random positions like in a glass. It is interesting to note that the tetratic ordering, which corresponds to M = 4 in our system, has been predicted and/or observed in a variety of physical systems with quite different constituents and underlying physical mechanisms [14,[30][31][32].
Because of the peculiarities of two dimensions, we expect fluctuations to affect positional localization [33,34]. It is known, in the context of 2d defect-mediated melting, that at finite temperature positional order can only be quasi-long ranged whereas orientational order can be truely long ranged [35]. A study of the corresponding phenomena for our positionally amorphous and orientationally ordered system is a very interesting direction for further investigation.
We point out that the finite bending rigidity is an essential ingredient of our model. It allows the effective decoupling of positional and orientational degrees of freedom close to the gelation transition. In the case of infinitely stiff polymers on the other hand, a rigid crosslink would automatically fix both position and orientation. The nature of the emerging network is an interesting problem which goes beyond the scope of this paper.
Another possible extension of our work concerns more elaborate variational Ansätze probing the appearance of combined M-fold and glassy orientational order. Here, the chains in the gel fraction are assumeed to be orien-tationally localized in preferential directions which vary from chain to chain but macroscopically average in an M-fold pattern.
The three-dimensional generalization of our model has to deal with the fact that a finite cross-linking angle between two wormlike chains prescribes a cone and not a plane. If we want to describe cross-links with torsional rigidity, we need to go beyond the simple WLC model and use the helical WLC [36].
In this work, we focused on orientational and glassy order in networks of semiflexible polymers mediated through appropriate cross-links. We assume that the excluded-volume interaction is such that it supports a macroscopically translationally and rotationally invariant liquid which, upon cross-linking, may give rise to ordered networks. Although this assumption is mathematically consistent and facilitates the analytical treatment of our model, it may be challenged in experimental realizations where the excluded volume may lead to lyotropic alignment in dense systems. In a future extension of our model, one may envisage adding a Maier-Saupe aligning pseudopotential as in Refs. [6,7] and studying its interplay with the cross-link induced interaction. The wormlike chain propagator G(ϕ, s 1 ; ϕ ′ , s 2 ) quantifies the probability that the tangential vector of monomer s 1 points into the direction ϕ provided that the tangential vector of monomer s 2 points into the direction ϕ ′ : Here, N denotes the normalization and ψ(s) the angle corresponding to the unity vector t(s). In principle, the path integral includes all the monomers from 0 to L, but the monomers that lie not between s 1 and s 2 do not affect the result and can be integrated out. In order to perform the remaining path integral we write down a discretized version of the above expression replacing the continuous degrees of freedom by a finite number l of them such that t 1 corresponds to t(s 1 ) and that t l corresponds to t(s 2 ). The distance ǫ between neighbors is determined by lǫ = |s 2 − s 1 |. Expressing the integral and the derivatives of the wormlike chain Hamiltonian by their discretized versions and calling the normalization constant for the discretized path integral N ǫ we arrive at [37] We want now to perform the integrations over the ϕ i . For that purpose it is convenient to decouple ϕ i and ϕ i+1 by means of where the I q (a) denote modified Bessel functions. Performing the integrations we arrive at and integrating over ϕ1 2π and ϕ l 2π we find that the normalization is given by N ǫ = exp(−κ/ǫ)I 0 (κ/ǫ) l−1 .
Last, we take the limit ǫ → 0 keeping lǫ = |s 1 − s 2 | constant. It is thus possible to express the modified Bessel functions by the asymptotic expansion I q (a) ∼ exp(a) √ 2πa 1 − 4q 2 −1 8a + . . . and the propagator converges to Calculating a general 2-point correlation function of the (real valued) observables O 1 and O 2 at positions s 1 and s 2 respectively by means of the above propagator we find whereÔ 1 andÔ 2 denote the Fourier transformation of the observables. We learn from this expression that only the Fourier components to the same q couple to each other and that the rate of the exponential decay of correlations along the filament scales with q 2 .
Appendix B: Mean-field replica free energy and the saddle point equations In order to obtain the disorder averaged free energy [F ] by means of the replica method we need to calculate the disorder averaged n-fold replicated partition function where r α i (s) and ψ α i (s) denote the position vector and angle of orientation of segment s belonging to polymer i inside the αth replica. For sensitive cross-links θ σ equals always the single crossing-angle θ, i.e. the normalization y = 1, and we can omit the summation over σ, but in the unsensitive case θ σ takes the two values θ 1 = θ and θ 2 = θ + π and so, y = 2.
Observing that the formula factorizes in the cross-link index e it is possible to perform the sum over the number of cross-links M that leads to an exponential function where for sensitive cross-links the function ∆ is defined as ∆ s (ψ, θ) ≡ e γ n α=1 cos(ψ α −θ) and in the unsensitive case on the other hand as ∆ u (ψ, θ) ≡ 1 2 e γ α cos(ψ α −θ) I n 0 (γ) + e γ α cos(ψ α −(θ+π)) I n 0 (γ) .
(B.3) After the disorder average, all the sites, i.e. all the polymer segments, are equivalent but still appear explicitly, coupled by the cross-linking constraint ∆ and the delta functions. Expressing the delta functions in Fourier space we can rewrite our formula in terms of the quantity Using this definition and writing H ev explicitly the replicated partition function reads In section II the parameter λ 2 was introduced in its most general form depending on |k| and m, but as it turns out a constant is sufficient for our purpose and allows in the following for a more compact notation. Because of the symmetry |Q(k,m)| 2 = |Q(−k, −m)| 2 it is only the real part of ∆m that contributes and we thus redefine the sensitive kernel ∆ s accordingly as For unsensitive cross-links ∆ u,m equals zero if the sum n α=1 m α is not even, but takes otherwise the values of the sensitive kernel in (B.4).
We are now going to rearrange the contributions of intra-polymer repulsion and Deam-Edwards distribution into contributions belonging to the following subsets of the space of (n + 1)-fold replicated vectorsk: • The 0-replica sector (0RS) consisting only of0 In the following we denote the sum over 0RS and HRS by k and the sum over the 1RS as˜ k . For the 1RS we obtain the new kernel UsingQ α (k,m) as shorthand notation for Q(k,m) when only the wave vector k α in replica α is non-zero the expression reads In order to decouple the sites we are now going to apply a Hubbard-Stratonovich transformation. The symmetry Q(k,m) = Q * (−k, −m) and the analogous relation for Q α (k,m) have to be reflected by their corresponding fields Ω(k,m) andΩ α (k,m). After the transformation [Z n ] can be written as where the integrals over the complex fields are meant to be integrations over real and imaginary parts separately.
The replica free energy F is given by Using the saddle point approximation we replace the integral (B.6) by its maximal contribution. The corresponding values of the order parameter fields have to fulfill the self-consistency equations that arise from Note that there are many ways to perform the HS transformation, each of them leading to a different replica free energy F . The resulting saddle-point replica free energy on the other hand is always the same as it can be checked easily. As a consequence, we are free to redefine the fields for later convenience by doing the replacements ℜΩ(k,m) → √ ∆m ℜΩ(k,m) and for the imaginary part accordingly. The advantage of this transformation is that the saddle-point fields Ω SP (k,m) are now directly related to an expectation value of Q(k,m) as presented in the main part, equation (12). This connection between Q and its corresponding field Ω is derived by means of an external field that couples to Q. Performing then the Hubbard-Stratonovich transformation and taking the logarithmic derivatives with respect to real and imaginary part of the field on both sides of the equation, i.e. for both the microscopic and the field theoretic representation of our theory, results in the desired relation between Q and Ω.
The 0RS/HRS part of the free energy in terms of the new fields reads then The 1RS part is treated in the next section.
Appendix C: Stability of the 1-replica sector The fields of the 1-replica sector,Ω α (k,m), can be subdivided into the fields wherem is such that only the entry m α is non-vanishing and those corresponding to the other possible values ofm. In the former case the fields describe M -fold orientationally symmetric density fluctuations with wave vector k in replica α and have thus a clear physical meaning. The other fields break replica symmetry in the sense that they describe density fluctuations e.g. in replica α accompanied by purely orientational fluctuations in other replicas. These fields are unphysical and need thus to equal zero.
In order to study the stability of the 1RS with respect to fluctuations in the physical fields ρ α (k, m) it is sufficient to study stability for one replica only. The expansion of the corresponding free energy reads up to second order with r x := r(s x ) and ψ x := ψ(s x ). Choosing λ 2 > µ 2 the kernel∆ m = λ 2 2 − µ 2 2 Im(γ) I0(γ) cos(mθ) is positive because Im(γ) I0(γ) ≤ 1. It is thus possible to redefine the fields by introducing ρ(k, m) := √ ∆ mρ (k, m) without changing the stability. The corresponding free energy is given by Let us now consider the two contributions to F 1RS separately: It is evident that the first term is a positive quadratic form provided that λ 2 is large enough. If we consider stiff rods as the limiting case of very unflexible wormlike chains the matrix of the second quadratic form can be diagonalized by switching back from the Fourier modes m to the real space variables ϕ. Writing r(s) = r 0 + st and integrating over r 0 and t we find for C s1,s2 and know thus that the second quadratic form is positive semi-definite. Altogether we have thus proved stability in the case of stiff rods, but we assume that this result will at least hold for wormlike chains with large L/L p , too.
Appendix D: Variational free energy: M-fold symmetric long ranged order In this Appendix, we present some details on the derivation of the variational free energy (27) of the Mfold orientationally ordered amorphous solid. We insert the corresponding replica order parameter Ansatz (24) into the general replica free energy (23) and expand it in powers the variational parameters 1 ξ 2 and η.

Gaussian part
The Gaussian part reads The spatial part is easily computed by replacing the sum over the replicated Fourier variablesk by an integral and representing the delta function in Fourier space. Performing the resulting integrations we arrive at For the orientational contribution f o we find For θ = k M 2π the cosine equals 1 and vanishes. In order to obtain the last line we replaced the three functions depending on the ϕ's by their Fourier representations. Performing then the integrations lead to Kronecker Deltas that cancel two of the three sums over Fourier modes.
Altogether we find for the Gaussian contribution in linear order in the replica index n Considering unsensitive crosslinker the corresponding expression reads (D.2) i.e. for odd M the contributions corresponding to odd q are projected out. For M even on the other hand, we obtain the same result as for the sensitive cross-links.
Note that for hard crosslinks (γ = ∞) there is an alternative expression for the f o . For sensitive cross-links it is given by and in the unsensitive case we have ) .

Log-trace contributions
Let us now turn to the log-trace contribution that we will expand up to the third order in the gel fraction Q: The first-order term gives a trivial contribution because of In the second step we user(s) =r(0) + s 0 dτt(τ ) and transform the original path integral D{r(s)} into an integral dr(0) V n+1 and a path integral D{ψ(s)} over angular variables. Performing ther(0) integration we get a Kronecker delta settingk to zero. The last line follows from the fact that the functions are normalized and the two integrations with respect toφ andψ(s) simply lead to 1.
The second-order contribution reads .
Before Taylor expanding the expression in the variational parameters, we need to perform the summations overk 1 andk 2 . We integrate over dr(0) as before and get a Kronecker Delta imposingk 1 = −k 2 . There is thus only one summation left. We replace this sum by an integral and represent the Kronecker delta in Fourier space. After the two Gaussian integrations we find cos(∆ϕ α 1 )+cos(∆ϕ α 2 ) {cos(M(ψ α 1 +∆ϕ α 1 ))+cos(M(ψ α 2 +∆ϕ α 2 ))} I 2n 0 (η) using the shorthand notations f α ≡ s2 s1 ds t α (s) and ψ α i := ψ α (s i ). This expression can be expanded up to the desired order and the remaining task then consists in calculating the corresponding correlation functions. Details are given in Appendix F.
For the calculation of the third-order term, we proceed the same way as before and need thus to perform three sums over wave vectors. With the abbreviation f α l ≡ s l s3 dτ t α (τ ) the result reads cos(∆ϕ α 1 )+cos(∆ϕ α 2 )+cos(∆ϕ α 3 ) .
Appendix E: Variational free energy: SIAS

Gaussian part
In the following, we present the derivation of the SIAS variational free energy. Let us first turn to the Gaussian contribution. It is given by where u α = (cos ϕ α , sin ϕ α ). The spatial part is treated the same way as for the M-fold case, but the orientational contribution is slightly more complicated because it does not factorize in the replica index. By means of the integral representation of the Bessel function, In order to obtain the last line, we simply switch from the angular variables to Fourier space. The expression is then easily expanded up to the desired order.

Log-trace contributions
The first-order contribution of the log-trace part gives the same (trivial) result as for the long range ordered case. As for the second and third order term, we again have to perform thek summations first, but this calculation is done exactly the same way as before.The corresponding results are obtained from those of the M-fold case by simply replacing the M-fold orientational distributions in (D.7) and (D.8) by the corresponding SIAS distributions.
Details on the calculation of the individual terms of the second order log-trace contribution are given in Appendix G.
Appendix F: Evaluation of the expectation values: M-fold In our expansion we include all the terms up to first order in 1 ξ 2 and up to sixth order in η.

Orientational contributions
As for the calculation of the purely orientational part we notice first that it is factorizing in the replica index α. It is thus possible to take the replica limit directly by means of the expansion . . . × q∈Z e − q 2 2κ |s1−s2| I 2 q (η) I 2 0 (η) e iq(∆ϕ1−∆ϕ2) .
The ∆ϕ-integrations lead to a Fourier transformation of the soft cross-link contributions and noticing that the resulting expression is symmetric in q we find The calculation of the expectation values in the second order contributions of the SIAS log-trace part can be done along the same lines as above using again the integral representation of the Bessel function (E.2). The calculation of the lowest order spatial contributions is the same as for the M-fold case. For the purely orientational we obtain f lt2,or = s1,s2 ϑ1,ϑ2 ln 1 + 2 ∞ q=1 I 2 q (η) I 2 0 (η) ×e − q 2 2κ |s1−s2| cos(q(ϑ 1 − ϑ 2 )) .
The lowest order coupling term is given by cos(ψ γ 1 −ϑ 1 +∆ϕ γ 1 )+cos(ψ γ 2 −ϑ 2 +∆ϕ γ 2 ) × cos(ψ δ 1 − ϑ 1 + ∆ϕ δ 1 ) + cos(ψ δ 2 − ϑ 2 + ∆ϕ δ 2 ) As in (F.3) the second term in the first line of the expression is already proportional to n and will thus finally lead to a contribution of O(n 2 ) that is not of interest for us. After integrating with respect to the ϑ-variables and then with respect to the ϕ-variables we obtain in the end −n η 2 8ξ 2 I 2 1 (γ) The results for the SIAS can be expressed in terms of the functions g, h and l that we already introduced in the preceding section.