Tube Concept for Entangled Stiff Fibers Predicts Their Dynamics in Space and Time
We study dynamically crowded solutions of stiff fibers deep in the semidilute regime, where the motion of a single constituent becomes increasingly confined to a narrow tube. The spatiotemporal dynamics for wave numbers resolving the motion in the confining tube becomes accessible in Brownian dynamics simulations upon employing a geometry-adapted neighbor list. We demonstrate that in such crowded environments the intermediate scattering function, characterizing the motion in space and time, can be predicted quantitatively by simulating a single freely diffusing phantom needle only, yet with very unusual diffusion coefficients.
pacs:87.15.hj, 87.15.H-, 66.10.C-
Long, stiff filaments are abundant in nature and nanotechnology and form entangled meshworks of remarkable mechanical response and complex dynamic behavior Bausch and Kroy (2006). Examples include biopolymers such as filamentous actin (F-actin) Liu et al. (2006); Koenderink et al. (2006); Wong et al. (2004), microtubuli Lin et al. (2007), fd viruses Lettinga and Grelet (2007); Addas et al. (2004), and xanthan Koenderink et al. (2004), as well as more recently synthetically fabricated carbon nanotubes, cellulose whiskers, and polymeric stiff rods Cassagnau et al. (2013). In the semidilute regime the excluded volume becomes irrelevant, such that the filaments in solution diffuse as infinitely thin needles and their peculiar transport properties emerge due to topological constraints imposed by their neighboring impenetrable filaments. As a consequence, the motion of a single constituent in the highly entangled regime is suppressed to a sliding back-and-forth movement in an effective tube composed of the surrounding filaments. This tube concept, pioneered by Doi and Edwards Doi (1975); Doi and Edwards (1978), thus reduces the complex many-body dynamics to an effective single-particle motion. Most prominently, it predicts a drastic suppression of the rotational diffusion coefficient Doi (1975); Doi et al. (1984); Tao et al. (2006); Tse and Andersen (2013) as well as of the translational diffusion coefficient perpendicular to the orientation of the stiff filament Teraoka and Hayakawa (1988); Szamel (1993).
Rod and needle systems have been studied in simulations and theory extensively Löwen (1994); Tucker and Hernandez (2010); Sussman and Schweizer (2011); Yamamoto and Schweizer (2015), differing in the underlying dynamics and aspect ratio of the constituents. For Newtonian dynamics, the transport coefficients behave rather differently Frenkel and Maguire (1981, 1983); Otto et al. (2006); Höfling et al. (2008a); Tucker and Hernandez (2011), characterized by an increase in the center of mass diffusion coefficient, first observed by Frenkel and Maguire Frenkel and Maguire (1981). A long-standing debate how a finite aspect ratio or a small flexibility of the stiff Brownian rods affects the Doi-Edwards scaling Doi et al. (1984); Fixman (1985), has been resolved only recently Cobb and Butler (2005); Tao et al. (2006); Tse and Andersen (2013); Zhao and Wang (2013), pointing out, that the densities considered so far just reach the onset of the asymptotic regime and corrections to scaling are relevant Tse and Andersen (2013). Within two-dimensional toy models Moreno and Kob (2004); Höfling et al. (2008b); Munk et al. (2009), higher densities can be achieved in the simulation and the Doi-Edwards picture has been shown to remain valid. Yet, extracting spatiotemporal information, in principle accessible in light-scattering experiments Doi and Edwards (1978); Aragón S. and Pecora (1985) or simulations, has not been an easy task and the fundamental implications for the motion in the tube have remained unexplored.
Here, we show that the tube concept can be elaborated to quantitatively predict the intermediate scattering function of highly entangled suspensions of stiff filaments at all time and length scales. We also compare to a needle Lorentz system, where a single tracer needle explores a quenched disordered array of needles, i.e. the proverbial needle in a haystack, and demonstrate that the dynamic rearrangement of the surrounding needles induces no qualitative changes. We corroborate the scaling behavior of the transport coefficients and show that the spatiotemporal dynamics can be rationalized in terms of a single phantom needle where transport is dominated by translation-rotation coupling. For very high entanglement analytic progress is made and we show that the intermediate scattering function can be evaluated in closed form for long times.
Transport coefficients.— In the semidilute regime the stiff filaments can be considered as infinitely thin needles characterized solely by their length , leading to the dimensionless number density as only structural control parameter. We investigate the dynamics of a suspension of needles relying on stochastic simulations where the needles undergo diffusion with short-time rotational diffusion coefficient and short-time translational diffusion coefficient for parallel and perpendicular motion with respect to the needle axis. The hard-core interaction between the needles is handled following a pseudo-Brownian scheme Scala et al. (2007); Frenkel and Maguire (1983) and a geometry-adapted neighbor list to speed up the collision detection by up to two orders of magnitude (see Supplemental Material).
High entanglement occurs deep in the semidilute regime and the motion of a single needle becomes strongly confined to a narrow tube and both the rotational diffusion coefficient Doi (1975); Doi et al. (1984); Tao et al. (2006); Tse and Andersen (2013) as well as the translational diffusion coefficient perpendicular to the orientation of the rod are anticipated to scale as Teraoka and Hayakawa (1988); Szamel (1993). The rotational diffusion coefficient extracted from simulations is shown in Fig. 1(a) deep in the semidilute regime. The transport coefficient is suppressed by a factor of up to densities of , the highest density which has been considered before Doi et al. (1984); Tse and Andersen (2013). Our data extend to densities of such that the rotational diffusion is suppressed by another factor of . The Doi-Edwards prediction is nicely followed for the high densities . The data for the needle Lorentz system are qualitatively similar, in the highly entangled regime the scaling prediction is nicely corroborated, however, the static environment reduces the rotational diffusion by another factor of .
The time-dependent correlation function for the unit vector of orientation of the needle slows down drastically as the density is increased [Fig. 1(b)]. The shape of the longtime relaxation is well represented by an exponential, , characterized by the rotational diffusion coefficient . We have also measured higher-order orientational correlation functions , , where denotes the Legendre polynomials, and found that they also decay as with the same transport coefficient. Therefore, we conclude that at long times the pure orientational motion corresponds to simple diffusion on a sphere. At intermediate times when the needle is confined to its initial tube with diameter , the orientational correlation function displays a plateau close to unity, which we use to extract the tilt angle displayed in Fig. 1(d).
Furthermore, we have measured the mean-square displacements in the body-fixed frame. Since a collision does not affect the dynamics parallel to the needle, the corresponding mean-square displacement is trivial, in particular, the parallel diffusion coefficient is density independent, . In contrast, the mean-square displacement in the body-fixed frame perpendicular to the needle axis becomes diffusive only at long times [Fig. 1(c)], from which we obtain the perpendicular diffusion coefficient . The diffusion coefficient is suppressed by up to three orders of magnitude as can be inferred from Fig. 1(a), following the scaling prediction Teraoka and Hayakawa (1988); Szamel (1993). For high needle densities, the data for the mean-square displacement [Fig. 1(c)] display an extended plateau which we use to measure the tube diameter shown in Fig. 1(d), scaling as with needle density. For the needle Lorentz system, the data for the rotational diffusion and perpendicular translation are qualitatively identical to the needle liquid.
Translation-rotation coupling.— As a second step we investigate the relaxation of the translation-rotation coupling in terms of the tube model. The translational dynamics is anticipated to be transiently highly anisotropic, since diffusion occurs essentially along the needle axis. To quantify the translation-rotation coupling we measure the displacements of the needle center projected along or perpendicular to its initial axis and average its square. We denote the corresponding time derivatives
as projected time-dependent diffusion coefficients. For short times, the projected diffusion coefficients display pure parallel diffusion , since the needle axis remains aligned with its initial orientation [see Fig. 2(a)]. In contrast, for long times the initial orientation is forgotten and the motion becomes isotropic diffusion . For the projected perpendicular diffusion, in addition to the short-time plateau for times , a new plateau emerges for intermediate times , reflecting the strong suppression of the perpendicular motion by the confining tube. The tube constraint is relaxed only at much longer times and the projected perpendicular dynamics reaches isotropic diffusion . The needle Lorentz system shown in Fig. 2(b) displays the same phenomenology. We compare the translation-rotation coupling induced by the tube to a phantom needle diffusing freely in three-dimensional space, relying on the transport coefficients measured at long times. Thus, the phantom needle undergoes anisotropic diffusion in a rather peculiar effective medium such that the perpendicular diffusion and the rotational motion are suppressed by orders of magnitude with respect to the bare parallel diffusion . The corresponding projected time-dependent diffusion coefficients, displayed in Fig. 2(c), reproduce quantitatively the translation-rotation coupling both for a needle liquid as well as for the needle Lorentz system, provided the intermediate plateau becomes manifest. Hence, we conclude that also the coupling between the orientational and translational motion is well accounted for in terms of effective anisotropic diffusion. The coupling for the phantom needle can be worked out analytically, for example, by extending the methods developed by Han et al. Han et al. (2006) for the free two-dimensional anisotropic diffusion of ellipsoidal particles. For the three-dimensional case, we obtain
where characterizes the diffusional anisotropy. The theoretical curves are included in Fig. 2 and are in excellent agreement with the simulation data.
Intermediate scattering function.— The success of the tube model for the scaling behavior of the transport coefficients and the translation-rotation coupling encourages us to explore the full ramifications of the tube model for the experimentally accessible intermediate scattering function , providing temporal information on the needle dynamics on a length scale of . For simplicity we discuss only the motion of the needle center
for different wave numbers . The data (Fig. 3) display a statistical accuracy down to and extend over up to five significant decades in time. The data for different wave numbers collapse for short diffusion times , since the needle merely diffuses along its initial tube. For increasing wave number the data approach an emergent master curve with a power-law decay for long times. In the highly entangled regime the data for the needle Lorentz system are indistinguishable from the phantom needle at all wave numbers shown.
Since the intermediate scattering function for the needle liquid and the needle Lorentz system is in perfect agreement with the simulation for the phantom needle, we also elaborate explicit analytic expressions for within the tube model. The stochastic dynamics of the phantom needle is completely encoded in the conditional probability for the needle to displace by and change from an initial orientation to in lag time . The time evolution of the conditional probability is then governed by a Smoluchowski-Perrin equation Doi and Edwards (1978, 1999) and its spatial Fourier transform evolves according to
with and the rotational operator . The intermediate scattering function is then obtained by averaging over all initial orientations and integrating over all final orientations :
Here, we rely on a numerical evaluation of in terms of eigenfunctions of the operator on the right-hand side of Eq. (6) which turn out to be spheroidal wave functions Aragón S. and Pecora (1985) (see Supplemental Material). The integral in Eq. (7) can be performed for each term in the expansion of and reduces to an evaluation of the projection of the zeroth order spheroidal wave functions of degree onto the zeroth degree Legendre polynomial, :
with spheroidal eigenvalue . In essence, the intermediate scattering function is a weighted average of relaxing exponentials only, i.e. a completely monotone function Feller (1970). In particular, the longtime relaxation is purely exponential with a rate constant provided by the lowest eigenvalue of Eq. (6).
In the highly entangled regime becomes small and for wave numbers fulfilling the orientational motion in Eq. (6) is dominated by translational diffusion. Then, one can show that the spheroidal eigenvalues approach and the spheroidal wave functions are approximated by the eigenfunctions of the harmonic oscillator for fixed Meixner and Schäfke (1954). Then, for times only terms where the approximation is valid contribute to the sum of Eq. (8) such that the sum can be evaluated in closed form (see Supplemental Material) to
We have checked that Eq. (9) is an accurate representation of the intermediate scattering function for , yet the full solution faithfully represents the simulation data already at smaller values . For times , our analytic expression simplifies to
which rationalizes the power-law tail observed in the simulation data (Fig. 3) as well as the data collapse for different wave numbers and times .
Summary and conclusion.— The spatiotemporal dynamics of a needle in a highly entangled suspension has been characterized in terms of the transport coefficients, the time-dependent translation-rotation coupling, and the intermediate scattering function for a broad range of wave numbers and times. Relying on a geometry-adapted neighbor list, collisions can be handled efficiently for dynamical crowded systems. In this regime, our simulation data collapse to the dynamics of a single phantom needle with renormalized transport coefficients, thus confirming for the first time the relevance of the tube model for the spatiotemporal dynamics.
Our data demonstrate that the notion of the tube concept entails predictive power for the intermediate scattering function for both needle Lorentz systems and needle liquids. In the latter case, the mobile constraints cause a dilation of the effective tube width and the tube dynamically reorganizes by constraint release processes Doi and Edwards (1999), which occur on the same time scale as the polymer disengages from its initial tube. Comparing our simulation results for the two cases confirms that constraint release only shifts the prefactors of the scaling law , thus speeding up the dynamics. Hence, the dynamics in the needle liquid is identical to that of the needle Lorentz albeit at a lower density. Yet, the reduction of the many-body problem to the single freely diffusing phantom needle is robust, as suggested by Doi and Edwards Doi and Edwards (1999).
We have provided explicit analytic expressions in terms of spheroidal wave functions and their asymptotic harmonic oscillator behavior including the long algebraic tail, a fingerprint of the sliding motion within the tube. Interestingly, the tube model becomes quantitative already at densities and wave numbers where the tail is not yet manifest. Correspondingly, to account for the dynamics for all times and wave numbers, one has to resort to the full solution.
In principle, in scattering experiments not only the center of the needle is relevant but the entire rod contributes to the scattering signal. In essence, a form factor has to be included Doi and Edwards (1999) which considerably complicates the analytic solution Aragón S. and Pecora (1985). To obtain closed expressions for the modified intermediate scattering functions, Doi and Edwards Doi and Edwards (1978) also rely on the analogy to the harmonic oscillator but use an additional technical approximation which is circumvented in our analysis (see Supplemental Material). In simulations, the modified intermediate scattering function can be evaluated easily in terms of the phantom needle, yet, the overall picture remains unchanged since the dynamics is dominated by the back and forth motion in the narrow confining tube.
In the context of semiflexible polymers in the highly entangled regime, the tube concept has also proved itself to be the key insight Odijk (1983); Semenov (1986), although its ramifications have been tested mostly for static properties Romanowska et al. (2009); Morse (2001). It is anticipated to still capture the dynamics at coarse-grained time and length scales Ramanathan and Morse (2007); Nam et al. (2010); Keshavarz et al. (2016), yet the tube displays additional tube-width fluctuations Glaser et al. (2010); Wang et al. (2010) and a more complex tube renewal. Therefore, in the semiflexible case constraint release and tube renewal might admit a faster terminal relaxation such that the algebraic tail in the intermediate scattering function may be difficult to observe. One may speculate that the most important change is again a modified scaling law for the transport coefficients, while the semiflexible polymer still diffuses in and out of its tube along its own backbone Käs et al. (1994); Fakhri et al. (2010).
The one-dimensional sliding motion should, in fact, also be the most important ingredient for the dynamics of suspensions of nanofibers with finite aspect ratio and flexibility, and we anticipate that these systems also become accessible in simulations relying on our novel algorithm. With the rapid technological advancement to fabricate stiff nanorods and nanofibers Cassagnau et al. (2013); Kasimov et al. (2016), the semidilute regime of highly entangled suspensions has become finally into experimental reach.
Acknowledgements.This work was supported by the Austrian Ministry of Science BMWF as part of the UniInfrastrukturprogramm of the Focal Point Scientific Computing at the University of Innsbruck. We acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) Contract No. FR1418/5-1.
- Bausch and Kroy (2006) A. R. Bausch and K. Kroy, Nat. Phys. 2, 231 (2006).
- Liu et al. (2006) J. Liu, M. L. Gardel, K. Kroy, E. Frey, B. D. Hoffman, J. C. Crocker, A. R. Bausch, and D. A. Weitz, Phys. Rev. Lett. 96, 118104 (2006).
- Koenderink et al. (2006) G. H. Koenderink, M. Atakhorrami, F. C. MacKintosh, and C. F. Schmidt, Phys. Rev. Lett. 96, 138307 (2006).
- Wong et al. (2004) I. Y. Wong, M. L. Gardel, D. R. Reichman, E. R. Weeks, M. T. Valentine, A. R. Bausch, and D. A. Weitz, Phys. Rev. Lett. 92, 178101 (2004).
- Lin et al. (2007) Y.-C. Lin, G. H. Koenderink, F. C. MacKintosh, and D. A. Weitz, Macromolecules 40, 7714 (2007).
- Lettinga and Grelet (2007) M. P. Lettinga and E. Grelet, Phys. Rev. Lett. 99, 197802 (2007).
- Addas et al. (2004) K. M. Addas, C. F. Schmidt, and J. X. Tang, Phys. Rev. E 70, 021503 (2004).
- Koenderink et al. (2004) G. H. Koenderink, S. Sacanna, D. G. A. L. Aarts, and A. P. Philipse, Phys. Rev. E 69, 021804 (2004).
- Cassagnau et al. (2013) P. Cassagnau, W. Zhang, and B. Charleux, Rheol. Acta 52, 815 (2013).
- Doi (1975) M. Doi, J. Phys. (France) 36, 607 (1975).
- Doi and Edwards (1978) M. Doi and S. F. Edwards, J. Chem. Soc., Faraday Trans. 2 74, 560 (1978).
- Doi et al. (1984) M. Doi, I. Yamamoto, and F. Kano, J. Phys. Soc. Jpn. 53, 3000 (1984).
- Tao et al. (2006) Y.-G. Tao, W. K. den Otter, J. K. G. Dhont, and W. J. Briels, J. Chem. Phys. 124, 134906 (2006).
- Tse and Andersen (2013) Y.-L. S. Tse and H. C. Andersen, J. Chem. Phys. 139, 044905 (2013).
- Teraoka and Hayakawa (1988) I. Teraoka and R. Hayakawa, J. Chem. Phys. 89, 6989 (1988).
- Szamel (1993) G. Szamel, Phys. Rev. Lett. 70, 3744 (1993).
- Löwen (1994) H. Löwen, Phys. Rev. E 50, 1232 (1994).
- Tucker and Hernandez (2010) A. K. Tucker and R. Hernandez, J. Phys. Chem. A 114, 9628 (2010).
- Sussman and Schweizer (2011) D. M. Sussman and K. S. Schweizer, Phys. Rev. Lett. 107, 078102 (2011).
- Yamamoto and Schweizer (2015) U. Yamamoto and K. S. Schweizer, ACS Macro Lett. 4, 53 (2015).
- Frenkel and Maguire (1981) D. Frenkel and J. F. Maguire, Phys. Rev. Lett. 47, 1025 (1981).
- Frenkel and Maguire (1983) D. Frenkel and J. F. Maguire, Mol. Phys. 49, 503 (1983).
- Otto et al. (2006) M. Otto, T. Aspelmeier, and A. Zippelius, J. Chem. Phys. 124, 154907 (2006).
- Höfling et al. (2008a) F. Höfling, E. Frey, and T. Franosch, Phys. Rev. Lett. 101, 120605 (2008a).
- Tucker and Hernandez (2011) A. K. Tucker and R. Hernandez, J. Phys. Chem. B 115, 4412 (2011).
- Fixman (1985) M. Fixman, Phys. Rev. Lett. 54, 337 (1985).
- Cobb and Butler (2005) P. D. Cobb and J. E. Butler, J. Chem. Phys. 123, 054908 (2005).
- Zhao and Wang (2013) T. Zhao and X. Wang, Polymer 54, 5241 (2013).
- Moreno and Kob (2004) A. J. Moreno and W. Kob, Europhys. Lett. 67, 820 (2004).
- Höfling et al. (2008b) F. Höfling, T. Munk, E. Frey, and T. Franosch, Phys. Rev. E 77, 060904 (2008b).
- Munk et al. (2009) T. Munk, F. Höfling, E. Frey, and T. Franosch, Europhys. Lett. 85, 30003 (2009).
- Aragón S. and Pecora (1985) S. R. Aragón S. and R. Pecora, J. Chem. Phys. 82, 5346 (1985).
- Scala et al. (2007) A. Scala, Th. Voigtmann, and C. De Michele, J. Chem. Phys. 126, 134109 (2007).
- Han et al. (2006) Y. Han, A. M. Alsayed, M. Nobili, J. Zhang, T. C. Lubensky, and A. G. Yodh, Science 314, 626 (2006).
- Doi and Edwards (1999) M. Doi and S. F. Edwards, The Theory of Polymer Dynamics (Oxford University Press, Oxford, 1999).
- Feller (1970) W. Feller, An Introduction to Probability Theory and Its Applications, 1st ed., Vol. 2 (John Wiley & Sons, Inc., New York, 1970).
- Meixner and Schäfke (1954) J. Meixner and F. W. Schäfke, Mathieusche Funktionen und Sphäroidfunktionen (Springer, Berlin, 1954).
- Odijk (1983) T. Odijk, Macromolecules 16, 1340 (1983).
- Semenov (1986) A. N. Semenov, J. Chem. Soc., Faraday Trans. 2 82, 317 (1986).
- Romanowska et al. (2009) M. Romanowska, H. Hinsch, N. Kirchgeßner, M. Giesen, M. Degawa, B. Hoffmann, E. Frey, and R. Merkel, Europhys. Lett. 86, 26003 (2009).
- Morse (2001) D. C. Morse, Phys. Rev. E 63, 031502 (2001).
- Ramanathan and Morse (2007) S. Ramanathan and D. C. Morse, Phys. Rev. E 76, 010501 (2007).
- Nam et al. (2010) G. Nam, A. Johner, and N.-K. Lee, J. Chem. Phys. 133, 044908 (2010).
- Keshavarz et al. (2016) M. Keshavarz, H. Engelkamp, J. Xu, E. Braeken, M. B. J. Otten, H. Uji-i, E. Schwartz, M. Koepf, A. Vananroye, J. Vermant, R. J. M. Nolte, F. D. Schryver, J. C. Maan, J. Hofkens, P. C. M. Christianen, and A. E. Rowan, ACS Nano 10, 1434 (2016).
- Glaser et al. (2010) J. Glaser, D. Chakraborty, K. Kroy, I. Lauter, M. Degawa, N. Kirchgeßner, B. Hoffmann, R. Merkel, and M. Giesen, Phys. Rev. Lett. 105, 037801 (2010).
- Wang et al. (2010) B. Wang, J. Guan, S. M. Anthony, S. C. Bae, K. S. Schweizer, and S. Granick, Phys. Rev. Lett. 104, 118301 (2010).
- Käs et al. (1994) J. Käs, H. Strey, and E. Sackmann, Nature (London) 368, 226 (1994).
- Fakhri et al. (2010) N. Fakhri, F. C. MacKintosh, B. Lounis, L. Cognet, and M. Pasquali, Science 330, 1804 (2010).
- Kasimov et al. (2016) D. Kasimov, T. Admon, and Y. Roichman, Phys. Rev. E 93, 050602 (2016).
- Chirikjian (2009) G. S. Chirikjian, Stochastic Models, Information Theory, and Lie Groups, Vol. 1 (Birkhäuser, Boston, 2009).
- Øksendal (2010) B. Øksendal, Stochastic Differential Equations: An Introduction with Applications, 6th ed. (Springer, Berlin, 2010).
- Donev et al. (2005) A. Donev, S. Torquato, and F. H. Stillinger, J. Comput. Phys. 202, 737 (2005).
- Hansen and Walster (2004) E. Hansen and G. W. Walster, Global Optimization Using Interval Analysis, 2nd ed. (Marcel Dekker, New York, 2004).
- Berne and Pecora (2000) B. J. Berne and R. Pecora, Dynamic Light Scattering: With Applications to Chemistry, Biology, and Physics (Dover Publications, Meniola, 2000).
- (55) DLMF, “NIST digital library of mathematical functions,” http://dlmf.nist.gov/, Release 1.0.10 of 2015-08-07, online companion to Olver et al. (2010).
- Olver et al. (2010) F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, eds., NIST Handbook of Mathematical Functions (Cambridge University Press, New York, NY, 2010) print companion to DLMF ().
- Hodge (1970) D. B. Hodge, J. Math. Phys. 11, 2308 (1970).
- Falloon et al. (2003) P. E. Falloon, P. C. Abbott, and J. B. Wang, J. Phys. A 36, 5477 (2003).
Appendix A Supplemental Material
a.1 Stochastic simulation
We consider the diffusion of an infinitely thin needle of length with the short-time rotational diffusion coefficient and the short-time translational diffusion coefficient for parallel and perpendicular motion. At every time , the state of the needle is completely described by the position of the center of the needle and the unit vector of the orientation . The starting point for the stochastic simulation are the following Langevin equations in the Itō interpretation for the change in orientation and position Chirikjian (2009); Øksendal (2010):
The independent Gaussian white noise processes and with zero mean and covariance determine the stochastic dynamics. The dyadic product is denoted by .
We implement the Langevin equations using a discrete fixed Brownian time step and assign to the needle the random pseudo-angular velocity and pseudo-translational velocity according to
The Gaussian white noise process are expressed via the independent and normal distributed random variables and with zero mean and unit variance.
During a Brownian time step , the needle evolves ballistically according to the propagation rules
We simulate both needle liquids as well as needle Lorentz systems. In the former case, we move every needle subsequently and fix all other needles, whereas in Lorentz systems only a single needle diffuses in a frozen array of other needles. For needle liquids we use a Brownian time-scale of and observe the dynamics over six decades in time, whereas for needle Lorentz systems we use and generate trajectories over ten decades in time.
The hard-core interaction of the moving needle with the other ones is handled following a pseudo-Brownian scheme Scala et al. (2007); Frenkel and Maguire (1983). Here we assume, that during a Brownian time step the needle propagates ballistically according to the propagation rules (13) in continuous time . Collisions between two needles are only possible if the distance of the intersection point of both needle centers is smaller than half the needle length. The number of possible candidates can be greatly reduced by keeping an efficient neighbor list adapted to the needle problem. It turns out that considering a surrounding cylinder is way more advantageous than a sphere conventionally used for spherical particles. Geometry-adapted neighbor lists have been considered before for ellipsoidal particles Donev et al. (2005), however its full advantage becomes manifest for very high aspect ratios, in particular needles leading to a speedup of up to two orders of magnitude for the highest densities considered. The efficiency is determined by the diameter of the cylinder which can be decreased upon considering smaller Brownian time-scales .
The fixed collision candidate needle with orientation intersects the rotational plane of the moving needle perpendicular to for times which fulfill
where the distance between the needle centers at time is denoted by . In addition to that, the intersection point of the obstacle needle with the rotational plane has to traverse the rotational disk of radius of the moving needle. This is described as solution of the quadratic equation
where we denote the distance of the center of the moving needle with the intersection point in the rotational plane at time as . As a result the time interval for a possible collision reduces to the interval . If the new interval is nonempty, collision times can be determined explicitly as solutions of the equation
via an interval Newton method Hansen and Walster (2004). Moreover, a correctly rounded math library has to be used for the transcendental functions in order not to miss collisions.
Upon colliding, the new velocities for translation and rotation are determined by energy, momentum and angular momentum conservation. We assume that the center of mass coincides with the geometric center of the needle. For smooth needles Frenkel and Maguire (1983); Otto et al. (2006), the momentum transfer onto the moving needle is directed perpendicular to both needle orientation and has a magnitude of
where the distance between the point of contact and the center of the moving needle has been introduced. In particular, the pseudo-velocity parallel to the needle does not change. The ratio of the mass to the moment of inertia remains a free parameter so far. We define a pseudo-temperature for the rotational degrees of freedom by , and similarly for the perpendicular translational degrees of freedom . From the equation for the pseudo-velocities (Eq. 12) we determine the averages and . In order to avoid an average flow of energy between the rotational and the perpendicular translational degrees for freedom we impose which yields
In general, the transport coefficients of the needle are not independent and their relation is established by a hydrodynamic calculation. For the simulation, we use the transport coefficients of a slender rod characterized by as well as Doi and Edwards (1999).
a.2 Intermediate scattering function of the phantom needle
The stochastic dynamics of the phantom needle is translationally invariant in time and space and described by the the propagator which represents the conditional probability for a displacement and change of orientation from to in lag time . The time evolution of the propagator with initial condition is determined by the Smoluchowski equation
The main quantity of interest for the dynamics is given by the self-intermediate scattering function
as an average over all initial orientations and integration over all final orientations of the spatial Fourier transform of the propagator :
The Smoluchowski equation for the spatial Fourier transform assumes the form
We choose a coordinate system such that the wave vector is along the -direction, . Then, the Smoluchowski-Perrin equation in spherical coordinates reads
where the relation establishes the connection to the polar angle .
This partial differential equation is solved by a separation of variables and the full solution (Ref. Aragón S. and Pecora (1985)) is obtained as an expansion in terms of spheroidal wave functions of degree and order :
where the real parameter and the characteristic decay constants depend only on the diffusion coefficients , , and and magnitude of the wave vector . The exact relations for both parameters are determined as solutions of the spheroidal wave equation
where in our case and the spheroidal eigenvalue determines the decay rate DLMF (); Olver et al. (2010). Only contributions of order remain after averaging over all initial orientations and integration over all final orientations:
Since the spheroidal wave functions display the symmetry property, , only functions with even degree contribute to the intermediate scattering function after integration over , we can restrict the sum to even . Expanding the spheroidal wave functions in terms of Legendre polynomials ,
and perform the remaining integration via the orthonormality condition of the Legendre polynomials:
Hence, the intermediate scattering function reduces to a sum over exponentially decaying functions with decay constants and amplitudes given by the first coefficient of the spheroidal wave function , the overlap between the zeroth order spheroidal wave function of degree and the zeroth degree Legendre polynomial:
Methods for the computation of the eigenvalues and the coefficients are readily available in the literature Hodge (1970); Falloon et al. (2003); DLMF (). The evaluation of the sum is truncated when the spheroidal eigenvalue becomes positive and the sum is close to unity for time zero. For the highest densities and wave numbers, , more than terms contribute and multi-precision floating point arithmetic is mandatory Falloon et al. (2003).
In the highly entangled regime and thus , the spheroidal eigenvalue approaches and the spheroidal wave functions for fixed converge in mean square to the eigenfunctions of the harmonic oscillator Meixner and Schäfke (1954):
where denote the Hermite polynomials. By analogy to the harmonic oscillator, the width of the eigenfunctions approaches , and the integral in Eq.(28) can be safely extended to the real line for . The remaining integral for even can be evaluated exactly by employing the generating function of the Hermite polynomials in terms of a generalized binomial coefficient
Thus, by Newton’s binomial series, the sum [Eq. (29)] can be evaluated in closed form and the intermediate scattering functions assumes the approximate form
for times , such that only terms contribute to the sum. In particular, one infers the terminal relaxation rate of the intermediate scattering function. Furthermore, for times the result can be simplified further
Thus in the regime the intermediate scattering function displays a characteristic power-law decay with exponent in a broad time window.
The origin of the tail can be easily understood as a consequence of a pure sliding motion. Since the tail becomes manifest for times where the needle barely rotates, , we approximate in this regime, such that the Smoluchowski equation simplifies to a first order differential equation:
The solution is directly obtained as , describing the motion of a needle where the initial orientation is preserved for all times. Then the intermediate scattering function can be evaluated in closed form
where denotes the complementary error function. Hence, for large , we recover Eq. (34) .
Let us compare to the solution strategy of Doi and Edwards (ignoring the form factor for the infinitely thin needle) Doi and Edwards (1978). There, the Smoluchowski-Perrin equation [Eq. (23)] is approximated by a time-dependent harmonic oscillator equation after integration over the azimuthal angle which is correct for polar angles close to the equator. Then, they obtain the intermediate scattering function
where the integral extends over all polar angles of the Green function of a harmonic oscillator
Using the generating function of the Hermite polynomials, the Green function can be expressed also as
which is obtained from the exact solution by replacing the eigenvalue and eigenfunctions by their harmonic oscillator analogue irrespective of the degree . However, as discussed above, this replacement is only valid for and correspondingly, the Green function is accurate only for times and polar angles close to the equator.
In the original work Doi and Edwards (1978), the integral in Eq. (37) can not be performed in closed form rather an additional approximation is introduced such that the needle keeps its initial polar angle. We note that one can formally obtain our result Eq. (33) by extending the integrals in Eq. (37) to the real line, which amounts to interchanging limits of the infinite sum with the integrals over the polar angles, irrespective of convergence.