Structure formation of surfactant membranes under shear flow
SStructure formation of surfactant membranes under shear flow
Hayato Shiba ∗ and Hiroshi Noguchi † Institute for Solid State Physics, University of Tokyo, Chiba 277-8581, Japan
Gerhard Gompper
Theoretical Soft Matter and Biophysics, Institute of Complex Systems and Institute for Advanced Simulation,Forschungszentrum J¨ulich, D-52425 J¨ulich, Germany (Dated: October 29, 2018)Shear-flow-induced structure formation in surfactant-water mixtures is investigated numericallyusing a meshless-membrane model in combination with a particle-based hydrodynamics simulationapproach for the solvent. At low shear rates, uni-lamellar vesicles and planar lamellae structures areformed at small and large membrane volume fractions, respectively. At high shear rates, lamellarstates exhibit an undulation instability, leading to rolled or cylindrical membrane shapes orientedin the flow direction. The spatial symmetry and structure factor of this rolled state agree withthose of intermediate states during lamellar-to-onion transition measured by time-resolved scattingexperiments. Structural evolution in time exhibits a moderate dependence on the initial condition.
PACS numbers: 83.50.Ax, 83.80.Qr, 87.16.dt, 87.15.A-
I. INTRODUCTION
Surfactant molecules in water self-assemble into vari-ous structures such as micelles and bilayer membranes,which display a rich variety of rheological properties un-der flow. Even if a basic structure remains to be a bi-layer membrane, its mesoscale structure can assume sev-eral different states, such as fluid L α or ripple P β phase.Under shear flow, lamellae can be oriented parallel orperpendicular to the shear-gradient direction. Diat andRoux first discovered 20 years ago closely-packed multi-lamellar vesicle (MLV) structures, so-called the onionphase, in nonionic surfactant-water mixtures under shearflow. In the last two decades, this onion structure hasbeen studied experimentally using light, neutron, and X-ray scattering, and also by freeze-fractureelectron microscopy, and the rheo-NMR method. Its rheology has been of large interest.
Typi-cally, a critical shear rate ˙ γ c separates the lamellae andonion phases, where the latter phase consists of mono-disperse onions containing hundreds of lamellar layers.The onion radius R ( ˙ γ ) is reversible and can be describedby a unique decreasing function of the shear rate ˙ γ . Time-resolving small-angle neutron scattering experi-ments have revealed that a two-dimensional intermediatestructure is formed during the lamellar-to-onion transi-tion with increasing shear rate. A cylindrical or wavylamellar structure was speculated to be the transient in-termediate structure, but could not be distinguished fromthe scattering pattern alone. Recent small-angle X-rayscattering experiments with increasing temperature atconstant shear rate also indicate a similar pattern aroundthe lamellar-to-onion transition. Thus, there are someexperimental evidences of a transient state, but its struc-ture is still under debate. An alternative experimentalapproach to gain insight into the structural changes isto characterize defects observed in the lamellar state formoderate shear rates, both in surfactant membranes and in thermotropic liquid crystals.
It is also worthmentioning that stable cylindrical structures on a ten-micrometer length scale are observed when strong shearflow is applied in the lamellar-sponge coexistence state. Several theoretical attempts have been made to tacklethis complex problem of structural evolution under shearflow, which consider either instability of the lamel-lar phase due to undulations or the break-up ofdroplets.
Recently, a “dynamical” free energy ofMLV under shear flow has been proposed, which takesinto account the slow modes induced by the solvent be-tween the membranes together with their bending andstretching forces. The scaling relations for the MLV sizeand the terminal shear rate are predicted in agreementwith the experiments of Refs. 2,3. In these theories, whilethe hydrodynamic effects of the solvent are taken intoaccount, the analysis is performed for geometrically sim-ple structures, asp spherical onions or planar lamellae.Thus, the kinetic process of the transformation from thelamellar to the onion phase could not be investigatedtheoretically so far. In this paper, we study the detailedstructural evolution of surfactant membranes under sim-ple shear flow using large-scale particle simulations.A few simulations have been performed for the for-mation of lamellar phases in shear flow, while onionformation has never been addressed so far. Orientedlamellae have been obtained in simulations of a coarse-grained molecular model for lipids, while defect dy-namics has been investigated in simulations of a phase-field model of a smectic-A system. Onion and in-termediate states have large-scale structures of the or-der of micrometers, which are beyond typical lengthscales accessible to molecular dynamics (MD) simula-tions of coarse-grained surfactant molecules. In ourstudy, we employ a meshless-membrane model, where a membrane particle represents not a surfactantmolecule but rather a patch of bilayer membrane –in order to capture the membrane dynamics on a mi- a r X i v : . [ c ond - m a t . s o f t ] A p r crometer length scale. This model is well suitable tostudy membrane dynamics accompanied by topologi-cal changes. Alternatively, membranes can be mod-eled as triangulated surfaces, which require how-ever discrete bond reconnections to describe topologicalchanges. After the first meshless-membrane modelwas proposed in 1991, several meshless-membranemodels have been developed. In contrast toother meshless-membrane approaches, our models are capable of separately controlling the membrane bend-ing rigidity κ and the line tension Γ of membraneedges. Previously, we have combined our meshless-membrane model with multi-particle collision (MPC)dynamics, a particle-based hydrodynamic simulationtechnique. With MPC, the hydrodynamic interactionsare properly taken into account, but due to the frictionalcoupling of membrane and solvent, solvent particles canpenetrate through the membrane. We here extend themeshless-membrane model into explicit solvent simula-tion model, in which the fluid particles interact with eachother and with the membrane via short-ranged repulsivepotentials, so that the solvent can hardly penetrate themembrane, and simulate it with dissipative particle dy-namics (DPD), another hydrodynamics simulation tech-nique.The simulation model and methods are introduced inSec. II. Then basic membrane properties, including bend-ing rigidity and line tension, are described in Sec. III. InSec. IV, structure formation in surfactant-water mixturesunder shear flow is investigated for variety of shear rates˙ γ and membrane volume fractions ϕ . At high ˙ γ and high ϕ , a novel structure of rolled-up membranes is found,which are oriented in the flow direction. A summary andsome perspectives are given in Sec. V. II. MODEL AND METHODSA. Coarse-Grained Model and InteractionPotentials
To simulate the structure formation in surfactant-membranes systems, we employ a meshless-membranemodel with explicit solvent. In this model, two typesof particles A and B are employed, which denote mem-brane and solvent particles, respectively. The number ofthese particles is N A and N B , which defines the parti-cle density φ = ( N A + N B ) /V – where V is the volumeof the simulation box – and the membrane volume frac-tion ϕ = N A / ( N A + N B ). The particles interact via apotential U , which consists of repulsive, attractive, andcurvature interactions, U rep , U att , and U α , respectively, Uk B T = (cid:88) i 25 ln { − ρ i − ρ ∗ )] } − C, (3)which is a function of the local density ρ i of the mem-brane particles defined by ρ i = (cid:88) i ∈A f cut ( r ij /σ ) . (4) C = 0 . 25 ln[1 + exp(4 ρ ∗ )] is chosen such that U att (0) =0. The cutoff function f cut in Eq. (4) is a C ∞ functionrepresented as f cut ( s ) = (cid:40) exp (cid:104) a (1 + | s | /s cut ) n − ) (cid:105) ( s < s cut )0 ( s ≥ s cut ) (5)where n = 12, a = ln(2) { ( s cut /s half ) n − } with s cut = 2 . s half = 1 . ρ ∗ = 6 to study a fluidmembrane in an explicit solvent.The curvature potential U α = (cid:88) i ∈A α pl ( r i ) (6)is introduced to incorporate the membrane bending rigid-ity. Here, the aplanarity α pl provides a measure for thedegree of deviation of membrane particle alignment froma planar reference state. It is defined by α pl = 9 D w T w M w = 9 λ λ λ ( λ + λ + λ )( λ λ + λ λ + λ λ ) , (7)where λ ≤ λ ≤ λ are the three eigenvalues ofthe local gyration tensor a αβ of the membrane nearparticle i , which is defined as a αβ = (cid:80) j ∈A ( α j − α G )( β j − β G ) w cv ( r ij ), with α, β ∈ { x, y, z } . Here, r G = (cid:80) j ∈A r j w cv ( r ij ) / (cid:80) j ∈A w cv ( r ij ) is the locally weightedcenter of mass, and w cv ( r ij ) is a truncated Gaussian func-tion w cv ( r ij ) = (cid:40) exp (cid:16) ( r ij /r ga ) ( r ij /r cc ) n − (cid:17) ( r ij < r cc )0 ( r ij ≥ r cc ) (8)which is smoothly cut off at r ij = r cc . The constants areset as n = 12 , r ga = 1 . σ , and r cc = 3 . σ .In all our simulations, the volume is chosen such thatthe number density φ of the particles is constant, φ = N/V = 0 . σ − , N = N A + N B . (9)For higher solvent densities, the system is closer to themelting point, and strong attractive interaction betweenmembranes sheets have been found (see Appendix A fordetails). The density φ = 0 . σ − is chosen in order toavoid such an attraction. B. Thermostats We simulated the membranes in the NVT ensembleunder shear flow. To keep the temperature constant,we employ the dissipative particle dynamics (DPD)thermostat, in which friction and noise forces are ap-plied to the relative velocities of pairs of neighboring par-ticles. Thus, linear and angular momentum is conserved,which implies that the system shows hydrodynamic be-havior on sufficiently large length and time scales. Theequation of motion for the i th particle is given by m d v i dt = − ∂U∂ r i + (cid:88) j (cid:54) = i {− w ij ( v i − v j ) · ˆ r ij + √ w ij ξ ij ( t ) } ˆ r ij , (10)where ˆ r ij = r ij /r ij . Here, the weight function w ij is w ij ( r ij ) = γθ ( Aσ − r ij ) with A = 2 . 7, where θ ( r )is the Heaviside step function. This type of weighthas also been used in the Lowe-Andersen thermostat. The scheme is discretized with Shardlow-S1 splittingmethod, whose time step ∆ t b = 0 . t is different fromthat for contributions from the molecular interactions,where ∆ t = 0 . t ( t = m/γ is simulation time unit).In the following, we use γ = √ mk B T /σ , where k B T isthe thermal energy unit. Although the DPD thermostatis usually employed for soft interaction potentials, itcan also be employed for systems with steeper potentials,such as the Weeks-Chandler-Andersen potential. Shear flow with velocity v x = ˙ γz in the x -direction andgradient in the z -direction is imposed by Lees-Edwardsboundary condition. The code is optimized for use on aparallel computer architecture by domain decomposition(see Appendix B). III. MODEL PROPERTIES The meshless-membrane model with explicit solventis expected to exhibit very similar equilibrium proper-ties as the original implicit-solvent version. We now -1 0 1 2 3 4 5 1.1 1.2 1.3 1.4 1.5 (cid:97) s A xy /N A k (cid:95) = 2.5 5.010.020.0 FIG. 1: Area dependence of surface tension γ s for the ex-plicit solvent meshless membrane model is plotted for k α =2 . , , , and 20 with (cid:15) = 4. confirm this dependence of surface tension, line tension,and bending rigidity on the control parameters for theexplicit-solvent model.First, a planar fluctuating membrane is simulated withthe particle numbers N A = 1600 (and N = 48 000 or64 000). For various membrane projected areas A xy = L x L y , the surface tension γ s = (cid:104) P zz − ( P xx + P yy ) / (cid:105) L z (11)is investigated. Here, P µν is the pressure tensor given by P µν = N (cid:88) i =1 (cid:18) mv µi v νi − µ i ∂U∂ν i (cid:19) (cid:46) V, (12)where { µ, ν } ∈ { x, y, z } , v i = ( v xi , v yi , v zi ), and the sum istaken over all the particles including the solvent compo-nent. In calculating P µν , the periodic image µ i + nL µ nearest to the other interacting particles is employed,when the potential interaction crosses one of the peri-odic boundaries.Figure 1 shows the dependence of γ s on the projectedmembrane area A xy for k α = 2 . , , 10, and 20, with (cid:15) = 4. For γ s (cid:39) 0, the intrinsic area A is larger than A xy due to the membrane undulations; buckling of themembrane occurs for A < A xy (the flat region at γ s < For tension-less membranes with N A = 1600,the projected area A xy is given by A xy = a xy N A with a xy = 1 . , . , . , . 35 for k α = 2 . , , , 20, re-spectively. The projected membrane area increases withincreasing k α , because both membrane bending fluctua-tions and protrusions are suppressed at larger k α .Figure 2 shows the bending rigidity κ as a functionof k α . Here, the bending rigidity is estimated from theheight spectrum (cid:104)| h ( q ) | (cid:105) = k B Tγ s q + κq . (13)of the tension-less membrane ( γ s = 0). In calculating (cid:104)| h ( q ) | (cid:105) , the raw positional data of the membrane par-ticles ( r i , i ∈ A ) is employed. Because of the slow (cid:103) / k B T k (cid:95) -1 -1 < | h ( q ) | > q/ (cid:47) q -4 FIG. 2: Dependence of the bending rigidity κ on k α is plottedfor (cid:15) = 4, estimated with the use of a tension-less planarmembrane at N A = 1600 through Eq. (13). As shown in theinset, the height spectrum of the membrane exhibits a q − spectrum. (cid:75)(cid:109) / k B T (cid:161) k (cid:95) =510 FIG. 3: (cid:15) dependence of the line tension Γ of explicit-solventmeshless-membrane model, calculated for k α = 5 and 10 withthe use of Eq. (14). dynamics of long-wavelength height fluctuations, it be-comes more time-consuming to obtain precise data thanfor the implicit-solvent model. Therefore, (cid:104)| h ( q ) | (cid:105) ismeasured using 16 independent runs for each k α , and theaveraged spectrum is then fitted to Eq. (13) in the range q < . σ − . As in the implicit model, the bending rigid-ity is found to be proportional to k α , which demonstratesthe controllability of κ in our model.In Fig. 3, the (cid:15) dependence of the line tension Γ isdisplayed for k α = 5 and 10 for a membrane strip withtwo edges with lengths equal to L x with N A = 1600.Here, Γ is determined via the relation Γ = (cid:104) ( P yy + P zz ) / − P xx (cid:105) L y L z / . (14)At a small (cid:15) , Γ is proportional to (cid:15) and almost inde-pendent of k α , similarly to the implicit-solvent meshlessmembrane model . While Γ increases linearly up to (cid:15) = 8 for k α = 10, it levels off and saturates at (cid:15) = 4 for k α = 5. When (cid:15) exceeds the value where Γ satu-rates — a value which becomes larger for larger k α —the membrane particles prefer to reside at the edges, be-cause the curvature force is not strong enough to avoidaggregation due to the stronger attractive forces.For our simulations of membrane ensembles in shearflow, we choose the model parameters k α = 5 and (cid:15) = 4,where the membrane has bending rigidity κ/k B T =11 ± 1. With the estimated values for the line tensionΓ σ/k B T = 4 . ± . 2, the relaxation time scale of struc-tural transitions can be characterized by τ = ηR c κ , (15)where η is the solvent viscosity, and R c is the criticalradius of a flat disk. Assuming a flat disk with radius R and a vesicle with the same membrane area, we obtainthe corresponding elastic free energies F d = 2 πR Γ and F s = 8 π ( κ + ¯ κ/ (cid:39) πκ ; thus, a disk exhibits transitionto closed vesicle shape if the radius is around R c = 2 κ/ Γ ∼ (5 . ± . σ. (16)With the solvent viscosity η = (2 . ± . × m ( σt ) − obtained from a solvent-only simulation, an estimationof the relaxation time yields τ = 28 . t . In the follow-ing, the time will be measured in unit of τ . Since themembrane thickness, typically around 5nm for non-ionicsurfactants like polyethylene-glycol-ethers C n E m , corre-sponds to the size σ of the membrane particles in oursimulation, τ is equivalent to about 0 . µ s (with the vis-cosity of water at 300K, η w (cid:39) . · s). IV. STRUCTURE FORMATION WITH ANDWITHOUT SHEAR FLOW We now employ the meshless-membrane modelwith explicit solvent to study structure formation insurfactant-water mixtures, both with and without shearflow. In all simulations, the total particle number is fixedat N = N A + N B = 960 000, and thus, the system sizeis a cubic box with side length L = 114 . σ . Simula-tions are performed for various membrane volume frac-tion ϕ = N A /N , with ϕ = 0 . . . . . . × MD steps, corresponding to2 . × τ . All the particles of both species are initiallydistributed randomly in the simulation box. Averagesare calculated over the last 2 × steps, where the sys-tem is assumed to have reached a stationary state. Afterbriefly explaining the structures obtained by equilibriumsimulations without shear in Sec. IV A, we present resultsfor the structure formation in a system under linear shearflow in Secs. IV B and IV C. FIG. 4: Snapshots of the configuration of membrane particles for ϕ = 0 . γ = 0) and (b) under shearflow with shear rate ˙ γτ = 0 . ϕ = 0 . 125 with ˙ γτ = 0 . ϕ = 0 . γτ = 0 . γτ = 0 . γτ = 0 . v = ˙ γz e x . Solvent particles are not displayed. A. Mesophases in Thermal Equilibrium At a low membrane volume fraction ϕ = 0 . N c = πR c /a (cid:39) 75. At a higher membrane volume frac-tion ϕ = 0 . 125 ( N A = 120 000), the membrane surfacesis found to percolate through the whole system via theperiodic boundaries. This behavior is not unexpected,because each disk with radius R c covers a region of vol-ume v c = (2 R c ) by rotational diffusion, and thereforethese disks should overlap and merge when more than n c ∼ V /v c = 1200 vesicles are present, which exceeds thenumber of vesicles with critical size ( N A /N c ).For ϕ = 0 . 25 and 0.3125, periodic lamellar states areformed owing to the repulsive interactions between themembranes, as shown in Fig. 5(a). The membranes arecurved to fill up space with random orientations. Inthermal equilibrium, these lamellar layers would prob-ably have a unique orientation throughout the system;this well-ordered state is difficult to reach in simulationsbecause the structural relaxation time well exceeds theaccessible simulation time scale. The 3D structure factor of the membrane density is calculated as S A ( q ) = (cid:90) d r e i q · r (cid:104) δ ˆ n A ( r ) δ ˆ n A ( ) (cid:105) , (17)where δ ˆ n A ( r ) = (cid:80) i ∈A σ δ ( r − r i ) − φ is the local devi-ation of the membrane particle density from its average.Figure 5(b) shows S A ( q ) for ϕ = 0 . q z = 0. Peaksarising from the inter-lamellar distance are observed at | q | = q = 1 . σ − and q = 3 . σ − = 2 q with heights9.8 and 0.78, respectively. The former corresponds tolength of L = 2 π/q = 3 . σ , which provides a preciseestimate of the interlayer distance. B. Dynamic Phase Diagram under Shear Flow In shear flow, the vesicle and lamellar states depend onthe concentration ϕ and the shear rate ˙ γ , as displayed inFig. 6. At ϕ = 0 . γτ < . 02. At a higher shear rate, themembrane particles tend to assemble into plate-like mem-brane disks, which then align parallel to the shear flowdirection, as demonstrated by the comparison of snap-shots in Figs. 4 (a) and (b). At ϕ = 0 . FIG. 5: (a) Snapshot of the configuration of membrane par-ticles for ϕ = 0 . γ = 0). To facili-tate visualization, only a thin planar slice is displayed. Sol-vent particles are not shown. (b) Structure factor S A [ q =(0 , q y , q z )] for ϕ = 0 . γ = 0. -2 -1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 (vesicle + plates) ˙ γ τ ϕ lam.+cyl und. lam. roll vesicle lamellae FIG. 6: Dynamic phase diagram of the explicit-solventmeshless-membrane model as a function of the volume frac-tion ϕ of the membrane component and the shear rate ˙ γ . Theparameters are k α = 5, (cid:15) = 4, φ = 0 . σ − , and N = 960 , ˙ γτ < . 06, as shown in Fig. 4(c). There is a larger varietyof phases at ϕ = 0 . γτ < . 1, thesystem is a mixture of lamellae and cylinders as shown inFig. 4(d); at ˙ γτ = 0 . γτ (cid:38) . ϕ , the lamellar states align perpendicularlyto the shear-gradient ( z ) direction in the regime of lowshear rates ( ˙ γτ < . γ , they exhibit an in-stability to a rolled-up shape whose axis is parallel to theflow direction, which will be investigated in more detail inSec. IV C below. At ϕ = 0 . 25, there is again a reentrantbehavior of lamellar states, i.e. nearly planar alignedlayers appear at large shear rates ˙ γτ ≥ . with increase in the shear rate at acertain membrane volume fraction; the mixture changesfrom lamellar state to onion state, and then after goingthrough the coexistence region it enters again into anoriented lamellar state again. Although the onion statewith densely packed MLVs has not been obtained in oursimulations, the reentrant behavior is qualitatively con-sistency with the experimental observations. C. Rolled-up Lamellar Structures 1. Structure Analysis Snapshots of membrane conformations are shown inFig. 7 for ϕ = 0 . γτ = 0 . . . ϕ ≥ . ϕ (cid:38) . 175 and ˙ γτ (cid:38) . n is calculated from the first-order mov-ing least-squares (MLS) method applied to theconfigurations of membrane particles. Using a weightedgyration tensor a αβ = (cid:80) j ( α (cid:48) j − α (cid:48) G )( β (cid:48) j − β (cid:48) G ) w cv ( r ij ),where α, β ∈ { x, y, z } , n is obtained as an eigenvectorcorresponding to the minimum eigenvalue of a αβ , whichtogether with the other two eigenvectors e and e con-stitutes an orthonormal basis. Here, the cut-off function w cv ( r ij ) in Eq. (8) is employed with the same cut-offlengths r (cid:48) cc = 3 . σ . As a quantitative measure for theundulation instability, we employ the orientational orderparameter S z = [2( n · ˆ e z ) − 1] = cos(2 θ ) (18) FIG. 7: Snapshots membrane conformations for volume fraction ϕ = 0 . γτ = 0 . γτ = 0 . γτ = 0 . x ) direction are shown in the upper panels. Corresponding cross-sectional views(with − . σ < x < . σ ) are shown in the lower panels. -3 -3 -3 -2 -1 -2 -1 ϕ (b) ˙ γτ (a) ! S z " ! H " σ FIG. 8: (a) Membrane orientational order parameter (cid:104) S z (cid:105) inEq. (18) for ϕ = 0 . . . 25, and 0 . (cid:104) H (cid:105) for ϕ = 0 . . 25, and0 . γτ . with the normalization for the two-dimensional order ofcylindrical and planar symmetry.The instability can also be characterized by calculat-ing the mean-square curvature. The second-order MLS method provides an estimate of the membrane curva-ture from the particle configurations in the following way.For each particle i , we perform a rotational transforma-tion into the principal coordinate system of the gyrationtensor of the neighbor particles j around i ’s weightedcenter of mass r G by X j Y j Z j = e e n ( r j − r G ) T , (19)and then employ the parabolic fit functionΛ ( r i ) = w (cid:80) j (cid:0) z + h x X j + h y Y j + h xx X j + h yy Y j + h xy X j Y j − Z j (cid:1) w cv ( r ij ) , (20)where the coefficients of the Taylor expansion z , h x , h y , h xx , h yy and h xy are fitting parameters. By a least-squares fit, the estimated value of the meancurvature H = ( C + C ) / i is then obtainedas H = (1 + h x ) h yy + (1 + h y ) h xx − h x h y h xy h x + h y ) / . (21)Figure 8 displays the results for the spatial average (cid:104) S z (cid:105) of the 2D orientational order parameter and forthe mean-square local curvature (cid:104) H (cid:105) . When the mem-branes are aligned with the x − y plane orthogonal tothe shear-gradient direction, S z becomes unity. As themembranes roll up, S z decreases. Here, (cid:104) S z (cid:105) = 0 forperfectly cylindrical state (where ( n · ˆ e x ) = 0 and FIG. 9: Color maps of the structure factor S A (0 , q y , q z ) for volume fraction ϕ = 0 . γτ = 0 . γτ = 0 . γτ = 0 . x ) direction, only data for q x = 0 are shown. ( n · ˆ e y ) = ( n · ˆ e z ) = 1 / (cid:104) S z (cid:105) = − y ) di-rection (where ( n · ˆ e z ) = 0). For all ϕ ≥ . (cid:104) S z (cid:105) and (cid:104) H (cid:105) provide evidence for rolled-upstructures at ˙ γτ = 0 . ϕ = 0 . (cid:104) S z (cid:105) increases again at ˙ γτ = 0 . ϕ = 0 . (cid:104) S z (cid:105) remains small (and (cid:104) H (cid:105) large), whichimplies that rolled-up structures exist also at ˙ γτ = 0 . ϕ = 0 . S A ( q ) at q x =0 are shown for the same set of data as in Fig. 7. Due tothe nearly complete alignment of membranes in the flow( x ) direction, all the structural features are reflected in S A ( q ) in the q y − q z wave-vector plane. Since lamellarlayers are nearly planar at the low shear rate ˙ γτ = 0 . q y , q z ) = (2 πq , γ .In the small-angle neutron and X-ray scattering ex-periments, the scattering beams can be injected from twodirections, either radial or tangential to the shear cell,which correspond to shear-gradient and flow directions,respectively. After the shear flow is applied, after a whilea Bragg peak of the radial beam develops in the vorticitydirection, while the scattering pattern becomes isotropicfor the tangential beam. This suggests 2D-isotropic un-dulations of the lamellar structure perpendicular to theflow direction. Later in time, the radial beam is scatteredisotropically in the tangential direction, indicating theformation of onion structures. The structure factor S ( q )in our simulation (Fig. 9(c)), agrees well with S ( q ) ofthis transient states in these scatting experiments. Mea-surements of solvent diffusion in a Rheo-NMR experi-ment of a non-ionic surfactant system show a diffu-sion anisotropy in the direction of shear flow in the inter-mediate state. These experimental results suggest that the membranes are aligned in the direction of shear flowin the intermediate state of the lamellae-to-onion tran-sition, although more detailed information on the struc-tural membrane arrangement could not be obtained. Therolled-up lamellar structures observed in our simulationsmatch the experimental evidence, so that they are a goodcandidate for the intermediate states. 2. Temporal Evolution of Membrane Structures In simulations of a large system, even if extended overa long time interval, the resultant structure often ex-hibits dependence on the initial conditions – owing to theslow processes involved in the dynamics. Therefore, wecompare here the time evolution at ϕ = 0 . . × MD steps at a small shear rate˙ γτ = 0 . (cid:104) S z (cid:105) forboth initial conditions are shown for ˙ γτ = 0 . . . (cid:104) S z (cid:105) exhibits an overshoot (see Fig. 10(c)), and fi-nally approaches a constant as the structure relaxes intoa (meta)stable state. The overshoot amplitude dependson initial states and random noises.Figs. 10(b) and (c) show the membrane conformationsafter an elapsed time of t = 2 . × τ (for ˙ γτ = 0 . (cid:104) S z (cid:105) in Fig. 10(a). When the random state istaken as initial condition, rolled-up structures are con- ˙ γτ t/τ ! S z " (b) (c)(a) (c)(b) FIG. 10: Comparison of structure formation from a randominitial configuration and from a lamellar state at ϕ = 0 . (cid:104) S z (cid:105) .The dotted and solid lines show the data starting from therandom and lamellar initial states, respectively. Blue, red,and black lines represent the shear rates ˙ γτ = 0 . . . γτ = 0 . 284 at t = 2 . × τ , as they have developedfrom (b) the random and (c) the lamellar initial states, arealso shown. siderably more pronounced; this may be traced back tothe presence of defects in the lamellar structure whichforms at short times t/τ (cid:39) γ with the random ini-tial configuration, less undulations take place with thelamellar initial configuration, as indicated by (cid:104) S z (cid:105) ∼ V. SUMMARY In this paper, we have constructed an explicit-solventmeshless-membrane model for surfactant-water mixtures.The model reproduces properties of an earlier implicit-solvent meshless membrane model, where membranebending rigidity and line tension can be independentlycontrolled to a large extent. The model enables large-scale simulations of structural changes, where dynamicaleffects of hydrodynamic interactions have to be takeninto account. At present, such a large simulation with asmany as one-million particles can be realized by paral-lelized molecular dynamics simulation methods.Our main results concern the effects of shear flow onthe structure formation of membrane ensembles. Variousstructures including vesicles, lamellae, and multi-lamellarstates with nearly cylindrical symmetry have been found,most of which are qualitatively consistent with experi-mental observations of non-ionic surfactant membranesunder shear flow. Especially, a cylindrical instability ofmulti-lamellar membrane is predicted to occur perpen-dicularly to the flow direction. The corresponding scat-tering patterns are in qualitatively agreement with theresults of small angle neutron (and X-ray) scattering ex-periments under shear deformation. The rolled-up lamel-lae are a good candidate for the intermediate structureson the way to the onion state, which are observed in theexperiments.Our simulations do not reproduce onion formation,which is ubiquitously observed in experiments on µ mlength scales. We speculate that it is due to the lim-ited system size, which can be overcome by larger-scalecalculations in the future. On the other hand, in theexperiments, strains larger than ˙ γt (cid:38) are necessaryto reach the onion state, which indicates that very longsimulation runs are required to obtain these states.The control of the physical parameters including theline tension Γ, the bending rigidity κ , and the Gaus-sian modulus ¯ κ is another challenge. While Γ and κ areeasily controlled in our model, ¯ κ is more difficult. TheGaussian modulus ¯ κ might also play an important rolein the structure formation, because it is directly relatedto topological changes happening on the way of onionformation. Appendix A: Solvent-Mediated Forces at HigherSolvent Density In the present simulation model, solvent particles havea similar size as membrane particles. Here, we discussthe finite-size effects of the solvent particles, if muchhigher solvent densities are used than in the presentstudy. When the solvent particles are densely packed,the system approaches a crystallization transition, andlocal crystalline order emerges to bring about interac-tions between closely spaced lamellar layers.As an example, we study the explicit-solvent meshless-0 z xy FIG. 11: Snapshot of the MLV state due to the solvent-mediated attractive interactions for the system S φ =0 . ϕ = 0 . 06 and N = 480 , 000 with different size ratio σ B /σ A = 1 . 2, for shear rate ˙ γτ = 0 . membrane model for higher number density (denoted sys-tem ( S φ = ( σ A N A + σ B N B ) /V = 0 . 72 with thetotal particle number N = N A + N B = 480 000). Theparticle radii of the two components are chosen such as σ B = 1 . σ A , where σ A and σ B denote the radii of mem-brane and solvent components, respectively. The cut-offlengths of the interactions are set to r rep c = 2 . σ A , r ga =1 . σ A , r cc = 3 . σ A , respectively. The repulsive in-verse twelfth-power potential exhibits melting transitionaround the volume fraction around 0.43 (correspondingto φ (cid:39) . Thus, our density φ = 0 . ϕ = N A σ A / ( N A σ A + N B σ B ) = 0 . 06 fora shear rate ˙ γτ = 0 . d (cid:39) . σ B , which seems to arise from the discrete-ness of the interstitial solvent particles – here it shouldbe noted that interlayer distance cannot be smaller than d (cid:39) . σ B because it is inside the range of U rep and U α ,which both act repulsively.To avoid this solvent-mediated force, a longer cut-offrange of the repulsive forces and a lower density are em-ployed in the simulation described in the main text (de-noted system ( S S S 2, we simulate three lamellarlayers in the following way. A periodic simulation boxwith lengths L x = L y = (cid:113) A xy /σ A , L z = V / ( L x L y )is set up so that a tension-less membrane (with particlenumber N A ) = 1600 is put along the xy -plane. Here,the projected areas are set to A xy = 1 . N A for S A xy = 1 . N A for S 2, respectively. With various inter-layer distance d ini , circular membrane disks with particlenumber N A = N A = 400 are put on both sides of the FIG. 12: Equilibrium lamellar distance d lam versus initiallamellar distance d ini between a tension-less planar membraneand two (smaller) membrane discs. For the systems S S 2, lengths are represented in units of σ and σ A , respectively.The dotted line indicates the line d ini = d lam . membrane. Solvent particles are then inserted to fill upthe system; they are relaxed for fixed membrane configu-ration for 10 MD steps. Afterward, a simulation of thefull system is performed for 1 . × MD steps to ob-tain an equilibrium interlayer distance d lam . As shown inFig. 12, while the interlayer distance increases with timein the S d lam = 3 . σ A and 4 . σ A in the S σ B and 4 σ B , respec-tively. Thus, in S 2, the effective attractive potentials leadto the stable multi-lamellar layers due to the attractiveinteractions. Appendix B: Numerics and Parallelization Numerical simulations have been carried out on mas-sively parallel supercomputers. On 256 CPUs of IntelXeon X5570 (2.93GHz) in SGI Altix ICE 8400EX atISSP, it costs 72 hours to perform a run of 1 . × simulation steps for ϕ = 0 . N = 960 , U rep , 30% for U α , around 20% for the construction ofthe buffer, and 20% for the communication between theMPI processes. The system is divided into cubic (or rect-angular) boxes, each of which is calculated by one MPIprocess. We also parallelize calculations of each processwith the use of OpenMP by performing calculations ofdifferent far away pairs at once on different threads. Thecode is optimized to achieve 15% performance comparedwith the theoretical limit on X5570 processors.Each process has separate cell lists and neighbor listsfor solvent and membrane particles. To apply shear, weemploy a simultaneous affine deformation of the total sys-tem, each MPI process box, and each cell for neighborsearch, so that a square is transformed into a parallel-ogram shape consistent with the shear deformation, as1 strain=0.5 convert tovv strain= -0.5 MPI synchronous communication strain=0 MPI MPI MPI MPIMPI MPI MPI MPIMPI MPI MPI MPIMPI MPI MPI MPI strain=0 FIG. 13: Schematic picture in two dimensions for the imple-mentation of Lees-Edwards boundary conditions in a programparallelized with MPI communication. illustrated in Fig. 13. When the strain reaches 0.5, these parallelograms are reflected to make the strain -0.5, seeFig. 13, which basically requires an all-to-all communi-cation of all the position and velocity data. Acknowledgments This work was supported by Grant-in-Aid for YoungScientists 24740285 from JSPS in Japan, Computa-tional Materials Science Initiative (CMSI) from MEXTin Japan, and the European Soft Matter Infrastructureproject ESMI (Grant No. 262348) in the EU. We wouldlike to thank S. Fujii, S. Komura, H. Watanabe, U.Schiller, M. Peltom¨aki, G.A. Vliegenthart, and D. Y. Lufor informative discussions. The numerical calculationswere carried out on SGI Altix ICE 8400EX and NEC SX-9 at ISSP in University of Tokyo (Japan), Fujitsu FX10at Information Technology Center in University of Tokyo(Japan), Hitachi SR16000 at YITP in Kyoto University(Japan), and JUROPA at J¨ulich Supercomputing Centerat Forschungszentrum J¨ulich (Germany). ∗ [email protected] † [email protected] R. G. Larson, The Structure and Rheology of Complex Flu-ids (Oxford University Press, New York, 1999). O. Diat, D. Roux, and F. Nallet, J. Phys. II France , 1427(1993). D. Roux, F. Nallet, and O. Diat, Europhys. Lett. , 53(1993). O. Diat, D. Roux, and F. Nallet, Phys. Rev. E , 3296(1994). L. Courbin, J. P. Delville, J. Rouch, and P. Panizza, Phys.Rev. Lett. , 148305 (2002). S. Koschoreck, S. Fujii, P. Lindner, and W. Richtering,Rheol. Acta. , 231 (2009). Y. Suganuma, M. Imai, T. Kato, U. Olsson, and T. Taka-hashi, Langmuir , 7988 (2010). S. Fujii, D. Mitsumasu, Y. Isono, and W. Richtering, SoftMatter , 5381 (2012). F. Nettesheim, J. Zipfel, U. Olsson, F. Renth, P. Lindner,and W. Richtering, Langmuir , 3603 (2003). M. Ito, Y. Kosaka, Y. Kawabata, and T. Kato, Langmuir , 7400 (2011). B. Medronho, S. Shafaei, R. Szopko, M. G. Miguel, U. Ols-son, and C. Schmidt, Langmuir , 6480 (2008). B. Medronho, C. Schmidt, U. Olsson, and M. G. Miguel,Langmuir , 1477 (2010). B. Medronho, M. Rodrigues, M. G. Miguel, U. Olsson, andC. Schmidt, Langmuir , 11304 (2010). W. Richtering, Curr. Opin. Colloid Interface Sci. , 446(2001). H. Miyazawa and H. Tanaka, Phys. Rev. E , 011513(2007). B. Medronho, S. Fujii, W. Richtering, M. G. Miguel, and U. Olsson, Colloid. Polym. Sci. , 317 (2005). B. Medronho, J. Brown, M. G. Miguel, C. Schmidt, U. Ols-son, and P. Galvosas, Soft Matter , 4938 (2011). S. Fujii, Y. Ishii, S. Komura, and C.-Y. D. Lu, EPL ,64001 (2010). S. Fujii, S. Komura, Y. Ishii, and C.-Y. D. Lu, J. Phys:Condens. Matter , 235105 (2011). A. G. Zilman and R. Granek, Eur. Phys. J. B , 593(1999). S. W. Marlow and P. D. Olmsted, Eur. Phys. J. E , 485(2002). E. van der Lindern, W. T. Hogervorst, and H. N. W.Lekkerkerker, Langmuir , 3127 (1996). E. van der Linden and J. H. M. Dr¨oge, Physica A , 439(1993). C.-Y. D. Lu, Phys. Rev. Lett. , 128304 (2012). H. Guo, K. Kremer, and T. Soddemann, Phys. Rev. E ,061503 (2002). T. Soddemann, G. K. Auernhammer, H. Guo, B. D¨unweg,and K. Kremer, Eur. Phys. J. E , 141 (2004). H. Guo and K. Kremer, J Chem. Phys. , 054902 (2007). O. Henrich, K. Stratford, D. Marenduzzo, P. V. Coveney,and M. E. Cates, Soft Matter , 3817 (2012). H. Noguchi, J. Phys. Soc. Jpn. , 041007 (2009). H. Noguchi and G. Gompper, Phys. Rev. E , 021903(2006). H. Noguchi and G. Gompper, J. Chem. Phys. , 164908(2006). H. Shiba and H. Noguchi, Phys. Rev. E , 031926 (2011). G. Gompper and D. M. Kroll, in Statistical Mechanics ofMembranes and Surfaces , edited by D. R. Nelson, T. Piran,and S. Weinberg (World Scientific, Singapore, 2004), 2nded. D. A. Fedosov, H. Noguchi, and G. Gompper, Biomech.Model. Mechanobiol., submitted. G. Gompper and D. M. Kroll, Phys. Rev. Lett. , 2284(1998). M. Peltom¨aki, G. Gompper, and D. M. Kroll, J. Chem.Phys. , 134708 [1 (2012). J.-M. Drouffe, A. C. Maggs, and S. Leibler, Science ,1353 (1991). M. G. Del P´opolo and P. Ballone, J. Chem. Phys. ,024705 (2008). T. Kohyama, Physica A , 3334 (2009). H. Yuan, C. Huang, and S. Zhang, Soft Matter , 4571(2010). A. Pasqua, L. Maibaum, G. Oster, D. A. Fletcher, andP. L. Geissler, J. Chem. Phys. , 154107 (2010). R. Kapral, Adv. Chem. Phys. , 89 (2008). G. Gompper, T. Ihle, D. M. Kroll, and R. G. Winkler,Adv. Polym. Sci. , 1 (2009). P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys.Lett. , 155 (1992). P. Espa˜nol and P. Warren, Europhys. Lett. , 191 (1995). R. D. Groot and P. B. Warren, J. Chem. Phys. , 4423(1997). H. Noguchi, N. Kikuchi, and G. Gompper, EPL , 10005(2007). H. Noguchi and G. Gompper, EPL , 2007 (2007). C. P. Lowe, Europhys. Lett. , 145 (1999). T. Shardlow, SIAM J. Sci. Comput. , 1267 (2003). S. J. Marrink, A. H. de Vries, and D. P. Tieleman, Biochim.Biophys. Acta , 149 (2009). M. Venturoli, M. M. Sperotto, M. Kranenburg, andB. Smit, Phys. Rep. , 1 (2006). M. Laradji and P. S. Kumar, Phys. Rev. Lett. , 198105(2004). H. Noguchi, Phys. Rev. E , 061919 (2011). W. Helfrich, Z. Naturforsch. C , 693 (1973). T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, andP. Krysl, Comput. Methods Appl. Mech. Eng. , 3(1996). P. Lancaster and K. Salkaskas, Math. Comput. , 141(1981). M. Hu, J. J. Briguglio, and M. Deserno, Biophys. J. ,1403 (2012). W. G. Hoover, M. Ross, K. W. Johnson, D. Henderson,and J. A. Barker, J. Chem. Phys.52