Conformational changes of polyelectrolyte chains in solvent mixtures †

We numerically investigate the behaviors of polyelectrolyte chains in solvent mixtures, taking into account the effects of the concentration inhomogeneity and the degree of the ionization. When changing the interaction parameters between the solvent components, we found a first order transition of the polymer conformation. In the mixing state far from the coexistence curve, the polymers behave as semi-flexible chains. In the phase-separated state, on the other hand, they show compact conformations included in the droplets. As the interaction parameters of the mixture are increased, an inhomogeneous concentration field develops around the polymer and induces critical Casimir attractive interactions among the monomers. The competition between the electrostatic interactions and the critical Casimir ones gives rise to drastic changes in the conformation.


Introduction
Polyelectrolytes play important roles in many systems from biology to industry.2][3][4][5][6] In comparison to electrically neutral polymers, polyelectrolytes are more rigid because of long-range electrostatic interactions among the monomers. 7,8They tend to behave like rod or semi-flexible polymers in good or theta solvents.In poor solvent conditions, on the other hand, effective attractive interactions work among the monomers.][11][12] One easily accessible way to control the interactions is to add salts.When salt is added to a polyelectrolyte solution, the electrostatic interactions among the monomers are screened and weakened.Then, the polyelectrolytes adopt compact conformations, if the molecular attractive interactions work among the monomers.4][15][16][17] They are strongly attached to the polymers, and some of them bridge the polymers, so that a small amount of the multivalent salt induces large conformational changes, compared to monovalent salts.
To control the interactions, solvent mixtures are also often useful. 13,14The electrostatic interaction is inversely proportional to the dielectric constant, dependent on the mixing fraction.The electrostatic interactions become stronger when a solvent with a small dielectric constant is mixed.The solvent mixtures also change the electrostatic interactions in a different way.9][20][21][22][23] Usually, polyelectrolytes are highly charged in water-like solutions, while they are poorly charged in less polar liquids.5][26][27][28][29] It is not so simple to use solvent mixtures to control the interactions in polyelectrolytes, even if the mixtures are homogeneously mixed.
1][32][33][34] The dimensions of the polymer are changed when the solvent mixture approaches its critical point.There, the concentration inhomogeneity in the solvent mixture induces effective attractive interactions among the monomers.On the other hand, there are only a few studies which consider the effects of the concentration inhomogeneity on the polyelectrolytes. 27,35To our knowledge, there is no study which considers both the contributions from the electrostatic interactions and the concentration inhomogeneity to the polymer conformation, in spite of their importance.Solvent mixtures have been importantly utilized in applications of polyelectrolytes in biology.For example, ethanol is added to aqueous solutions of polyelectrolytes, such as DNA and proteins, in order to purify and concentrate them. 14DNA collapses into a compact state in aqueous solutions of poly(ethylene glycol). 13,14,24,27,36So, it is very important to understand the behaviors of polyelectrolytes in solvent mixtures.
8][39] In particular, we focus on the effects of the ionization degree and the attractive interactions, caused by the concentration inhomogeneity, on the conformation of the polymer chain.Here, we follow previous studies on charged colloids in binary mixtures. 40It was indicated that the inhomogeneous concentration profiles near the colloidal surfaces induce attractive interactions among the particles.The ionization degree and the affinity between the colloid and the solvent depend on the interaction parameters between the solvent components.

Free energy
We consider polyelectrolyte chains immersed in a mixture of water-like (W) and organic (O) solvents, which are both good solvents for the polymer.The W-solvent is more polar than the O-solvent.Salt was added to the mixture and was dissociated into monovalent ions.We describe the system with five coarsegrained variables.The polymers are described with y (see below), and the compositions of the W-and O-solvents are (1 À y)f and (1 À y)(1 À f), respectively, where 0 r f r 1.The cation and anion fractions are given by c c and c a , respectively.
Polyelectrolyte chains are described by a particle-based model, namely, the bead-spring model.Each chain consists of N beads and the i-th bead is expressed by a shape function with a smooth profile y i (r) in a lattice space as 37,41 y i ðrÞ ¼ y 0 h i ðrÞ; where r is the lattice coordinate, a is the bead radius and d is the thickness of the bead interface.R i is the center of the mass of the i-th bead and it is defined in an off-lattice space.In the limit of a c d, y i (r) approaches y 0 inside the bead, while it is zero in the exterior.Each bead corresponds to the Kuhn monomer 4 consisting of n repeat units, and each repeat unit has one ionizable group.y 0 is the density of the repeat units inside the bead and is given by y 0 = nv 0 /V i , where v 0 is the volume of the repeat unit (and solvent molecules) and we assume v 0 = d 3 .
is the effective volume of the i-th bead.
The distribution of the beads is defined as where N p is the number of chains in the simulation boxes.The free energy functional consists of the following five terms, The first term in the right hand side of eqn ( 2) is related to the bead-spring model.It is given by where K is the constant of the spring connecting the adjacent beads.b is the natural spring length corresponding to the Kuhn length 4 and we set b = nd.R i, j = |R i À R j | is the distance between the i-th and j-th beads.U(R i, j ) is a direct interaction between the beads.We assume short-ranged repulsive interactions of U(R i, j ) since both of the solvents are good solvents.It is given by U 0 (>0) is the strength of the interaction.With appropriate K, b and U 0 , the polymer behaves as a self-avoiding flexible chain.Fig. 1(a) shows a schematic picture of our bead-spring model.The circles with solid lines indicate the beads, which correspond to the Kuhn monomers.The small circles with broken lines represent the repeat units.In Fig. 1(a), n = 4 repeat units compose a bead.We note that the kinetics of each repeat unit is not considered in our simulations.The size of a cell grid corresponds to those of the solvent molecules and the repeat units (d).
Since each bead contains n ionizable groups, we treat the ionization degree with a continuous variable a, ranging from zero to unity.Depending on the local environment, it obeys the free energy of the ion dissociation given by 23,40  where b = (k B T) À1 is the inverse of the temperature T with k B being the Boltzmann constant.D 0 and D 1 control the dissociation.If D 1 > 0, the ionizable groups are dissociated more in the W-solvent.
Upon the ionization, the polymer is dissociated into a negatively charged polymer and mobile cations.We assume that the counterions and the cations from the salt are the same species, so that they are described by a single variable c c .On the other hand, the small anions originate only from the salt.The conservation of the ions and the charge neutrality conditions give where c s is the salt concentration and V t is the total system volume.
We assume the free energy of the ions as 40,[42][43][44] ignoring the direct interactions between the polymer and the ions, and the volume effect of the ions. 45g k represents the interaction parameters between the W-solvent and the ions (k = c, a).The affinity of the ion to a solvent is often discussed with its solvation energy in the Born model. 43Most ions prefer to be dissolved in water than organic solvents because the dielectric constant of water is usually larger than those of organic solvents.The Born model is sometimes not capable of explaining some properties of the actual ion affinity. 42,46,47For example, the correlation of the fluctuating ion concentration modifies the solvation energy and makes it dependent on the ion concentration. 47When the water fraction is sufficiently small, the formation of a solvation shell of the water molecules gives rise to a non-linear dependence of the solvation energy on the water composition. 42Furthermore, it is known that some organic ions show hydrophobicity. 42,43Although the solvation energy of the ions is difficult to model, we employ constant solvation energies g k in this model, for simplicity.3][44] We need to introduce more specific interactions of ions with solvents when our model is applied to specific compounds of polyelectrolyte solutions.F ele is the electrostatic energy given by where F is the electrostatic potential and e is the permittivity.
We assume e depends linearly on f and the polymer fraction as with e being the average of the permittivity.e f and e y correspond to the dielectric constants of the W-solvent and the polymer.The electrostatic potential is obtained by solving the Poisson equation, where e is the elementary electric charge.Z c and Z a are the valencies of the cations and anions.It is known that the polymer conformation is drastically changed when multivalent ions are added.In this study, however, we consider only monovalent ions of Z c = ÀZ a = 1, since our coarse-grained description of the ions is insufficient to consider some of the effects of multivalent ions, such as the ion-bridging effect. 15In order to consider them, the particle-based description of the ions should be considered.The right hand side of eqn (10) is the local electric charge density. 23Its last term Àay represents the electric charge of the polyelectrolytes, which are negatively charged.
The solvent free energy is given in terms of f by 48 where w is the interaction between the two solvents, and g p represents the affinity between the polymer and the W-solvent.
C is a constant coefficient of the gradient energy of f.The phase diagram is indicated in Fig. 1(b).

Kinetics
We treat only f(r) and {R i } (or y(r)) as slow variables.In other words, the ion distributions and the ionization degree are calculated by minimizing F with respect to c k and a at any time t.
We readily obtained where l k is a Lagrange multiplier for eqn (6).To determine l k (k = c, a), we solve eqn ( 6) numerically (see Appendix).Substituting the solutions of eqn (13) into eqn (10), we get the Poisson-Boltzmann equation.When the polymer chain is stretched as a rod and it is highly charged, counterion condensation will occur. 1,4or the ionization degree, we also obtained We note that this ionization degree depends on not only the concentration f, but also the local electrostatic potential F and the salt concentration c s via l c .The concentration field obeys 48 The first term of the right hand side stands for the convection by the flow u(r), and D f is the diffusion constant of the solvents.The hydrodynamic flow is obtained by solving where r is the mass density.
$ P is the mechanical stress tensor and is given by , where Z(r) is the viscosity.In fluid particle dynamics, it depends on the polymer distribution as ZðrÞ ¼ Z 0 þ DZ P i h i ðrÞ. 37,38Z 0 is the solvent viscosity and Z 0 + DZ is the viscosity of the polymer bead.p is the pressure, which imposes the incompressible condition This journal is © The Royal Society of Chemistry 2016 rÁu = 0.
$ s R represents the random stress noise, which satisfies the fluctuation-dissipation relation at T. 49 The first term is due to the concentration inhomogeneity and includes the effect of the interface tension. 48The second term of eqn ( 16) represents the force acting on the beads.F i is given by F i = ÀqF/qR i . 41,50he kinetics of the beads are solved in the off-lattice space by means of a fluid particle dynamics method.The bead is transported by the hydrodynamic flow averaged inside the bead as 37 Here we note the roles of the random noises in our simulation.Our model is a hybrid one, employing the field theoretic and particle-based descriptions.The ion distributions and the ionization degree are determined by minimizing the free energy adiabatically.Also, we do not impose the random noises to f, c k and a.Thus, our model is not beyond the mean-field description when we concern the resulting interactions among the beads.Although the correlations among the electric charges would lead to interactions in the polymers, 2,51 such charge correlation effects are not considered in this work.
On the other hand, we consider the conformation of the polymer chains with the particle-based description and we introduce the random noises for the chain motions via hydrodynamic flow in eqn (16).Then, the kinetics of the polymer conformations are given by a Langevin-type equation.As we will see below, the conformational change occurs like a first-order transition, so the random noises for the particle motions play important roles.Such explicit conformational changes cannot be treated in the field theoretic description of the polymer chains. 52,53

Simulation parameters
In numerical simulations, we assume the molecular size d(=v 0 1/3 ) = 3 Å and the bead radius a = 2d.Thus, each bead contains n(=2a/d) = 4 repeat units and the repeat unit density is set to y 0 (Dnv 0 /V i ) = 0.12.In F pol , we set bKd 2 = 10 and bU 0 = 1, with which we ensured |R i,i+1 À b| { d.The average Bjerrum length is c B = e 2 /(4p ek B T) = 3d, which corresponds to e D 62 at T = 300 K.The ionization parameters are set D 0 = 2 and D 1 = 10, which give a 0 D 0.999 and a 0 D 0.119 for f = 1 and f = 0, respectively.Here, a 0 is obtained from dF dis /da = 0 as a 0 = {1 + e D0ÀD1f } À1 .In F sol , we use C = d 2 .A mean field theory gives the interface tension as g = k B T|w À 2| 3/2 /(2d 2 ), 48 which estimates g = 4.76 mN m À1 at w = 2.35.The ion interaction parameters are g c = g a = 2, both of which correspond to the Gibbs energy of transfer of 4.99 kJ mol À1 .The ions prefer the W-solvent (hydrophilic).For simplicity, we also set g p = 0, with which the polymer backbone is neutral in the affinity to the solvent components.In this article, we fix the spatial average of the W-solvent composition at f = 0.2.This choice of f is because we expect that attractive interactions caused by the inhomogeneous concentration field are strengthened when the fraction of the W-solvent is small, as in colloidal systems. 54,55The binodal and spinodal points for f = 0.2 are w t = 2.31 and w s = 3.125, respectively.We add the salts, the concentrations of which are set to c s = 1 Â 10 À5 and 4 Â 10 À5 .They correspond to the Debye parameter ka = (8pc B c s /v 0 ) 1/2 a E 0.055 and 0.11, respectively.
It is known that the coexistence curve of the binary mixture shifts when impurities such as salts are dissolved.The degree of the change depends on the affinity of the impurity and the concentration. 43For our parameters of c s and g k , the shift of the phase separation point is negligibly small, so that we do not consider the change of the phase separation point in this study.The volume fraction of the polymer beads is 7.9 Â 10 À5 .We set the salt concentration to c s = 4 Â 10 À5 .The colors of the beads represent the degree of ionization.The interfaces between the phase-separated domains are also drawn with isosurfaces of f = 0.5.In the mixed state, far from the coexistence curve at w = 1.8 as in Fig. 1(b), the polymers behave as flexible chains.Here sharp interfaces are not formed.At w = 2.35, on the other hand, they are compartmented in the droplets of the phase-separated domains.Here, the relative fraction of the W-solvent rich phase is 0.022 and this minority phase forms the droplets.Since the droplet size is smaller than the contour length of the polymers, the polymer chains are strongly distorted by the interface tension in the droplets.

Results and discussion
Since the polymer volume fraction is small and the relative fraction of the W-solvent rich phase is also small, each droplet contains a single chain and the droplets are dispersed well.But, some of the polymer chains can occasionally be contained in the same droplet.In order to see detailed behaviors of the conformational changes, it is better to avoid many body effects.Then, we consider a single polyelectrolyte chain in a cuboid box, the size of which is V t = 256d Â (128d). 2 The polymer is composed of N p = 40 beads and its volume fraction is 3.2 Â 10 À4 .Since the periodic boundary conditions are employed, the simulated polymer corresponds to one of the chains in a dilute polymer solution, not an isolated one in an infinite cell.
Fig. 3(a) and (b) show their snapshots, for which w is changed between 2.1 and 2.35 (ESI †).The salt concentration is c s = 1 Â 10 À5 .Assuming a lower critical solution temperature (LCST)-type mixture, we consider two quenching processes; heating from a mixed state with increasing w in Fig. 3(a) and cooling from a phase-separated state with decreasing w in Fig. 3(b).In Fig. 4, we plot their gyration radii R g with respect to w.This is calculated as and averaged in a time interval (282a 2 /(6D p ) r t r 564a 2 /(6D p )), where D p is the diffusion constant of the particles.The simulation of a neutral polymer of N = 40 in a solvent of f = 0 estimates its gyration radius as R g E 4.85b.R g of a neutral polymer is indicated by the broken line in Fig. 4. The gyration radii in the solvent mixture with c s = 4 Â 10 À5 are also plotted.
In Fig. 3(a) and (b), we observe drastic changes in the polymer conformations.As w is increased, an expanded chain collapses into a compact state and vice versa.In both cases, the polymer behaves as a flexible chain at w = 2.1(ow t ).As indicated in Fig. 4, it swells larger than the neutral polymer of the same N because of the electrostatic interactions.At w = 2.35(>w t ), on the other hand, the bulk mixture phase-separates and then the polyelectrolyte is covered by the minority W-rich droplet.Since the droplet radius is smaller than the gyration radius, the contributions of the electrostatic interactions and the interface tension compete in determining the droplet shape and the polymer conformation.In the cases of w = 2.35 in Fig. 3(a) and (b), the interface tension is so large that the droplet remains spherical.In the intermediate range of w (2.2 r w r 2.25), the structures depend on the quenching path.This difference of the conformations is also seen in the hysteresis loop of R g in Fig. 4, suggesting the nature of a first order transition. 27,28It is interesting that the drastic conformational changes occur in the one-phase state of the bulk mixture (w o w t ).
In the case of c s = 4 Â 10 À5 , the polymer chains tend to collapse at smaller w, in comparison with that of c s = 1 Â 10 À5 .This is because the salt not only weakens the electrostatic repulsions, but also enhances the affinity of the charged polymers to the W-solvent. 44As we will see below, the affinity to one of the solvent components gives rise to the attractive interactions among the Kuhn monomers.

Ionization and adsorption
Fig. 5(a) plots the degrees of ionization as functions of w in the two processes of Fig. 3(a) and (b).They are averaged inside the beads as a ¼ Ð draðrÞyðrÞ Ð dryðrÞ.At low w, the ionization degrees are small and then they increase with w.The increase of a accompanies that of the electric charges.The slight increase of R g for w r 2.25 in the heating process (see Fig. 4) reflects the resultant enhancement of the electrostatic repulsions.We attribute the change of a to that of the concentration field.The effective affinity between the polymer and the W-solvent is given by Even though g p = 0 as in the present simulations, the polymer tends to adsorb the W-solvent owing to D 1 a in g ˜p.In Fig. 5(b), we plot the adsorption amount G, which is defined as f drfðrÞ, with w.As w is increased, the polymer beads adsorb the W-solvent more, since the susceptibility of the concentration field to the external field (here g ˜p) is increased.Since the local composition of the W-solvent is enriched around the beads, the ionizations of the beads occur more widely.Furthermore, the progress of the ionization makes the polymer more preferred to the W-solvent. 23hus, the ionization and adsorption proceed with w cooperatively.Fig. 5(b) also shows that the adsorption amount in the cooling process is larger than that in the heating one.It indicates that the compact conformations can adsorb the W-solvent more than the expanded conformations.The ionization degree in the cooling process is large compared to that in the heating process, too.This is because the larger amount of the W-solvent wets the chain in the cooling process.

Pairwise interaction between two beads
The inhomogeneous concentration field induces attractive interactions among the beads.For colloidal systems, they are sometimes called critical Casimir interactions in mixed solvents, or capillary ones in phase-separated solvents.Here, we consider a bead pair (N = 2) to see the pairwise interaction more clearly.In Fig. 6, we plot the force F 1,2 acting on the bead pair with respect to the particle separation R 1,2 .It is estimated by the spring force F 1,2 = K(R 1,2 À b) by changing the natural spring length b.In Fig. 6, we set $ s R ¼ 0. In the equilibrium states, this force balances with the sum of the other interactions, such that this spring force gives an estimation of the force between the nonconnected beads.Its positive and negative values mean repulsive and attractive interactions, respectively.The salt concentration is set to c s = 1 Â 10 À5 .Fig. 6 shows that the force consists of shortranged attraction and a long-ranged repulsion as reported for colloidal systems. 40The former attraction is due to the inhomogeneous concentration field and the latter repulsion comes from the electrostatic interaction.As w is increased, the repulsive interaction becomes strong for large R 1,2 .This is because the ionization degree is increased with w as in Fig. 5(a).Also, the range showing the attractive force expands, and the corresponding force becomes negatively large with increasing w.The critical Casimir interaction (or the capillary interaction) works between the beads when the enriched profiles of the W-solvent around the beads overlap each other.As in Fig. 5(b), the adsorption layer of the W-solvent develops with w, so that the inhomogeneous concentration profiles can overlap easily when w is large, even for large separations.The changes of the force curves leading to the conformational changes seen in Fig. 3.However, the force curves themselves cannot explain the hysteresis behaviors of R g in Fig. 4. The combination of the total interactions among the beads and the conformation entropy would give rise to the first order transition suggested in Fig. 3 and 4. 27,28 4 Conclusions and future works In this article, we numerically studied the conformations of polyelectrolytes in solvent mixtures.We take into account the electrostatic interaction, the effective interaction induced by the concentration inhomogeneity and the ionization degree, simultaneously.We observed drastic changes in the conformation and a hysteresis loop of R g when changing the interaction parameters between the solvent components.Approaching the phase separation point, the inhomogeneous concentration field develops around the beads and it induces attractive interactions among the beads.With these attractive interactions, the polymer collapses into a compact conformation against the electrostatic repulsions.We also found the ionization degree and the inhomogeneous concentration field develop with the interaction parameters in a cooperative way.
Our simulations demonstrate the collapsed state of polyelectrolytes in the mixed state of a solvent mixture.They may  suggest an aggregation and precipitation of polyelectrolyte chains in solvent mixtures, which do not show phase separation, as observed in the so-called ethanol precipitation of DNA.Addition of an organic liquid reduces the dielectric constant of the solvent and strengthens the electrostatic interactions among the ions and polymers.The stable ionic bonds between the polymers and the counterions are formed, and then the polymers become uncharged.The uncharged chains will aggregate and precipitate.Usually, the solvent mixture is considered to be well mixed, since it does not show the phase separation at any accessible temperature.
Addition of salt into a mixed solvent sometimes leads to phase separation, that is, salting out.We consider that the phase separation might be induced locally around the polyelectrolytes and the interface tensions work to aggregate them, even if the mixture appears to be mixed macroscopically.In order to confirm this mechanism in actual polyelectrolyte solutions, some experimental observations are desired.However, it would not be easy to measure the interface tension of such droplets of a molecular scale size.If we can measure the local contents of the W-solvent and the ions in the aggregates, and they are larger than those in the supernatant, they would support our mechanism.
We have many parameters in our model.The results we presented here were limited in the choice of the parameters so we could not discuss the behaviors in a wide range of the parameters.Other choices of parameters, such as f, c s and g k , would give richer and more interesting behaviors of polyelectrolytes in solvent mixtures.We will enumerate what we should do below and we hope to report some results about them in the near future.
(i) The composition of the mixture can be controlled easily.In this study, we consider only an asymmetric mixture ( f = 0.2), with which we expect the interaction caused by the inhomogeneous concentration field is strengthened as in colloidal systems.However, other conformations will be observed when the mixing fraction is changed.When approaching the critical point in a symmetric mixture ( f = 0.5), the interface tension is reduced 48 and the effective attractive interaction would also be weakened.Our preliminary simulations indicate that non-spherical domains containing semi-flexible chains, which are expanded by the electrostatic repulsion, are formed (not shown here).
In the present simulations, both of the relative fractions of the W-solvent rich phase and the polymer fractions are small even in the two phase region (w > w t ).Then, we observed in Fig. 2 that each polymer chain is compartmented in a droplet of the W-solvent rich phase.If we increase them, we will see that some polyelectrolytes are confined in a same droplet.With changing the interaction parameter w, aggregates of the polymer chains will also be observed.The aggregation process of the polyelectrolytes via the mechanisms discussed above is an interesting problem.
In Fig. 2 and 3, the polymer forms a compact conformation in the phase-separated state (w > w t ).This is because the diameters of the droplets are shorter than the polymer contour length.Since the periodic boundary conditions are employed, the maximum droplet size is uniquely determined by the relative fraction of the phase-separated domain and the simulation box size.If we consider an isolated chain in a larger cell or an infinite one, the droplet of the minority phase will grow larger than the polymer size and the polymer will adopt a more relaxed conformation in the droplet.In this sense, the gyration radius of the polyelectrolyte will depend on the system size.Also, we fix the total amount of salt as constant in the simulations.If we are interested in the infinite cell, we had better control its chemical potential of the ions, not their total amounts.This difference may lead to some quantitatively different results.
(ii) The salt concentration is also an important parameter, which is easy to change.Although only two salt concentrations are employed and the difference between them is not large, a quantitatively large difference in the conformational changes is observed (Fig. 4).If we change the salt concentration in a wider range, more drastic conformational changes will be observed.In particular, it is expected that the polyelectrolyte chain will behave as a neutral polymer when a sufficiently large amount of salt is added.However, our present simulation scheme does not apply to such a high salinity limit as noted in Appendix A. It is desired to improve the simulation method.
Also, we treat the ions by means of a coarse-grained model, namely, the Poisson-Boltzmann equation.Then, we do not consider the multivalent ions [15][16][17] and the correlation effect of the fluctuating ions. 2,51We have to treat the ions by a particle-based model, in order to see these effects.
(iii) Although the dependence of the conformation on the chain length is quite important, 4 our polymer chains consist only of N = 40 beads.Our simulations with shorter chains (N = 10 and 20) (not presented here) give essentially the same behaviors that we reported in this article.But, some qualitatively different behaviors may be observed when longer chains are used.For example, when a long charged polymer collapses into a compact conformation, the electrostatic repulsions among the monomers become enlarged.When the droplet of a finite interface tension no longer confines the polymer, Coulomb explosion locally occurs and some complex structure may be formed as a necklace. 4In the solvent mixtures, the ions are strongly confined in the W-solvent rich domain because of the solvation energy.Then, the entropic gain of the counterions is possibly different from that in a one-component solvent.To consider it, another scaling theory should be developed.
Also from the viewpoint of phase separation, the usage of long polymer chains will be interesting.As shown in Fig. 2 and  3, the polyelectrolyte chains behave as nuclei for the phase separation near w = w t .It is known that the interface tension suppresses the formation of the small droplets of the minority phase.So, when long polyelectrolytes are dissolved, they may induce phase separation more efficiently.Although the usage of long polymers is important and interesting, it is not easy to simulate them.When we use longer chains, we have to carry out simulations with much longer annealing times to equilibrate the system.Simulation improvements are desired.
(vi) In this work, we do not change the material parameters, such as g k , D 0 and D 1 .In particular, the simulated solvation This journal is © The Royal Society of Chemistry 2016 energy g k is rather small, in comparison with actual ions and solvents.Not only does the addition of a salt with large g k change the phase diagram largely, but it will also enhance the interaction between the charged polymer and the solvent. 44lthough our present simulation cannot adopt large g k , we have to consider its influence seriously.Also, D 1 plays an important role in our mechanism.It is also needed to consider the dependence of the conformation on it.

A Simulation scheme
Our simulation model deals with lattice and off-lattice spaces.The solvent composition f, the ion concentration c k , the local ionization degree a, the electrostatic potential F and the hydrodynamic flow u are solved in the lattice space.The grid size of the lattice space is set to d as shown in Fig. 1(a), and the periodic boundary conditions are employed.Here we describe our simulation scheme.
(i) First, we determined c k , F and a recursively, since they are coupled in a complicated manner.At this stage, we fixed y(t) and f(t).The electrostatic potential F (p) was obtained from eqn (10) with c k = c (p) k and a = a (p) by means of the Crank-Nikolson method.Here, the suffix p(=0, 1, 2,. ..) indicates the iteration index.(As the initial conditions for the iterations at t = 0, we set c (0) k = c s and a (0) = {1 + e D0ÀD1 f } À1 .Otherwise, we set c l ðpÞ a ¼ ln from eqn (6).Then, we had a (p+1) and c (p+1) k .Solving eqn (10) with them, we obtained the electrostatic potential at the next iteration F (p+1) .We iterated the above calculations until l (p) k converged to a certain value, and then we obtained c k (t + Dt) and a(t + Dt).
We should note that the range of the salt concentration is limited in our scheme.To solve the electrostatic potential in the lattice space, the Debye screening length should be longer than the simulation grid length d.Since the ions prefer to be dissolved in the W-solvent, their concentrations in the droplets of the W-solvent rich phase becomes much larger than their spatial averages.Then, we have to keep the salt concentration so small that the Debye screening length remains longer than the grid size even in the W-solvent rich phase.On the other hand, the concentration of the counterions in the simulation box is about 2 Â 10 À5 when the polyelectrolyte chain is highly charged, so that the solution of c s = 1 Â 10 À5 corresponds to a low salinity solution.
(ii) We calculated the hydrodynamic flow u by means of the Maker and cell method. 56The flow vector field was defined in a staggered lattice.Since the Reynolds number is considered negligibly small, we iterated to develop eqn (16) without updating f and y until Ð drjr@u=@tj Ð drj À frdF=df þ P i h i F i V i j o 10 À3 .The fluid particle dynamics method was originally developed to consider colloidal suspensions in a simple liquid.In the limits of DZ/Z 0 -N and d/a -0, this method is able to describe solid spherical particles with hydrodynamic interactions with high accuracy.Since we do not need to consider the boundary conditions at the colloidal surfaces, we can calculate the hydrodynamic flow efficiently, and then apply this method to some systems.The details of this numerical scheme are given elsewhere. 38In this study, the beads do not represent colloid particles, so we use a moderate value of DZ/Z 0 .
(iii) Now we knew u, c k and a.We updated the bead positions R i and the concentration field by means of the explicit Euler scheme with the same time increment Dt = 10 À2 d 2 /D f .The new polymer field y is given by eqn (1).Then we returned to (i) using y(t + Dt) and f(t + Dt), instead of y(t) and f(t).

Fig. 1
Fig.1 (a)A schematic picture of our bead-spring model.The circles with solid lines and the broken lines indicate the beads corresponding to the Kuhn monomers and the repeat units, respectively.a and b stand for the bead radius and the Kuhn length.The center of the mass of the i-th bead is given by R i and is defined in the off-lattice space.On the other hand, the solvent composition f, the ion concentration c k , the electrostatic potential F and the hydrodynamic flow u are calculated in the lattice space.(b) The phase diagram of the binary mixture under consideration.f is the solvent composition and w is the interaction parameter.The critical point is at f = 0.5 and w = 2.0.The solid and broken lines represent the coexistence and spinodal curves, respectively.In the simulations, we fix the averaged composition to f = 0.2, and two paths of w (cooling and heating) are considered (see below).With the employed parameters (c s r 4 Â 10 À5 and g k = 2.0), the influence of the salts on the phase diagram is negligibly small.
The permittivity parameters are set to e f = 1.4 e and e y = e.The viscosity parameters are Z 0 = 0.1k B T/(D f d) and DZ = 10Z 0 .Our results presented below are essentially insensitive to our choices of Z parameters.

3. 1
Fig.2(a) and (b) show the snapshots of the polyelectrolyte chains in the solvent mixtures at w = 1.8 and 2.35, respectively (ESI †).We solve N p = 10 chains, each of which consists of N = 40 beads, in a cubic box of V t = (256d).3The volume fraction of the polymer beads is 7.9 Â 10 À5 .We set the salt concentration to c s = 4 Â 10 À5 .The colors of the beads represent the degree of ionization.The interfaces between the phase-separated domains are also drawn with isosurfaces of f = 0.5.In the mixed state, far from the coexistence curve at w = 1.8 as in Fig.1(b), the polymers behave as flexible chains.Here sharp interfaces are not formed.At w = 2.35, on the other hand, they are compartmented in the droplets of the phase-separated domains.Here, the relative fraction of the W-solvent rich phase is 0.022 and this minority phase forms the droplets.Since the droplet size is smaller than the contour length of the polymers, the polymer chains are strongly distorted by the interface tension in the droplets.Since the polymer volume fraction is small and the relative fraction of the W-solvent rich phase is also small, each droplet

Fig. 2
Fig. 2 The snapshots of the polymer conformations in the solvent mixtures.The interaction parameters are set to w = 1.8 in (a) and 2.35 in (b).N p = 10 polymer chains, each of which consists of N = 40 beads, are dispersed in the cubic box with periodic boundary conditions.The salt concentration is c s = 4 Â 10 À5 .The interfaces of the phase-separated domains are shown with isosurfaces of f = 0.5.The colors of the beads represent the ionization.

Fig. 3
Fig. 3 The snapshots of the polymer conformation in the solvent mixture from w = 2.1 to w = 2.35.The spheres represent the beads and their colors represent the degree of ionization.The interfaces between the phase-separated domains (grey) are also depicted with isosurfaces of f = 0.5.Assuming an LCST mixture, the conformational changes in the (a) heating and (b) cooling processes are shown.The fraction of the W-solvent is f = 0.2 and the salt concentration is c s = 1 Â 10 À5 , which corresponds to ka = 0.055.

Fig. 4
Fig. 4 The changes in the gyration radii R g during heating (red circle) and cooling (blue square) are plotted with respect to w.The closed and open symbols are the gyration radii for c s = 1 Â 10 À5 and 4 Â 10 À5 , respectively.The broken line is the gyration radius of a neutral polymer of the same N.The mixture of f = 0.2 in the bulk phase separates in the hatched region (w > w t ).

Fig. 5
Fig. 5 Plots of the (a) the ionization degree a and (b) the adsorption amount G as functions of w.The fraction of the W-solvent is f = 0.2 and the salt concentration is c s = 1 Â 10 À5 .The red circles and blue squares correspond to a and G in the heating and cooling processes, respectively.

Fig. 6
Fig. 6 The sum of the electrostatic repulsive force and the attractive force caused by the inhomogeneous concentration field for the bead pair (N = 2) is plotted with the particle separation R 1,2 .It is estimated by the spring force F 1,2 = K(R 1,2 À b) by changing b.The positive and negative values indicate the repulsive and attractive interaction work between the two beads, respectively.The interaction parameter is changed from w = 2.0 to 2.2.