{'fraction_non_alphanumeric': 0.0386762360446571, 'fraction_numerical': 0.02196049564470617, 'mean_word_length': 4.4459662602305, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structure formation of surfactant membranes under shear flow

Shear-flow-induced structure formation in surfactant-water mixtures is investigated numerically using a meshless-membrane model in combination with a particle-based hydrodynamics simulation approach for the solvent. At low shear rates, uni-lamellar vesicles and planar lamellae structures are formed at small and large membrane volume fractions, respectively. At high shear rates, lamellar states exhibit an undulation instability, leading to rolled or cylindrical membrane shapes oriented in the flow direction. The spatial symmetry and structure factor of this rolled state agree with those of intermediate states during lamellar-to-onion transition measured by time-resolved scatting experiments. Structural evolution in time exhibits a moderate dependence on the initial condition.


I. INTRODUCTION
Surfactant molecules in water self-assemble into various structures such as micelles and bilayer membranes, which display a rich variety of rheological properties under flow. 1 Even if a basic structure remains to be a bilayer membrane, its mesoscale structure can assume several different states, such as fluid L α or ripple P β phase. Under shear flow, lamellae can be oriented parallel or perpendicular to the shear-gradient direction. Diat and Roux first discovered 20 years ago closely-packed multilamellar vesicle (MLV) structures, so-called the onion phase, in nonionic surfactant-water mixtures under shear flow. [2][3][4] In the last two decades, this onion structure has been studied experimentally using light, 2,4-8 neutron, 2,4,9 and X-ray scattering, 7,8,10 and also by freeze-fracture electron microscopy, 5 and the rheo-NMR method. [11][12][13] Its rheology has been of large interest. 3,6,8,9,14,15 Typically, a critical shear rateγ c separates the lamellae and onion phases, where the latter phase consists of monodisperse onions containing hundreds of lamellar layers. The onion radius R(γ) is reversible and can be described by a unique decreasing function of the shear rateγ. 16 Time-resolving small-angle neutron scattering experiments have revealed that a two-dimensional intermediate structure is formed during the lamellar-to-onion transition with increasing shear rate. 9 A cylindrical or wavy lamellar structure was speculated to be the transient intermediate structure, but could not be distinguished from the scattering pattern alone. Recent small-angle X-ray scattering experiments with increasing temperature at constant shear rate also indicate a similar pattern around the lamellar-to-onion transition. 10 Thus, there are some experimental evidences of a transient state, but its structure is still under debate. An alternative experimental approach to gain insight into the structural changes is to characterize defects observed in the lamellar state for moderate shear rates, both in surfactant membranes 13,17 and in thermotropic liquid crystals. 18,19 It is also worth mentioning that stable cylindrical structures on a tenmicrometer length scale are observed when strong shear flow is applied in the lamellar-sponge coexistence state. 15 Several theoretical attempts have been made to tackle this complex problem of structural evolution under shear flow, which consider either instability of the lamellar phase due to undulations 20,21 or the break-up of droplets. 22,23 Recently, a "dynamical" free energy of MLV under shear flow has been proposed, 24 which takes into account the slow modes induced by the solvent between the membranes together with their bending and stretching forces. The scaling relations for the MLV size and the terminal shear rate are predicted in agreement with the experiments of Refs. 2,3. In these theories, while the hydrodynamic effects of the solvent are taken into account, the analysis is performed for geometrically simple structures, asp spherical onions or planar lamellae. Thus, the kinetic process of the transformation from the lamellar to the onion phase could not be investigated theoretically so far. In this paper, we study the detailed structural evolution of surfactant membranes under simple shear flow using large-scale particle simulations.
A few simulations have been performed for the formation of lamellar phases in shear flow, while onion formation has never been addressed so far. Oriented lamellae have been obtained in simulations of a coarsegrained molecular model for lipids, [25][26][27] while defect dynamics has been investigated in simulations of a phasefield model of a smectic-A system. 28 Onion and intermediate states have large-scale structures of the order of micrometers, which are beyond typical length scales accessible to molecular dynamics (MD) simulations of coarse-grained surfactant molecules. In our study, we employ a meshless-membrane model, [29][30][31][32] where a membrane particle represents not a surfactant molecule but rather a patch of bilayer membranein order to capture the membrane dynamics on a mi-crometer length scale. This model is well suitable to study membrane dynamics accompanied by topological changes. Alternatively, membranes can be modeled as triangulated surfaces, 29,33,34 which require however discrete bond reconnections to describe topological changes. 33,35,36 After the first meshless-membrane model was proposed in 1991, 37 several meshless-membrane models have been developed. [30][31][32][38][39][40][41] In contrast to other meshless-membrane approaches, our models [30][31][32] are capable of separately controlling the membrane bending rigidity κ and the line tension Γ of membrane edges. Previously, 29,31 we have combined our meshlessmembrane model with multi-particle collision (MPC) dynamics, 42,43 a particle-based hydrodynamic simulation technique. With MPC, the hydrodynamic interactions are properly taken into account, but due to the frictional coupling of membrane and solvent, solvent particles can penetrate through the membrane. We here extend the meshless-membrane model into explicit solvent simulation model, in which the fluid particles interact with each other and with the membrane via short-ranged repulsive potentials, so that the solvent can hardly penetrate the membrane, and simulate it with dissipative particle dynamics (DPD), another hydrodynamics simulation technique.
The simulation model and methods are introduced in Sec. II. Then basic membrane properties, including bending rigidity and line tension, are described in Sec. III. In Sec. IV, structure formation in surfactant-water mixtures under shear flow is investigated for variety of shear rateṡ γ 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 and some perspectives are given in Sec. V.

A. Coarse-Grained Model and Interaction Potentials
To simulate the structure formation in surfactantmembranes systems, we employ a meshless-membrane model with explicit solvent. In this model, two types of particles A and B are employed, which denote membrane and solvent particles, respectively. The number of these particles is N A and N B , which defines the particle density φ = (N A + N B )/V -where V is the volume of the simulation box -and the membrane volume fraction ϕ = N A /(N A + N B ). The particles interact via a potential U , which consists of repulsive, attractive, and curvature interactions, U rep , U att , and U α , respectively, Here, the former sum for repulsive interactions is taken over all pairs of particles, the latter only over the membrane particles.
All neighbor particle pairs interact via the shortranged repulsive potential where r ij is the distance between particles i and j. This potential is cut off at a distance r rep c = 3.2σ. The length σ, representing the particle diameter, is employed as the length unit. The constant B is chosen such as to ensure the continuity of the potential at r = r rep c . The (dimensionless) potential strength is set to c = 4.
To favor the assembly of membrane particles into smoothly curved sheets in three-dimensional (3D) space, membrane particles interact via the additional potentials U att and U α , which have been introduced in the implicitsolvent version of the model previously. [30][31][32] With these potentials, the membranes particles self-assemble into a single-layer sheet, which is a model representation of a bilayer membrane. Here, the attractive interaction is given by which is a function of the local density ρ i of the membrane particles defined by where n = 12, a = ln(2){(s cut /s half ) n −1} with s cut = 2.1 and s half = 1.8 are used. We set ρ * = 6 to study a fluid membrane in an explicit solvent. The curvature potential is introduced to incorporate the membrane bending rigidity. Here, the aplanarity α pl provides a measure for the degree of deviation of membrane particle alignment from a planar reference state. It is defined by where λ 1 ≤ λ 2 ≤ λ 3 are the three eigenvalues of the local gyration tensor a αβ of the membrane near particle i, which is defined as a αβ = j∈A (α j − α G )(β j − β G )w cv (r ij ), with α, β ∈ {x, y, z}. Here, r G = j∈A r j w cv (r ij )/ j∈A w cv (r ij ) is the locally weighted center of mass, and w cv (r ij ) is a truncated Gaussian function which is smoothly cut off at r ij = r cc . The constants are set as n = 12, r ga = 1.5σ, and r cc = 3.2σ. In all our simulations, the volume is chosen such that the number density φ of the particles is constant, For higher solvent densities, the system is closer to the melting point, and strong attractive interaction between membranes sheets have been found (see Appendix A for details). The density φ = 0.64σ −3 is chosen in order to avoid such an attraction.

B. Thermostats
We simulated the membranes in the NVT ensemble under shear flow. To keep the temperature constant, we employ the dissipative particle dynamics (DPD) thermostat, [44][45][46][47][48] in which friction and noise forces are applied to the relative velocities of pairs of neighboring particles. Thus, linear and angular momentum is conserved, which implies that the system shows hydrodynamic behavior on sufficiently large length and time scales. The equation of motion for the ith particle is given by is the Heaviside step function. This type of weight has also been used in the Lowe-Andersen thermostat. 49 The scheme is discretized with Shardlow-S1 splitting method, 50 whose time step ∆t b = 0.2t 0 is different from that for contributions from the molecular interactions, where ∆t = 0.005t 0 (t 0 = m/γ is simulation time unit). In the following, we use γ = √ mk B T /σ, where k B T is the thermal energy unit. Although the DPD thermostat is usually employed for soft interaction potentials, 51-53 it can also be employed for systems with steeper potentials, such as the Weeks-Chandler-Andersen potential. 25,26 Shear flow with velocity v x =γz in the x-direction and gradient in the z-direction is imposed by Lees-Edwards boundary condition. The code is optimized for use on a parallel computer architecture by domain decomposition (see Appendix B).

III. MODEL PROPERTIES
The meshless-membrane model with explicit solvent is expected to exhibit very similar equilibrium properties as the original implicit-solvent version. [30][31][32] We now confirm this dependence of surface tension, line tension, and bending rigidity on the control parameters for the explicit-solvent model.
First, a planar fluctuating membrane is simulated with the particle numbers N A = 1600 (and N = 48 000 or 64 000). For various membrane projected areas A xy = L x L y , the surface tension is investigated. Here, P µν is the pressure tensor given by where , and the sum is taken over all the particles including the solvent component. In calculating P µν , the periodic image µ i + nL µ nearest to the other interacting particles is employed, when the potential interaction crosses one of the periodic boundaries. Figure 1 shows the dependence of γ s on the projected membrane area A xy for k α = 2.5, 5, 10, and 20, with = 4. For γ s 0, the intrinsic area A is larger than A xy due to the membrane undulations; buckling of the membrane occurs for A < A xy (the flat region at γ s < 0 in Fig. 1). 54 For tension-less membranes with N A = 1600, the projected area A 0 xy is given by A 0 xy = a 0 xy N A with a 0 xy = 1.1, 1.27, 1.33, 1.35 for k α = 2.5, 5, 10, 20, respectively. The projected membrane area increases with increasing k α , because both membrane bending fluctuations and protrusions are suppressed at larger k α . Figure 2 shows the bending rigidity κ as a function of k α . Here, the bending rigidity is estimated from the height spectrum 55 of the tension-less membrane (γ s = 0). In calculating |h(q)| 2 , the raw positional data of the membrane particles (r i , i ∈ A) is employed. 30,32 Because of the slow dynamics of long-wavelength height fluctuations, it becomes more time-consuming to obtain precise data than for the implicit-solvent model. 30,32 Therefore, |h(q)| 2 is measured using 16 independent runs for each k α , and the averaged spectrum is then fitted to Eq. (13) in the range q < 1.2σ −1 . As in the implicit model, the bending rigidity is found to be proportional to k α , which demonstrates the controllability of κ in our model. In Fig. 3, the dependence of the line tension Γ is displayed for k α = 5 and 10 for a membrane strip with two edges with lengths equal to L x with N A = 1600. Here, Γ is determined via the relation 32 At a small , Γ is proportional to and almost independent of k α , similarly to the implicit-solvent meshless membrane model 30 . While Γ increases linearly up to = 8 for k α = 10, it levels off and saturates at = 4 for k α = 5. When exceeds the value where Γ saturates -a value which becomes larger for larger k αthe membrane particles prefer to reside at the edges, because the curvature force is not strong enough to avoid aggregation due to the stronger attractive forces. For our simulations of membrane ensembles in shear flow, we choose the model parameters k α = 5 and = 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 ± 0.2, the relaxation time scale of structural transitions can be characterized by where η is the solvent viscosity, and R c is the critical radius of a flat disk. Assuming a flat disk with radius R and a vesicle with the same membrane area, we obtain the corresponding elastic free energies F d = 2πRΓ and F s = 8π(κ +κ/2) 4πκ; thus, a disk exhibits transition to closed vesicle shape if the radius is around With the solvent viscosity η = (2.1 ± 0.1) × m(σt 0 ) −1 obtained from a solvent-only simulation, an estimation of the relaxation time yields τ = 28.4t 0 . In the following, the time will be measured in unit of τ . Since the membrane thickness, typically around 5nm for non-ionic surfactants like polyethylene-glycol-ethers C n E m , corresponds to the size σ of the membrane particles in our simulation, τ is equivalent to about 0.36µs (with the viscosity of water at 300K, η w 0.8mPa · s).

IV. STRUCTURE FORMATION WITH AND WITHOUT SHEAR FLOW
We now employ the meshless-membrane model with explicit solvent to study structure formation in surfactant-water mixtures, both with and without shear flow. In all simulations, the total particle number is fixed at N = N A + N B = 960 000, and thus, the system size is a cubic box with side length L = 114.47σ. Simulations are performed for various membrane volume fraction ϕ = N A /N , with ϕ = 0.0625, 0.125, 0.1875, 0.25, and 0.3125. The dynamical evolution is integrated over a total time interval of 1.2×10 7 MD steps, corresponding to 2.11 × 10 3 τ . All the particles of both species are initially distributed randomly in the simulation box. Averages are calculated over the last 2 × 10 6 steps, where the system is assumed to have reached a stationary state. After briefly explaining the structures obtained by equilibrium simulations without shear in Sec. IV A, we present results for the structure formation in a system under linear shear flow in Secs. IV B and IV C.

A. Mesophases in Thermal Equilibrium
At a low membrane volume fraction ϕ = 0.0625, membrane particles self-assemble into vesicles, each of which is composed of around 100 membrane particles (see Fig. 4(a)). This result is consistent with Eq. (16), which predicts the critical particle number of a vesicle to be N c = πR 2 c /a 0 75. At a higher membrane volume fraction ϕ = 0.125 (N A = 120 000), the membrane surfaces is found to percolate through the whole system via the periodic boundaries. This behavior is not unexpected, because each disk with radius R c covers a region of volume v c = (2R c ) 3 by rotational diffusion, and therefore these disks should overlap and merge when more than n c ∼ V /v c = 1200 vesicles are present, which exceeds the number of vesicles with critical size (N A /N c ).
For ϕ = 0.25 and 0.3125, periodic lamellar states are formed owing to the repulsive interactions between the membranes, as shown in Fig. 5(a). The membranes are curved to fill up space with random orientations. In thermal equilibrium, these lamellar layers would probably have a unique orientation throughout the system; this well-ordered state is difficult to reach in simulations because the structural relaxation time well exceeds the accessible simulation time scale. The 3D structure factor of the membrane density is calculated as where δn A (r) = i∈A σ 3 δ(r − r i ) − φ is the local deviation of the membrane particle density from its average. Figure 5(b) shows S A (q) for ϕ = 0.3125. The scattering intensity is spherically symmetric, as demonstrated by a two-dimensional (2D) color-map for q z = 0. Peaks arising from the inter-lamellar distance are observed at |q| = q 1 = 1.74σ −1 and q 2 = 3.49σ −1 = 2q 1 with heights 9.8 and 0.78, respectively. The former corresponds to length of L = 2π/q 1 = 3.61σ, which provides a precise estimate of the interlayer distance.

B. Dynamic Phase Diagram under Shear Flow
In shear flow, the vesicle and lamellar states depend on the concentration ϕ and the shear rateγ, as displayed in Fig. 6. At ϕ = 0.0625, assemblies of uni-lamellar vesicles are observed forγτ < 0.02. At a higher shear rate, the membrane particles tend to assemble into plate-like membrane disks, which then align parallel to the shear flow direction, as demonstrated by the comparison of snapshots in Figs. 4 (a) and (b). At ϕ = 0.125, lamellar layers with non-uniform lamellar distances are observed for   Fig. 4(c). There is a larger variety of phases at ϕ = 0.1875; at low shear ratesγτ < 0.1, the system is a mixture of lamellae and cylinders as shown in Fig. 4(d); atγτ = 0.142, the lamellae are rolled up collectively, see Fig. 4(e); finally, at large shear ratesγτ 0.3, the system exhibits a reentrant lamellar state, as shown in Fig. 4(f).
At higher ϕ, the lamellar states align perpendicularly to the shear-gradient (z) direction in the regime of low shear rates (γτ < 0.1). At largerγ, they exhibit an instability to a rolled-up shape whose axis is parallel to the flow direction, which will be investigated in more detail in Sec. IV C below. At ϕ = 0.25, there is again a reentrant behavior of lamellar states, i.e. nearly planar aligned layers appear at large shear ratesγτ ≥ 0.284. Note that experimental phase diagrams often exhibit reentrant behaviors 2,3 with increase in the shear rate at a certain membrane volume fraction; the mixture changes from lamellar state to onion state, and then after going through the coexistence region it enters again into an oriented lamellar state again. Although the onion state with densely packed MLVs has not been obtained in our simulations, the reentrant behavior is qualitatively consistency with the experimental observations. The structural changes accompanying this instability can be characterized by the average orientation and the mean-square local curvature of the membrane. The normal unit vector n is calculated from the first-order moving least-squares (MLS) method 30,56,57 applied to the configurations of membrane particles. Using a weighted gyration tensor a αβ = j (α j − α G )(β j − β G )w cv (r ij ), where α, β ∈ {x, y, z}, n is obtained as an eigenvector corresponding to the minimum eigenvalue of a αβ , which together with the other two eigenvectors e 1 and e 2 constitutes an orthonormal basis. Here, the cut-off function w cv (r ij ) in Eq. (8) is employed with the same cut-off lengths r cc = 3.2σ. As a quantitative measure for the undulation instability, we employ the orientational order parameter with the normalization for the two-dimensional order of cylindrical and planar symmetry. The instability can also be characterized by calculating the mean-square curvature. The second-order MLS method 30 provides an estimate of the membrane curvature from the particle configurations in the following way. For each particle i, we perform a rotational transformation into the principal coordinate system of the gyration tensor of the neighbor particles j around i's weighted center of mass r G by and then employ the parabolic fit function where the coefficients of the Taylor expansion z 0 , h x , h y , h xx , h yy and h xy are fitting parameters. 30 By a least-squares fit, the estimated value of the mean curvature H = (C 1 + C 2 )/2 for particle i is then obtained as In Fig. 9, results for the structure factor S A (q) at q x = 0 are shown for the same set of data as in Fig. 7. Due to the 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 lamellar layers are nearly planar at the low shear rateγτ = 0.0284, a sharp peak is observed around (q y , q z ) = (2πq 1 , 0). Because of the undulation instability, the pattern becomes more circular at higherγ.
In the small-angle neutron 9 and X-ray 10 scattering experiments, the scattering beams can be injected from two directions, 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 while a Bragg peak of the radial beam develops in the vorticity direction, while the scattering pattern becomes isotropic for the tangential beam. This suggests 2D-isotropic undulations of the lamellar structure perpendicular to the flow direction. Later in time, the radial beam is scattered isotropically in the tangential direction, indicating the formation of onion structures. The structure factor S(q) in our simulation (Fig. 9(c)), agrees well with S(q) of this transient states in these scatting experiments. Measurements of solvent diffusion in a Rheo-NMR experiment of a non-ionic surfactant system 11 show a diffusion anisotropy in the direction of shear flow in the intermediate state. These experimental results suggest that the membranes are aligned in the direction of shear flow in the intermediate state of the lamellae-to-onion transition, although more detailed information on the structural membrane arrangement could not be obtained. The rolled-up lamellar structures observed in our simulations match the experimental evidence, so that they are a good candidate for the intermediate states.

Temporal Evolution of Membrane Structures
In simulations of a large system, even if extended over a long time interval, the resultant structure often exhibits dependence on the initial conditions -owing to the slow processes involved in the dynamics. Therefore, we compare here the time evolution at ϕ = 0.3125 starting from both, a lamellar state and a random distribution of membrane particles (the latter corresponding to the simulations described in the previous subsections). For an initial lamellar state, the final configuration of a simulation run with 1.2 × 10 7 MD steps at a small shear ratė γτ = 0.0284 (see Fig. 7(a)), is employed.
In Fig. 10(a), the orientation order parameter S z for both initial conditions are shown forγτ = 0.568, 0.284, and 0.142. For the case of a random initial distribution, small disks merge into randomly oriented surfaces, which then align in the shear flow to become a nearly perfect lamellar stack with some defects. Afterwards, lamellae roll up into slightly larger rolls. During rolling up, S z exhibits an overshoot (see Fig. 10(c)), and finally approaches a constant as the structure relaxes into a (meta)stable state. The overshoot amplitude depends on initial states and random noises.
Figs. 10(b) and (c) show the membrane conformations after an elapsed time of t = 2.11 × 10 3 τ (forγτ = 0.284) for the two types of initial conditions, and explain the origin of the substantially different values of the order parameter S z in Fig. 10(a). When the random state is taken as initial condition, rolled-up structures are con- siderably more pronounced; this may be traced back to the presence of defects in the lamellar structure which forms at short times t/τ 100.
For the case of a lamellar initial configuration, the undulation instability becomes more conspicuous when the applied shear flow is stronger. Moreover, while strong undulations are observed at lowγ with the random initial configuration, less undulations take place with the lamellar initial configuration, as indicated by S z ∼ 1. Thus, the initial conditions play an important role in the selection of the transition path and structure formation. It may depend not only on the shear rate and relaxation time but also on the distribution of structural defects in the lamellar layers. Thus, more systematic studies are required in the future to clarify the hysteresis of these systems.

V. SUMMARY
In this paper, we have constructed an explicit-solvent meshless-membrane model for surfactant-water mixtures. The model reproduces properties of an earlier implicitsolvent meshless membrane model, where membrane bending rigidity and line tension can be independently controlled to a large extent. The model enables largescale simulations of structural changes, where dynamical effects of hydrodynamic interactions have to be taken into account. At present, such a large simulation with as many as one-million particles can be realized by parallelized molecular dynamics simulation methods.
Our main results concern the effects of shear flow on the structure formation of membrane ensembles. Various structures including vesicles, lamellae, and multi-lamellar states with nearly cylindrical symmetry have been found, most of which are qualitatively consistent with experimental observations of non-ionic surfactant membranes under shear flow. Especially, a cylindrical instability of multi-lamellar membrane is predicted to occur perpendicularly to the flow direction. The corresponding scattering patterns are in qualitatively agreement with the results of small angle neutron (and X-ray) scattering experiments under shear deformation. The rolled-up lamellae are a good candidate for the intermediate structures on the way to the onion state, which are observed in the experiments.
Our simulations do not reproduce onion formation, which is ubiquitously observed in experiments on µm length scales. We speculate that it is due to the limited system size, which can be overcome by larger-scale calculations in the future. On the other hand, in the experiments, strains larger thanγt 10 4 are necessary to reach the onion state, which indicates that very long simulation runs are required to obtain these states.
The control of the physical parameters including the line tension Γ, the bending rigidity κ, and the Gaussian modulusκ is another challenge. While Γ and κ are easily controlled in our model,κ is more difficult. The Gaussian modulusκ might also play an important role in the structure formation, because it is directly related to topological changes happening on the way of onion formation. 58

Appendix A: Solvent-Mediated Forces at Higher Solvent Density
In the present simulation model, solvent particles have a similar size as membrane particles. Here, we discuss the finite-size effects of the solvent particles, if much higher solvent densities are used than in the present study. When the solvent particles are densely packed, the system approaches a crystallization transition, and local crystalline order emerges to bring about interactions between closely spaced lamellar layers.
As an example, we study the explicit-solvent meshless- membrane model for higher number density (denoted system (S2), where φ = (σ 3 A N A + σ 3 B N B )/V = 0.72 with the total particle number N = N A + N B = 480 000). The particle radii of the two components are chosen such as σ B = 1.2σ A , where σ A and σ B denote the radii of membrane and solvent components, respectively. The cut-off lengths of the interactions are set to r rep c = 2.7σ A , r ga = 1.5σ A , r cc = 3.0σ A , respectively. The repulsive inverse twelfth-power potential exhibits melting transition around the volume fraction around 0.43 (corresponding to φ 0.8 in our definition). 59 Thus, our density φ = 0.72 is close to the crystallization line.
As an example, Fig. 11 shows a snapshot for membrane volume fraction ϕ = N A σ 3 A /(N A σ 3 A + N B σ 3 B ) = 0.06 for a shear rateγτ = 0.0142. After an initial relaxation from a random configuration of the membrane particles, the membrane layers assemble into stacks, and then form MLVs after gathering a certain amount of the membrane patches. In the snapshot, the layers are stacked with distances d 3.0σ B , which seems to arise from the discreteness of the interstitial solvent particles -here it should be noted that interlayer distance cannot be smaller than d 2.0σ 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-off range of the repulsive forces and a lower density are employed in the simulation described in the main text (denoted system (S1)). To demonstrate the difference between the systems S1 and S2, we simulate three lamellar layers in the following way. A periodic simulation box with lengths L x = L y = A 0 xy /σ A , L z = V /(L x L y ) is set up so that a tension-less membrane (with particle number N A0 ) = 1600 is put along the xy-plane. Here, the projected areas are set to A 0 xy = 1.267N A0 for S1 and A 0 xy = 1.335N A0 for S2, respectively. With various interlayer distance d ini , circular membrane disks with particle number N A1 = N A2 = 400 are put on both sides of the membrane. Solvent particles are then inserted to fill up the system; they are relaxed for fixed membrane configuration for 10 5 MD steps. Afterward, a simulation of the full system is performed for 1.5 × 10 6 MD steps to obtain an equilibrium interlayer distance d lam . As shown in Fig. 12, while the interlayer distance increases with time in the S1 case, there are stable lamellar distances around d lam = 3.5σ A and 4.7σ A in the S2 case with higher solvent density, which correspond to 3σ B and 4σ B , respectively. Thus, in S2, the effective attractive potentials lead to the stable multi-lamellar layers due to the attractive interactions.

Appendix B: Numerics and Parallelization
Numerical simulations have been carried out on massively parallel supercomputers. On 256 CPUs of Intel Xeon X5570 (2.93GHz) in SGI Altix ICE 8400EX at ISSP, it costs 72 hours to perform a run of 1.2 × 10 7 simulation steps for ϕ = 0.3125 and N = 960, 000. Here, about 10% of the total time is for the force calculation of U rep , 30% for U α , around 20% for the construction of the buffer, and 20% for the communication between the MPI processes. The system is divided into cubic (or rectangular) boxes, each of which is calculated by one MPI process. We also parallelize calculations of each process with the use of OpenMP by performing calculations of different far away pairs at once on different threads. The code is optimized to achieve 15% performance compared with the theoretical limit on X5570 processors.
Each process has separate cell lists and neighbor lists for solvent and membrane particles. To apply shear, we employ a simultaneous affine deformation of the total system, each MPI process box, and each cell for neighbor search, so that a square is transformed into a parallelogram shape consistent with the shear deformation, as illustrated in Fig. 13. When the strain reaches 0.5, these parallelograms are reflected to make the strain -0.5, see Fig. 13, which basically requires an all-to-all communication of all the position and velocity data.