Semiflexible Polymers in the Bulk and Confined by Planar Walls

Semiflexible polymers in solution under good solvent conditions can undergo an isotropic-nematic transition. This transition is reminiscent of the well-known entropically-driven transition of hard rods described by Onsager’s theory, but the flexibility of the macromolecules causes specific differences in behavior, such as anomalous long wavelength fluctuations in the ordered phase, which can be understood by the concept of the deflection length. A brief review of the recent progress in the understanding of these problems is given, summarizing results obtained by large-scale molecular dynamics simulations and density functional theory. These results include also the interaction of semiflexible polymers with hard walls and the wall-induced nematic order, which can give rise to capillary nematization in thin film geometry. Various earlier theoretical approaches to these problems are briefly mentioned, and an outlook on the status of experiments is given. It is argued that in many cases of interest, it is not possible to describe the scaled densities at the isotropic-nematic transition as functions of the ratio of the contour length and the persistence length alone, but the dependence on the ratio of chain diameter and persistence length also needs to be considered.


Introduction
There exist many macromolecules where local stiffness is relatively large, i.e., they have a persistence length ( p ) much larger than the diameter (d) of the effective monomeric units. Such semiflexible polymers may exhibit liquid-crystalline order [1,2], are of interest as building blocks of various complex soft materials [3,4] and also occur as ingredients of biological matter [5][6][7][8][9][10]. Here, we shall disregard all of these very interesting applications, focusing only on the generic case where these macromolecules exist in semidilute or concentrated solutions under good solvent conditions. However, applying a coarse-grained description, solvent molecules are not explicitly considered; monomeric units then exhibit an effective repulsive interaction, and the monomer concentration ρ, as well as the contour length L of the chains, are the basic parameters to control the properties of such systems. A characteristic feature, observed for large enough values of the ratio p /d, is the onset of nematic long range order [11][12][13][14][15][16][17][18][19][20][21][22]. The simplest case of such lyotropic liquid crystalline systems is solutions of hard rods, such as tobacco mosaic viruses in water [23], a problem that inspired Onsager [24] to develop his famous theory of purely entropically-driven phase transitions [25,26], resulting from the competition of translational and the orientational entropy contributions of the rods. The present problem reduces to this hard rod limit if L/d 1 and p is taken to infinity. It then suffices to consider the free energy of the system on the level of an approximation where only the second virial coefficient is kept [24]. Describing semiflexible polymers in the limit of the Kratky-Porod worm-like chain model [27], this Onsager-type treatment has been carried over to the description of the isotropic to nematic transition by Khokhlov and Semenov [12][13][14], Odijk [15,16], Chen [20] and others. In the isotropic phase, it is implied that the mean-squared end-to-end distance R 2 e of the chains is described by the Kratky-Porod model [27], which reduces to R 2 e = L 2 in the rod-limit (L p ), but implies Gaussian chain behavior ( R 2 e = 2 p L = l K L, l K being the Kuhn length [28]) in the opposite limit. Obviously, the swelling of the coils due to the excluded volume interaction between the monomers [29] is neglected, although excluded volume matters in the nematic phase. The monomer concentration ρ i , where the nematic order starts, is predicted to scale as follows [12][13][14][15][16]20]: Note that the isotropic to nematic transition is weakly of first order, but the width of the two-phase coexistence region (ρ i < ρ < ρ n , only for ρ ≥ ρ n , the system exhibits uniform nematic long-range order) is predicted to be rather narrow [12][13][14][15][16]20].
In many cases of practical interest, however, the system is not in the limit d/ p → 0, then, ρ i (and ρ n ) are not extremely small, and the assumption that only the second virial coefficient matters is not justified. Various attempts have been made to estimate the free energy going beyond the second virial coefficient [17][18][19]21], but unfortunately, many of these treatments require somewhat ad hoc assumptions and/or uncontrolled approximations. As a result, various approaches have led to different results not in agreement with each other, and hence, no general consensus on how to go beyond Equation (2) has emerged.
The present authors have taken a new approach towards these issues, combining large-scale molecular dynamics (MD) simulations [30,31] (feasible via the use of graphics processing units (GPUs) [32,33]) with classical density functional theory (DFT). The latter approach is considered in general the most powerful version of mean field theory to describe ordering phenomena in condensed matter and has seen broad and significant progress in recent years [34][35][36][37], but for the most part, this technique is concerned with simple liquids (described by point-like particles interacting via isotropic potentials). Both the generalization to anisotropic particles, whose interactions depend on their mutual orientations (see, e.g., [38]) and the generalization to flexible (e.g., [39]) and semiflexible polymers (e.g., [40,41]) are highly nontrivial and are still subjects of ongoing research. While for simple fluids, the basic object of the theory is the spatially nonuniform density ρ(r), for semiflexible macromolecules, one needs to operate with a function ρ mol (r, ω), which depends not only on the particle position r, but also on the local orientation ω (ω is a shorthand notation for the polar angles θ, φ of the molecular bonds). Thus, our recent work [42][43][44][45] not only goes far beyond previous simulation approaches (e.g., [46][47][48][49][50]) by simulating much larger systems (up to 700,000 monomers) and varying parameters, such as L and p over much wider ranges than were accessible in the earlier work, but we also extend the DFT methodology for semiflexible polymers by considering simultaneously both of their orientational distribution and spatial inhomogeneity. The comparison with MD, in turn, provides a stringent test of the conditions under which these extensions are accurate and elucidates the reasons for its limitations. A third new ingredient of the work [42][43][44][45] is the use of the concepts of the cylindrical confinement of semiflexible polymers [51][52][53][54] to interpret the anomalous fluctuations of the semiflexible chains, on the length scale of the deflection length λ along the contour. These long wavelength collective deflections of bundles of neighboring chains are soft modes that are not implicit in the DFT framework. At this point, we also note that these fluctuations are also not correctly included when one models the semiflexible chains by lattice models (which otherwise may be advantageous from the computational point of view [55][56][57][58][59][60][61][62]). The present off-lattice models thus have also a particular advantage when one wishes to test theories of the interaction of semiflexible polymers with confining walls (e.g., [63][64][65][66]).
The outline of the remainder of this review is as follows. In Section 2, we summarize the models on which the MD and DFT work is based, outline the basic aspects of the computational approaches that are used and remind the reader about the basic aspects of earlier work. In Section 3, we summarize the main features of the behavior of semiflexible polymers in the bulk, considering the variation of persistence length, contour length and monomer density, and also present some comparisons with earlier theories and experiments. Section 4 gives some results on the interaction of semiflexible chains with repulsive walls, including a discussion of capillary nematization in thin films. Section 5 gives a brief summary and an outlook.

Molecular Dynamics
There is much recent effort (e.g., [67][68][69][70]) in trying to construct simulation models of macromolecules that address specific effects of their chemical structure. This is not our focus here; we are rather concerned with generic models that address only the general features of lyotropic solutions of semiflexible macromolecules. Thus, the MD work is based on the standard bead-spring model [71,72] amended by a bond-angle potential (see, e.g., [73]), where consecutive effective monomers along the backbone of the chain (the beads) interact with the finitely-extensible nonlinear elastic (FENE) potential [71,72] U FENE (r) plus the repulsive part of the Lennard-Jones potential (Weeks-Chandler-Andersen [74] potential U WCA (r)). Explicitly, these potentials are defined as follows: while U FENE (r > r 0 ) ≡ 0, and: and U WCA (r > r c ) ≡ 0. The parameters of the FENE potential are taken as r 0 = 1.5σ and k = 30ε/σ 2 as usual, and then, the total binding potential U FENE (r) + U WCA (r) yields a rather sharp minimum at the effective bond length b ≈ 0.97σ. The contour length of a chain containing N beads is We stress that the potential U WCA (r) is also applied between any two effective monomers in the system, not only bonded ones. Note that the solvent molecules are not explicitly considered; hence, the model corresponds to very good solvent conditions, and the effective monomer diameter d can be taken as d = σ. Finally, the bond-bending potential is: where the energy parameter b controls the chain stiffness and the angle θ ijk is the bond angle between the two subsequent bond vectors, a i = r j − r i and a j = r k − r j . Choosing units such that σ = 1 sets the scale of length and ε = 1 sets the scale of energy (the thermal energy being chosen k B T = 1 throughout), the parameters chain length N and bending energy b , together with the monomer density ρ, are the control parameters of the model. We also note that b is simply related [75] to the persistence length p via the average of cos θ ijk : The definition in Equation (6) implies (via Equation (5)) that the persistence length plays the role of a coupling constant in the effective Hamiltonian of worm-like chains (continuous version of the Kratky-Porod model). Other definitions of the persistence length exist, e.g., [76], but have some clear disadvantages, e.g., they are not applicable in d = 2 dimensions [75,77].
Note that the approximate equalities in Equations (5) and (6) hold for b ≥ 2, and then, U bend (θ ijk ) is harmonic in θ ijk ; hence U bend (θ ijk ) = k B T/2 = 1/2, and hence, θ 2 ijk = 1/ b ; and thus, p / b = b (in our units). We stress that for our problem, the notion of persistence length makes sense only as a description of the local chain stiffness [75], and it must not be related to the asymptotic decay of bond-angle correlations via the textbook formula [28]: While Equation (7) reduces to Equation (6) for s = 1, it would fail for large s in general, apart from non-interacting phantom chains [28,75]. Due to excluded volume effects, the asymptotic decay of cos θ i,i+s is of power law form, as discussed extensively elsewhere [75], for p / b < s N, in the isotropic phase of the solution. In the nematic phase, we expect that cos θ i,i+s develops a plateau independent of s, related to the nematic order parameter.
MD simulations are done at constant density, choosing N chains of length N in a box of volume V box , and hence: In order to study the phase behavior in the bulk, one chooses a box of simple cubic shape (linear dimension L box = V 1/3 box ) and periodic boundary conditions in the x, y and z directions. When we are interested in the effect of repulsive walls, one chooses V box = L 2 x L z , with two repulsive walls at z = 0 and z = L z , and periodic boundary conditions only in x and y directions. The repulsive potential due to the walls then simply is: When one chooses L z very large, in particular L z L, for densities corresponding to the density of the bulk isotropic phase ρ b , the two walls can be treated as non-interacting, and thus, one can relate the behavior observed in such simulations to wall effects on semi-infinite solutions. In the case where the distance L z between the walls and the contour length L are comparable, this is not the case, and the problem of "capillary nematization" (well-known for hard-rod fluids; see, e.g., [78]) must be addressed.
In the MD simulations, configurations of the chains in the simulated volume are generated by solving numerically Newton's equation of motion for the beads, using the velocity Verlet algorithm [30,31], and applying a "Langevin thermostat" to fix k B T = 1: where m(= 1) is the mass of an effective monomer having position r n , t is the time along the generated system trajectory through phase space, F tot (r n ) describes the total force obtained as the gradient of the total potential, due to U FENE , U WCA , U bend (and U wall in the presence of repulsive walls) acting on the considered bead. The friction coefficient γ(= 0.25) is related to the random force F rand n (t) by the fluctuation-dissipation theorem: For the present choice of parameters, the MD time unit is τ MD = √ mσ 2 / = 1, and to actually solve Equation (10), discrete time increments δt = 0.01 are used, employing the HooMD-blue software package [32,33] on various GPUs. From the obtained trajectories, various quantities of interest, such as pressure, chain linear dimensions and orientational order parameters, can be deduced.

Density Functional Theory
Starting with the seminal work of Onsager [24], DFT has been the basic theoretical approach to deal with entropically-driven phase transitions. We choose here the nomenclature to reserve the notation DFT for the implementations developed by us recently [42][43][44][45] and refer to the earlier DFT theories by the names of the respective authors. For implementing DFT, it is more convenient to describe the polymer molecules as necklaces of tangent hard spheres of diameter σ, but one can choose the same bending potential as done in MD, Equation (5). Thus, it is also implied that non-bonded monomers interact with the simple hard-sphere potential; and when one includes in the model the interaction of monomers with repulsive walls, the latter are also taken to be hard walls: However, we first discuss the DFT implementation used for studying the phase behavior in the bulk [42,43]. The density ρ mol (r, ω) is then taken as the product of the chain density ρ mol = N /V and the orientational distribution function f (ω) [41,79,80], where f (ω) = 1/4π in the isotropic phase. The Helmholtz free energy then is decomposed, as usual, into ideal and excess terms: For a bulk system, the ideal term is, after the integration over ω: Obviously, the last term (the ideal orientational entropy) is zero in the isotropic phase, but nontrivial in the nematic phase.
When infinitely long rigid rods are considered, Onsager's theory [24] readily yields the second term in Equation (12), based on the second virial coefficient term in the virial expansion: Here, V exc (ω, ω ) is the excluded volume for two rods with orientations ω and ω . However, here, we wish to study semiflexible polymers of finite length rather than rigid rods of infinite length. To go beyond the second virial approximation, a common way is to introduce a rescaling prefactor [81,82] a resc : Egorov et al. [42,43] applied two different choices for a resc . The first one, leading to a version of DFT called DFT-CS (DFT-Carnahan-Starling), is based [81,82] on the Carnahan-Starling [83] equation of state for the simple hard-sphere fluid: Note that η is just the monomer packing fraction, and ρ b = Nρ mol . Equation (16) completely disregards chain connectivity, but has the merit that it reduces to the Onsager limit for η → 0.
The second choice [84] is based on an (approximate) form for the excess free energy F iso exc (ρ mol ) of the polymeric fluid in the isotropic state (but does not reproduce the Onsager limit for η → 0). While various expressions (see, e.g., [40,85,86]) are known for F iso exc (ρ mol ) for flexible polymers, F iso exc (ρ mol ) is not known for semiflexible polymers. Therefore, the generalized Flory dimer (GFD) equation of state [87] was used [42,43] to compute the rescaling factor a resc in this formulation [84], to obtain a version of DFT that was termed DFT-GFD.
However, both versions require the knowledge of V exc (ω, ω ) for semiflexible polymers, which is not known analytically, of course. Therefore, Egorov et al. [42,43] used an empirical expression proposed by Fynewever and Yethiraj [41] obtained by fitting the data from Monte Carlo simulations of a system containing just two semiflexible chains. Of course, from this description, it is evident that in spite of the rigorous foundation [37] of DFT in terms of a variational principle for the grand potential as a functional of the average particle density, the practical implementation of DFT in the present application is hampered by various more or less empirical assumptions and approximations whose accuracy it is difficult to assess a priori. The main motivation for the two versions of DFT that are presented here is that they are validated (at least qualitatively) over a fairly broad range of parameters by comparison with the MD results. We emphasize that in this treatment, we have effectively accounted for the effects of higher-order terms in the virial expansion, at least approximately.
Of course, when F(ρ) has been computed for the isotropic and nematic phases, the chemical potential and pressure for both phases can be easily obtained, and the phase diagram for the isotropic-nematic transition can be constructed. The order parameter in the nematic phase follows from the molecular orientational distribution function, recalling that ω is a shorthand notation for (θ, φ): When one now considers the extension of the theory to confinement by planar walls [44,45], the starting point is still Equation (12), but ρ mol (r, ω) can no longer be taken as ρ mol (r, ω) = ρ mol f (ω); rather, we must have: restricting attention to the case where the solution far from the walls is in the isotropic phase. We note here that while in the bulk DFT calculations (which do not resolve individual monomers on the chain), ω describes the orientation of the entire chain [42,43], the DFT calculations under planar confinement are performed on a monomer-resolved level, and f (z, ω) for the molecule is obtained by averaging the corresponding orientational distribution functions for the individual bonds over all of the bonds in the chain [44,45]. The task is now the minimization of the grand potential Ω: where µ is the polymer chemical potential and V mol ext (z, ω) is the external potential due to the two hard walls acting on the polymer molecules. The ideal term in Equation (12) for this case is still known exactly [44], while the excess term is split into an isotropic part F iso exc (ρ iso (z)) and an orientational part. As a generalization of Equation (14), one needs a model for the excluded volume term V exc (r, r , ω, ω ), which is both spatially and orientationally dependent. However, this term is known explicitly only for two rigid rods under planar confinement [88]; for two semiflexible polymers, we make a heuristic approximation using again the same approximation for V exc (ω, ω ) as used for the bulk [41]. We note here that a better approximation would have been obtained by averaging the corresponding Mayer function over x and y variables. However, we do not pursue this approach here, because even within the crude approximation of Equation (20), the minimization of Equation (19) with respect to the two functions ρ iso (z) and f (z, ω) is delicate and requires substantial numerical effort (see [44], for details). The generalization of Equation (17) then is: to obtain the order parameter as a function of distance from the walls. Here, θ is understood as a polar angle with the z-axis perpendicular to the walls (unlike Equation (17), where nematic order is assumed and θ is measured relative to the director), since we focus here on the orientational order induced by the walls. Thus, S(z) = 0 corresponds to random chain orientation (recall that the orientational distribution function f (z, ω) is defined as an average over all of the bonds in the molecule), while S(z) = −0.5 corresponds to perfect alignment of the chain parallel to the wall. Note that a particular bonus of DFT is that it yields the free energy explicitly, which will depend on L z due to the wall excess contributions. The dimensionless surface tension hence can be obtained when the corresponding bulk term L z f bulk is subtracted ( f bulk is the bulk free energy density, and it is assumed that F is normalized per unit area of the walls): where Ω and Ω bulk are the grand potential density and its bulk value, respectively.

A Brief Review of Earlier Theories
Onsager's theory for the isotropic-nematic transition of rigid rods [24] was extended to semiflexible polymers by Khokhlov and Semenov [12][13][14] and Odijk [15,16]. Some approximations made by these authors are avoided in the treatment of Chen [20], and hence, we only sketch this latter treatment here briefly. The free energy per chain of the system is written as a sum of three terms (c = N L 2 d/V box is a dimensionless number density): Here, Q is the partition function of a semiflexible chain, and hence, the first two terms on the right-hand side represent the entropy of a semiflexible polymer. The last term represents the excluded volume interaction between the two chains; α is the angle between two unit vectors pointing at ω and ω . The mean field U MF (ω) represents the averaged orienting effect of the neighboring chains on the considered chain and needs to be determined self-consistently. Thus, from Equation (23), it is clear that excluded volume is only dealt with on the level of the second virial coefficient, and fluctuations beyond the mean-field approximation are neglected.
The explicit computation of F is based on the Kratky-Porod [27] model for semiflexible chains, using a functional integral approach [89], where the semiflexible polymer is described by a continuous space curve, specified by its tangent unit vector n(t), 0 ≤ t ≤ 1 (t =0, 1 correspond to the free ends of the polymer chain). The statistical probability P{n(t)} of such a chain configuration is [89]: There is no excluded volume considered in Equation (24), and thus, Equation (24) can be shown to yield Equation (1). One introduces a partition function q(t, ω) of a chain of length tL that has the final end-point pointing at orientation ω, which satisfies the equation [90]: with q(t = 0, ω) = 1, and is related to f (ω) via: where Q = dωq(1, ω). Of course, this set of coupled self-consistent nonlinear equations cannot be solved analytically, but requires a numerical iteration procedure [20]. Khokhlov and Semenov [12][13][14] and Odijk [15,16] have avoided this problem by using an approximate variational method to minimize F, employing trial functions f (ω) with a single variational parameter. We shall briefly discuss the differences between these results when we compare them to the results of DFT and MD methods [42,43]. Here, we also mention the extension of the Khokhlov-Semenov theory due to Hentschke [17], who attempted to improve the treatment of the orientational free energy of semiflexible polymers by accounting more carefully for the excluded volume effects, while DuPré and Yang [19] attempted to go beyond the second virial coefficient in their treatment of the excluded volume interaction, in order to enable the study of liquid-like densities. A related goal was addressed by Sato and Teramoto [18,21], who extended the scaled particle theory [91,92] to renormalize the strength of the excluded volume interaction in Equation (23), modifying also the ideal gas-like term. With respect to f (ω), Sato and Teramoto [18,21] continued to use trial functions similar as done by Khokhlov and Semenov [12][13][14]. Unfortunately, the accuracy of these various extensions of the Khokhlov-Semenov theory is hard to assess a priori, and hence, this problem is outside of the scope of the present brief review.

The Isotropic-Nematic Transition and Its Dependence on p , L, and d
We discuss here the variation of system properties as a function of the density ρ of the effective monomers, Equation (8). Since the isotropic-nematic transition is of the first order, one expects to find a two-phase coexistence region from ρ = ρ i to ρ = ρ n , at a coexistence pressure p = p coex . The order parameter S (Equation (17)) should vary linearly from S(ρ i ) = 0 to S(ρ n ) = S c , due to the lever rule of two-phase coexistence. While this behavior is consistent with the DFT calculations, where one treats the isotropic and nematic phases separately and locates the transition from the condition that both the chemical potentials of coexisting phases and their pressures must be equal, MD work in the canonical constant density ensemble is hampered by finite size effects, and so, in the pressure vs. density isotherms (Figure 1a,b), the transition only shows up as a small wiggle, rather than as a strictly horizontal line. Furthermore, in the variation of the order parameter S with density finite size, effects cause a "finite size tail" [93] in the disordered phase, a well-known effect from simulation studies of other phase transitions. Therefore, the kink singularities in the S(ρ) vs. ρ curve at ρ = ρ i and ρ = ρ n are not seen due to the finite size rounding. An interesting observation is the fact that the order parameter predicted by DFT reaches saturation much faster than found in MD (Figure 1c,d). In this context, it is interesting to note that for the nematic order parameter of rod-like molecules described by the hard Gaussian overlap fluid, also an overestimation of the order parameter by DFT in comparison with simulation results has been observed [94]. We shall discuss the behavior of the order parameter and its relation to the deflection length in Section 3.2 below.
Testing now the quantitative accuracy of the two versions of DFT introduced in Section 2.2 for the equation of state (Figure 1a,b) by comparing the p vs. ρ isotherms to the corresponding MD data, a clear trend as to which version works better does not emerge. For large values of b , the width of the I-N (Isotropic-Nematic) coexistence region is so small, that on the scale of Figure 1b, it is hardly resolved. Of course, a perfect quantitative agreement between MD and DFT should not even be expected, since the underlying chain models differ slightly. From Figure 1c,d, we see that the theory of Chen [20] predicts the I-N coexistence region rather well as long as the transition densities are small enough (ρ < 0.3), but becomes more and more inaccurate the larger the transition densities get. Such a failure is expected, of course, since the theories [12][13][14][15][16]20] describe the chain interactions only on the level of the second virial coefficient, as discussed in Section 2.2.  The theory of Chen [20] only gives information on the variation of the transition densities ρ i , ρ n as a function of the parameters p and L and does not discuss the properties of the system in the nematically ordered phase further. Thus, we proceed next to a comparison of these predictions [20] to the MD and DFT results [42]; Figure 2. It is seen that for very stiff chains, both the theory of Chen [20] and DFT [42] are in reasonable agreement with the MD results, and nematic order still occurs for chains as short as N = 8. For relatively flexible chains ( b = 8), DFT and MD agree rather well, while the theory based on the second virial approximation [20] is far off, as expected, because the transition densities are rather large. At the same time, we see that DFT-CS predicts a spurious upturn of the transition densities with N for intermediate values of b (specifically, b = 16 and 32), which is not confirmed by MD [42]. As stated already in Equation (2) and can be deduced from Equation (23), noting that c ≈ ρ if the distinction between L/d and N is neglected (strictly speaking, L = (N − 1) b , and we take b ≈ d ≈ σ here), the scaled transition densities ρ i p /d, ρ n p /d are functions of the ratio p /L only. Hence, when one studies the variation of these transition densities with the parameter d/ p keeping the ratio p /L fixed, one expects simply to find horizontal straight lines. However, Figure 3 shows that this is not the case: the scaled transition densities distinctly decrease with increasing d/ p . This behavior is rather well described by DFT-CS for N/ p = 0.25, 0.5 and 1, while for N/ p = 2, the agreement is only qualitative. Replotting this latter case choosing a logarithmic rather than linear abscissa scale, we have also included predictions from the theories of DuPré and Yang [19] and Sato and Teramoto [18,21] for comparison. As has been discussed in Section 2.3, these theories renormalize the prefactor of the second virial term in Equation (23), to allow an extension of the description to higher densities and, thus, predict now a dependence of the transition densities on the parameter d/ p even when L/ p is constant. Figure 3b shows, however, that the variation predicted by the scaled particle theory [18,21] is far too strong, at least for the shown example. In this case, the theory of DuPré and Yang [19] is rather close to both DFT and MD results [42]. It is also shown that the differences between the theories of Khokhlov and Semenov [12][13][14], Odijk [15,16] and Chen [20] do not seem important, in comparison to the shortcoming that they do not yield any dependence on the parameter d/ p at all. There is also some disagreement between the two versions of DFT introduced here (Section 2.2); in view of the rather crude approximations that were necessary to introduce in the DFT framework, these discrepancies are not really surprising. We also mention that many of the corresponding experimental data (see [22] for a review) are believed to correspond to the regime where 0.005 < d/ p < 0.03. In this range, the deviations between the MD results [42] and Chen's theory [20] are typically less that 15%. Since experimental data are hampered by polydispersity and by significant uncertainty regarding both parameters d and p [22], very good quantitative agreement between theory and experiment cannot be expected. However, in those cases, when d/ p is not so small, the dependence of the transition densities on this parameter found in [42,43] should be relevant.  [43]; while the horizontal shaded stripes show the I-N coexistence regions from [20]. Note that ρL/d = (ρ p /d)L/ p is plotted here rather than ρ p /d discussed in the text, to avoid cluttering the figure, and the factor π/4 accounts for the normalization of the density with the cylinder volume p d 2 π/4 as in the Onsager theory; (b) Transition density ρ i π p /(4d) plotted vs −1 p , comparing MD data to the theories of Chen [20], Odijk [16], Khokhlov and Semenov [12,13], scaled particle theory (SPT) [18,21], DFT-CS, DFT-generalized Flory dimer (GFD) and the DuPré-Yang theory [19], as indicated in the key. Reproduced from [43] with permission from the Royal Society of Chemistry.

Nematic Order Described as an Effective Cylindrical Confinement
Here, we return to Figure 1c, where the nematic order parameter was plotted as a function of density. A remarkable feature is that the increase of S(ρ) from S c towards the ultimate saturation value S = 1 predicted by DFT is much more rapid than according to MD, irrespective of how well the transition densities (ρ i , ρ n ) are predicted. A clue to this qualitative discrepancy is obtained from an examination of snapshot pictures of the chain configurations in the nematic phase ( Figure 4). Even when the nematic order parameter is already large, the chains still exhibit considerable orientational disorder on large length scales. Superimposed on these large wavelength deflections, there are also some more or less random small-scale orientational fluctuations of the bond orientations relative to such a coarse-grained contour, which can be taken schematically as the axis of the bent tube of diameter 2r ρ (see Equation (27)) shown doubly shaded in Figure 4c. Since the excluded volume interactions between monomers from different chains are strictly respected, we can define these bent tubes such that they contain monomers from the considered chain only. In the coordinate system along the tube contour, we then have (for N 1): Since for large S, the deflections of the tubes away from the z-axis (which we orient along the nematic director) can involve only small polar angles θ, we shall have r eff ≈ λθ, where the deflection length λ is a characteristic length of the problem. (c) Schematic description of nematic order as effective cylindrical confinement: each chain has its own cylindrical bent tube of diameter 2r ρ defined such that it contains only monomers from the considered chain. This tube roughly follows the contour of this macromolecule, which shows long wavelength undulations with a typical wavelength given by the deflection length. The typical amplitude of these deflections is of the order r eff defining a cylinder (the straight axis of this cylinder is oriented along the director of the nematic phase). This cylinder contains not only a single bent tube, but rather is densely filled by a whole bundle of neighboring tubes whose deflections are strongly correlated. Reproduced from [42] with permission from the APS.
The deflection length λ is a concept well known for the problem of confinement of a single semiflexible chain in a cylinder with repulsive walls [51][52][53][54]. The deflection length concept means that the orientation correlation function cos θ(s) along a chain reaches a plateau as a function of s, because a chain confined in a cylinder is "deflected back" when it reaches the cylinder walls, and then, the angular mean-square displacement θ 2 (s) cannot increase any further. If one considers a regime where hairpin formation [95,96] can still be neglected, one can still apply Equation (7) for distances s b < λ, and we use cos θ(s) ≈ 1 − 1 2 θ 2 (s) to conclude: and putting, s b = λ one finds λ = ( p r 2 eff /2) 1/3 . Since S = 3 2 cos 2 θ − 1 2 ≈ 1 − 3 θ 2 /2 = 1 − 3/2(r eff /λ) 2 = 1 − (3/2 1/3 )(r eff / p ) 2/3 , one concludes that 1 − S ≈ (3/2)(2r eff / p ) 2/3 . By a similar argument, one can estimate the z-component of the end-to-end vector of the (strongly stretched) chain in the cylindrical tube. One can consider the macromolecule as a sequence of essentially straight pieces of length λ with an average misorientation given by the factor cos θ(s) ≈ 1 − 1 2 θ 2 (s) . Roughly, the mean-squared end-to-end distance R 2 e will be reduced relative to L 2 by a corresponding factor, putting R 2 e ≈ |R ez | 2 ≈ L 2 (1 − θ 2 ). Of course, such an order of magnitude estimates can be substantiated by more accurate theories [52][53][54] to yield: Equation (29) was derived for a single semiflexible chain confined in a cylinder of radius r eff with repulsive walls. The key idea of Odijk [15,16] and Egorov et al. [42,43] has been to postulate that the effect of the "nematic mean field" orienting the considered semiflexible macromolecule in the nematic phase can be described by a confinement in an effective cylinder. Thus, Equation (28) for s b = λ implies a relation between the deflection length and the order parameter reduction, namely: Odijk [16] also suggested that the concept of the deflection length implies chain-end effects on the local orientational order S(s) along the chain: near s ≈ 0 or s ≈ L, with S ∞ being the local order far from the chain ends, for L p . Equations (29)-(31) were first tested by MD simulations by Egorov et al. [42,43], and we reproduce their key results in Figures 5 and 6. Indeed, the description developed above is confirmed, at least qualitatively. Equation (31) is compatible with the data, and the resulting estimates of λ are roughly proportional to p (but somewhat smaller than predicted by Equation (30)). Figure 5b shows that r eff exceeds r ρ by far, for small densities, confirming the qualitative picture ( Figure 4c); when S decreases (with decreasing density), the deflection length increases, and hence, also r eff increases.
The inset of Figure 5b shows that even for N = 128, the chains are too short to clearly display a pronounced horizontal plateau in the plot of r 2 ⊥ (i) vs. i in the middle part of the chains, and hence, the accuracy with which estimates for r eff and λ can be extracted here is still limited.  Figure 4c) is included for comparison. The inset shows a plot of the transverse mean-square displacements r 2 ⊥ (i) , relative to the end-to-end vector of the chain, as a function of i, for b = 128 at two densities, as indicated. Equation (28) implies that r 2 ⊥ (i) increases linearly with i and saturates at i ≈ λ with r 2 ⊥ (i) ≈ (r eff /λ) 2 . The data points for r eff are extracted from these maximum transverse displacements, and the estimates extracted from Equation (29) are included for comparison. Reproduced from [43] with permission from the Royal Society of Chemistry.  (29). Reproduced from [43] with permission from the Royal Society of Chemistry.
A particularly interesting comparison is seen in Figure 6, which shows that for 1 − S 1, Equation (29) is quantitatively satisfied, for all N and b values shown there, without any adjustable parameters. This implies that in the well-ordered nematic state, the picture of the ordering as effective cylindrical confinement is self-consistent. It is clear that for 1 − S ≥ 0.3, this picture gradually breaks down; we have used θ 2 1 throughout, and this no longer holds then. Therefore, the description in terms of cylindrical confinement is no longer accurate in the nematic phase near the I-N transition (where 1 − S ≈ 0.5).
We conclude this picture by stressing that Figure 4c (which emphasizes r eff r ρ , as verified in Figure 5b) implies that there must occur collective deflection fluctuations of many neighboring chains, since the arrangement of the bent tubes with the diameter 2r ρ is space filling and the tubes must not overlap each other. Such collective modes (with rather large wavelengths) are missing in the DFT descriptions of the ordered phase, of course: the picture is analogous to the molecular field theory of (classical) isotropic magnets, where the Langevin function also predicts a much faster saturation of order, rather than theories taking into account the long wavelength spin waves.
An intriguing question is to connect this description to the long wavelength description of fluctuations in the nematic in terms of the Frank elastic constants. Gemünden and Daoulas [97] succeeded in estimating the latter for a discretized worm-like chain model where bonds interact with a soft anisotropic potential.
It would also be very desirable if one would have experimental data on these issues to compare. X-ray scattering experiments from polymer nematic liquid crystals (such as poly-γ-benzyl glutamate in dioxane) have been performed and reveal very interesting information on the anisotropic character of density fluctuations [98], but cannot elucidate the nontrivial interplay with the fluctuations of orientational order.

Semiflexible Chains Confined by Repulsive Walls
The effect of repulsive walls on solutions of semiflexible polymers is somewhat subtle: the chains have reduced translational freedom near a wall, but their orientation parallel to a wall gets enhanced. The first effect dominates as long as the solution is relatively dilute, while the second effect will lead to local nematic order near the wall when enhancement of density enforces enough chains to be located close to the walls.
These qualitative expectations are substantiated by the MD and DFT calculations of Egorov et al. [44]. Figure 7 gives an example for relatively short chains (N = 32) at ρ = 0.10, where a choice of L z = 40 is enough to ensure that the system at distances near z = L z /2 still exhibits bulk-like behavior. In this case, the MD simulations suffer from the obvious problem that MD is performed at a constant monomer density ρ, which is chosen beforehand, but in general differs from the bulk density ρ b (which is seen eventually in the middle of the slit, provided L z is large enough) by a correction of order 1/L z , due to the wall excess density. Of course, the choice of L z is somewhat arbitrary, and there clearly is interest in data that are not affected by a dependence on such a parameter (which drops out in the limit of a semi-infinite system, but the latter is not accessible by simulations). For a meaningful quantitative comparison between MD and DFT (Figure 7a), where the density ρ b corresponds to that of a bulk system in the grand canonical ensemble, Egorov et al. [44] chose the chemical potential of the DFT calculation, such that ρ b coincides with ρ middle as observed in the MD simulation. One can see from Figure 7a that ρ b exceeds ρ by about 14% due to the negative surface excess density at the repulsive walls. However, when this readjustment is made, nearly perfect agreement between DFT and MD is noted. We also mention that MD cannot use the definition given by Equation (22) to estimate the surface tension, but rather one exploits the anisotropy of the pressure tensor (e.g., [99]) p αβ (z) [α, β = x, y, z]: A nontrivial issue is the fact that the implementation of Equation (32) cannot follow the standard prescriptions [100] since p αβ (z) near the walls is strongly affected by the three-body forces deriving from the bending potential given by Equation (5) [101] even though these forces cancel out in the bulk. We also note (Figure 7a) that the bending potential, that distinguishes the semiflexible polymers from the flexible ones, leads to a much larger range of z near the walls over which the monomer density is depleted. Note also that for ρ = 0.1, there is no indication yet of the oscillatory density profile near the wall ("layering") that occurs at higher densities.  For simplicity, in Figure 7b,c, the distinction between ρ b and ρ is neglected, and since the chain models of MD and DFT differ slightly, it would be unrealistic to expect perfect agreement between the two methods (such as noted in Figure 7a) in general. However, both MD and DFT show the same qualitative trends, namely an increase of γ with either b or N. Note that the observed agreement between MD and DFT can only be achieved provided that the DFT takes into account both the spatial density variations (contributing already to the isotropic part of γ) and the orientational effects ( Figure 7c).
The bonus of the DFT calculations is that very precise results can also be obtained for the density profiles of individual monomers, such as the densities of chain-ends or mid-monomers, respectively ( Figure 8). These data are again in fair agreement with the corresponding MD results [44] (not shown here), but the latter suffer from strong statistical scatter.
For the low densities shown in Figure 8, there is not yet any trace of the familiar density oscillations ("layering") near the repulsive walls that is found at distinctly higher densities, rather there is a clear depression of the density near the walls, and this depression extends more and more towards the film center when the chain stiffness increases. Thus, for low densities, indeed, there is no strong tendency in favor of local nematic order near the walls, but instead, there is a significant depletion of the density ρ mid of the middle monomers of the chains near the walls over a range that increases with stiffness, followed by a maximum in their density profile further away from the walls.  This behavior changes when we study densities close to the transition density ρ tr of the bulk (Figure 9a). While S(z) (Equation (21)) in the isotropic phase far from the transition decays to zero rather fast, this decay gradually becomes slower as ρ approaches ρ tr , and very close to it, the decay occurs in two steps, indicating the formation of nematically-ordered layers attached to the walls. Right at ρ = ρ tr , the system is already well-ordered throughout the film (this "capillary nematization" [78] effect will be discussed below in more detail). This surface-induced ordering can be interpreted as the unbinding of an interface [102] between the nematic phase and the isotropic phase from the wall, and for short-range forces, between the wall and the interface (which are implied here because of Equation (9)), one predicts a logarithmic growth of the thickness of the nematically ordered layer. This thickness can be measured by the surface excess order parameter due to a wall, as long as S(z ≈ L z /2) = 0, so one still has the isotropic phase in the middle of the film. However, the isotropic-nematic interface is a mesoscopic, slowly fluctuating object, and hence, the statistical fluctuations of Ψ s when sampled from MD [44] or from Monte Carlo (MC) simulations [60][61][62] are huge. Both the MD data shown in Figure 9b and the earlier MC work [60,61] are compatible with the theoretically-expected variation [102]: but more precise data clearly would be desirable. The MC work [60][61][62] was based on the bond fluctuation model [103] for polymers on the simple cubic lattice, and this work suffers from the additional problem that the orientation of the director in the nematic layer can be only along the x or y axes of the lattice. dzS(z), vs. ρ tr − ρ, for the case N = 16, where ρ tr = 0.30, using data for L z = 40 and L z = 100 (to check for finite size effects) and displaying data for Ψ s extracted from the range from z = 0 to L z /2, as well as from z = L z /2 to L z , to illustrate the large statistical scatter. The straight line illustrates the fit to a logarithmic variation, |Ψ s | = 0.36 − 0.79 ln(ρ tr − ρ). Reproduced from [43] with the permission of AIP Publishing.
The order parameter S(z) defined in Equation (21), where θ is the polar angle with respect to the z-axis, which was used in Figure 9, measures only the extent to which the bonds are oriented parallel to the wall. It does not measure the extent to which the bonds are aligned with a director. Using thus the local tensor Q α,β n,i characterizing the orientation of a unit vector u n,i along the bond r i+1,n − r i,n connecting monomers i, i + 1 of the n-th chain, one finds a local director in a slice of width δz around z and averages Q α,β n,i only over all of the bonds for which r i,cm = (r i+1,n + r i,n )/2 falls inside the slice. The largest eigenvalue λ + (z) then characterizes the proper local order (and the associated eigenvector is the local director). Figure 10 shows an example of the local directors in typical configurations of the system at three densities. One can see that at the lowest density ρ = 0.1 for this case of rather stiff chains ( p = 32, N = 32), whose end-to-end distance in the bulk is R 2 e 1/2 = 25.82, and hence, of the same order is L z , the arrows are more or less parallel to the walls, but in the xy-plane, their orientation is still rather random: although S(z) as defined from Equation (17) is close to −1/2 throughout the film, this clearly is not a nematically-ordered state. For ρ = 0.2, on the other hand, there are several layers close to both walls where the directors are oriented parallel to each other, and only in the center of the film there still occurs a misalignment, as a remainder of a phase, which is still isotropic in d = 2 dimensions. For ρ = 0.3 (which corresponds to ρ tr in the bulk), there is already a very high degree of order. This picture is substantiated when we record a thermally-averaged order parameter profile λ + (z); Figure 11. While for ρ ≤ 0.2, there is only rather little order at the walls, for ρ = 0.25, both walls are clearly coated by ordered layers. For ρ = 0.27, the order parameter in the bulk is already about 0.5: capillary nematization has occurred. Near the bulk transition (ρ = 0.3), the order in the film is clearly larger than it would be in the corresponding bulk system (Figure 1c). . Thermally-averaged order parameter profile λ + (z) vs. distance z for the case N = 32, b = 32, L z = 40, N = 1500 and various densities, as indicated. Note that in the bulk, the I-N transition occurs at ρ tr ≈ 0.30. Reproduced from [45] with the permission of Wiley-VCH.
We add a caveat: the observation of quasi-two-dimensional long-range nematic order implied by this discussion of Figure 11 could be just a finite size effect. As is well-known, the existence of nematic order in d = 2 dimensions is a controversial issue (analogous to the Kosterlitz-Thouless transition [104] of two-dimensional XY-ferromagnets, there could be a state with a power law decay of orientational correlation functions, rather than true long-range order; see, e.g., [105,106]). In contrast, the lattice Monte Carlo simulations of Ivanov et al. [60][61][62] did find a well-defined order-disorder transition of the quasi-two-dimensional wall-attached layers, but due to the only two discrete director orientations, this is an Ising model-like transition and not a faithful representation of the nematic order of real semiflexible polymers.
Finally, Figure 12 shows corresponding DFT results. One can clearly see that for thin films, S is never strictly zero, due to the wall-induced order. For small values of L z , such as L z = 30 or L z = 40, the variation of S with µ is clearly nonsingular, increasing µ (or increasing the density) simply causes a gradual onset of order, without a sharp transition. Egorov et al. [45] concluded that near L z = L * z ≈ 55, a capillary critical point occurs, and for L z > L * z , there is a first-order transition (from a state with a smaller nematic order parameter to a state with a larger order parameter). Of course, due to its mean-field nature, DFT cannot answer the questions about the character of two-dimensional long-range nematic order.  It is encouraging that first experimental results on capillary nematization of colloidal rods under confinement have very recently become available [107], and we hope that the present article will stimulate corresponding work on confined semiflexible polymers.

Summary
In this review, recent theoretical and simulation work on the isotropic-nematic transition and the nematic order of semiflexible polymers was reviewed, considering both MD simulations and DFT calculations, but earlier theoretical work by Khokhlov and Semenov, Odijk, Chen and others [12][13][14][15][16][17][18][19][20][21], was also briefly included in the discussion. Only semiflexible polymers in lyotropic solutions were discussed, where the effective interactions between the monomeric units are short-ranged and repulsive in character, representing effectively the excluded volume of these units. Solvent molecules are not explicitly considered, and as in previous theoretical work [12][13][14][15][16][17][18][19][20][21], no attention was payed to the detailed atomistic structure and the corresponding potentials (e.g., torsional potentials etc.); and also, electrostatic forces were disregarded throughout. Thus, the stiffness of the semiflexible polymers was dealt with on a coarse-grained level, namely via the local persistence length p (or the corresponding energy parameter b of the bending potential between subsequent effective bonds along the chain, cf. Equations (5) and (6)). The repulsive monomer-monomer interaction is characterized by another characteristic length, the effective diameter d of the resulting worm-like chain, or (equivalently) the diameter σ = d of the effective monomeric units, which were modeled via the potential U WCA (r) (see Equation (4)) in the MD framework, or by hard spheres in the DFT framework. Chain connectivity was modeled by anharmonic springs (described by U FENE (r) + U WCA (r)) in the MD work and by requiring the hard spheres to be tangent in the DFT work, while much of the earlier theories [12][13][14][15][16][17][18][19][20][21] were based on the continuum version (Equation (24)) of the Kratky-Porod model, where chain interactions then are described like in Onsager's theory for long and thin hard rods [24] via the second virial coefficient [12][13][14][15][16]20] or modifications thereof [17][18][19]21]. We have not reviewed here the early simulation work (e.g., [46][47][48][49][50]), since most of this work dealt with comparatively short chains and relatively small systems, less suitable to characterize the I-N transition and the character of the nematic phase, but we have included results from early work (e.g., [49]) where appropriate. Both the behavior in the bulk solution and the effect of confinement by repulsive planar walls was considered.
It was shown that the scaled density where the I-N transition takes place shows a distinct dependence on both dimensionless parameters L/ p and d/ p , and not only on L/ p alone, as the theories of Khokhlov and Semenov, Odijk and Chen [12][13][14][15][16]20] imply. DFT can account for this additional dependence on d/ p rather well, for typical cases that were studied, while the scaled particle theory [18,21] strongly overestimates this dependence on d/ p . As expected, for d/ p less than 0.01, this dependence on d/ p becomes weak, and the approximations based on the second virial truncation [12][13][14][15][16]20] then become reasonably accurate. On the other hand, DFT predicts a spurious upturn of the transition density studied as a function of L for intermediate values of p , which is not confirmed by the simulations. Particularly drastic approximations are needed when one develops DFT to study the behavior of semiflexible chains near hard walls: it is necessary to neglect the coupling between the spatially inhomogeneous density profile and the orientational interaction term (Equation (20)), which needs to be inferred from the bulk (and there it needs to be taken from dedicated two-chain MC simulations [41]). Obviously, the DFT formulation for semiflexible polymers cannot yet be cast in the form of a completely self-contained analytical theory, it is still based on somewhat heuristic approximations. Nevertheless, the comparison with the MD results shows that particularly at high densities, it performs significantly better than the other theories. However, in the nematically-ordered phase, it predicts too large values of the nematic order parameter. This problem can be attributed to the neglect of collective deflection modes of the chains relative to the director. These deflections can be understood in terms of the analogy between semiflexible chains confined in cylindrical tubes and in the nematic phase.
A key point of our description is (Figure 5b) that the deflection length λ is not related to the perpendicular monomer displacements of order r ρ (nearest monomer distance in the plane perpendicular to the director), but a much larger length r eff due to collective chain bending, contrary to naive expectations [16,98].
When the effect of repulsive walls is considered, one finds that for small monomer densities, the conformations of chains located near the wall are strongly deformed, leading to an enhancement of chain-end densities near the walls, while the middle monomers are depleted near the wall and are more likely to be found away from it. At larger densities, one finds that a wall-induced nematic order sets in. Of course, only a rough estimate could be given of the capillary nematization critical point, and a more complete variation of all of the parameters ( p , L) for this problem is still lacking. No attempt to clarify the nature of the quasi-two-dimensional phase in thin film geometry could be made.
Egorov et al. [42,43] also have suggested that the results reviewed here may also help to better understand the available experimental data on the I-N transition of semiflexible polymers, if one takes the predicted dependence on the parameter d/ p into account. Unfortunately, these parameters d and p are known only rather imprecisely, and hence, it is difficult to draw firm conclusions on the experiments. Thus, it remains a challenge to improve the database on which a really conclusive comparison between experiment and theory can be based. Future work could also consider the effect of confinement by planar walls with attractive interactions, spherical confining surfaces, etc. Finally, we note that the focus of the present review is on the isotropic-nematic phase behavior of semiflexible polymers. In this regard, it is important to mention that the effect of stiffness on the nematic-smectic phase transition of semiflexible polymers is another interesting problem, which has been actively studied both experimentally and theoretically [108][109][110][111]. Thus, calculating more complete phase diagrams by including also the smectic phase would be of interest and could be considered in future research.