Name Regimes of electrostatic collapse of a highly charged polyelectrolyte in a poor solvent . †

We perform extensive molecular dynamics simulations of a highly charged, collapsed, flexible polyelectrolyte chain in a poor solvent for the case when the electrostatic interactions, characterized by the reduced Bjerrum length `B, are strong. We find the existence of several sub-regimes in the dependence of the gyration radius of the chain Rg on `B characterized by Rg ∼ `−γ B . In contrast to a good solvent, the exponent γ for a poor solvent crucially depends on the size and valency of counterions. To explain the different sub-regimes, we generalize the existing counterion fluctuation theory by including a more complete account of all possible volume interactions in the free energy of the polyelectrolyte chain. We also show that the presence of the condensed counterions modifies the effective attraction among the chain monomers and modulates the sign of the second virial coefficient in poor solvent conditions.


Introduction.
Charged polymers in solution, or polyelectrolytes (PEs), are ubiquitous in nature and play a significant role in our everyday life.Common examples of PEs include biologically important molecules like DNA, RNA and proteins [1][2][3] .PEs also find applications in industries such as chemical [4][5][6][7] , pharmaceutical [8][9][10][11] , food 12,13 etc.The mechanical and chemical properties of a PE depend on its conformational state, which could vary from being linear and extended to compact and collapsed.The conformational state is determined essentially by three properties: strength of electrostatic interactions in the system, entropy of the counterions and quality of the solvent.Determining the precise role of these properties in the conformational state of a simple model PE is fundamental to understand the physics of more realistic PE systems.
The strength of the electrostatic interactions in the system depends on the charge density of the PE chain and is quantified by the reduced Bjerrum length where ε is the dielectric permittivity of the solvent, k B is the Boltzmann constant, T is temperature, β = (k B T ) −1 , and a is the distance between neighbouring charged monomers of the PE chain.B quantifies the ratio of the electrostatic interaction energy between the neighbouring charged monomers and thermal energy: larger the value of B , stronger are the electrostatic interactions in the system.
For small B , a PE behaves like a neutral polymer and the counterions are dispersed away from the PE to occupy all accessible volume, resulting in a state with high entropy.With increasing B , the electrostatic interaction energy becomes comparable to the thermal energy and counterions begin to condense onto the PE, renormalizing its charge density 14 .The condensed counterions being in a close vicinity of the PE implies a lower entropy of the system.On the other hand, if B is not large the conformation of the PE is dictated primarily by the solvent quality.
The solvent quality, in turn, is determined by the relative strength of the attractive interactions between chain monomers and between monomers and solvent particles.In a good or thetasolvent, these attractive interactions are stronger for monomersolvent particle pairs, while in a poor solvent monomer-monomer attractive interactions dominate 15 .
Several experiments and simulations have shown that at large enough B , a like-charged PE chain undergoes a transition from extended to collapsed conformation regardless of the solvent quality [16][17][18][19][20][21][22][23][24][25] .This counterintuitive transition is driven by the condensation of counterions onto the chain, reducing the effective charge density.The nature of the effective attractive interactions driving the transition is not well understood and there are competing theories explaining their origin (see below).For the collapsed state, these theories predict that the gyration radius, R g , of a PE has the scaling form R g ∼ N 1/3 −γ B , where N is the degree of polymerization of the PE, and the exponent γ can po-tentially depend on system parameters.In the literature, there are three theoretical approaches 19,[26][27][28][29][30][31] , based on different physical models, to account for the electrostatics-driven collapse of PEs.All the theories predict a single collapsed regime, but differ in their prediction of the exponent γ characterizing the dependence of R g on l B .In the first approach 28 the system of charged, collapsed PE and the corresponding counterions is modeled as an amorphous ionic solid.In this theory, the system achieves minimum free energy when the PE is in a collapsed state and the R g is independent of the charge density of the PE, thus leading to γ = 0.In the second theory, pairs of condensed counterions and PE monomers are treated as fluctuating dipoles, and the attractive dipole-dipole interaction energy is identified as the driving force behind the collapse of the PE in both good as well as the poor solvents 26,29-31 * .The dipole theory proposes a scaling of R g on l B as R g ∼ −2/3 B , leading to a value of γ = 2/3 for both good and poor solvents.The third theory, referred to as counterion fluctuation theory, proposes that attractive interactions driving the collapse transition is a result of the density fluctuations of the condensed counterions leading to a negative pressure.This theory predicts γ = 1/2 for good solvent conditions.
In a recent molecular dynamics (MD) study of a flexible PE in a good solvent, we showed that a collapsed PE conformation demonstrates at least two different sub-regimes, one with γ = 1/2, and the other at larger B with γ = 1/5.The exponent γ in both regimes was found to be independent of the valency of counterions and PE chain length 34 .All inter-particle interactions, other than electrostatic, were repulsive in these systems.The counterion-fluctuation theory, originally developed for good solvent with a single exponent γ = 1/2, has been generalized by us 34 through the inclusion of higher order virial coefficients to reproduce both the collapsed regimes seen in our MD simulations.
It is more challenging, however, to study a collapsed state of a flexible PE in a poor solvent [35][36][37][38][39][40][41] since, unlike in a good solvent, there exist additional attractive interactions between monomers which compete with the repulsive part of electrostatic interactions.The valency of counterions determines the number of counterions condensed inside the collapsed globule and hence affects the monomer-monomer distance and thereby modifies the attractive volume interactions.Hence, the valency of the counterions which played no role in determining the exponent γ for a good solvent becomes significant for a poor solvent and possibly influences the exponent γ.In the present study, we report MD simulations for the collapsed phase of a strongly charged flexible PE in a poor solvent and find several novel collapse sub-regimes.The conformational behavior of a PE in a poor solvent, observed in MD simulations, can be explained by extending the modified counterion-fluctuation theory 19,34 with a more complete account of the volume interactions between all species in the system.We show that the new theory can describe the MD results for both good and poor solvents.
The rest of the paper is organized as follows.In Sec. 2, we define the model and give details of the MD simulations.In Sec. 3, we present our theory of the electrostatic collapse in a poor solvent and compare the predictions with data from MD simulations.A summary and discussion of our results are given in Sec. 4.

Model and Methods
We model a flexible PE chain as N monomers of charge e connected by harmonic springs of energy where k is the spring constant, a is the equilibrium bond length, and r is the instantaneous distance between the bonded monomers.The PE chain and N c = N/Z neutralizing counterions with a valency Z are placed in a box of linear size L. Pairs of all non-bonded particles (counterions and monomers) separated by a distance r i j interact through the volume (or van der Waals) interactions.Here we model these interactions by the Lennard Jones (LJ) 6 − 12 potential with a cutoff of r c : The values of ε LJ and r c are varied to model solvents of different quality (see below).The electrostatic energy between charges q i and q j separated by r i j is given by U c (r i j ) = q i q j εr i j .
In the simulations, we use a = 1.12σ , k = 500.0kB T /σ 2 , L = 370σ , N=204.All energies are measured in units of k B T , and we maintain temperature at 1 through a Nosé-Hoover thermostat.All distances are measured in terms of σ which we set to 1.The long-range Coulomb interactions are evaluated using particleparticle/particle-mesh (PPPM) technique.
The equations of motion are integrated in time using molecular dynamics simulation package LAMMPS 42 with a time step of 0.001.All the systems are equilibrated for 5 × 10 6 timesteps and the data presented in this paper are averaged over 5 × 10 6 timesteps of production runs.
We model different poor solvent conditions by choosing the following combinations of the LJ energy parameters ε LJ and cutoff distance r c for the monomer-monomer interactions: (i) ε LJ = 1 and r c = 2.5 and (ii) ε LJ = 1.5 and r c = 2.5.For all other volume interactions among counterions and between monomers and counterions, the LJ interactions are repulsive, ε LJ = 1.0 and r c = 1.0.These conditions are chosen in such a way that when the charge on the monomers is zero (neutral polymer), a PE chain adopts a collapsed conformation, mimicking poor solvent conditions.We also performed additional simulations in which the counterion size is varied.We note that all simulations reported in this paper have been performed for values of B , where the equilibrium configuration of a PE is a collapsed state with R g ∼ N 1/3 .

Results
To develop a generalized theory for electrostatic-driven collapse of a PE in a poor solvent, we start with the modified counterionfluctuation theory for a good solvent and retain the electrostatic term.As mentioned in the Introduction, the modified counterion fluctuation theory is able to explain the collapsed sub-regimes (with correct exponent γ) that were observed in MD simulations of a flexible PE in a good solvent 34 .However, in the case of a poor solvent we anticipate that the attractive volume interactions between chain monomers would cause even stronger collapse compared to that in a good solvent.In the following subsections, we generalize the theory of Ref. 34 develop a more complete account of the volume interactions and show that the counterionfluctuation theory is applicable to both good and poor solvent conditions.

Free energy of a PE system
We write an expression for free energy of the system F(R g ), taking into account the contributions from various entropic and energetic terms, and minimize it with respect to R g to obtain the dependence of R g on B .The different contributions to the free energy may be written as: where F id.ch , F en , F el , and F vol account for the free energy of an ideal chain, entropy of the counterions, the electrostatic interactions between all the charged particles, and the volume interactions between all the species respectively.F id.ch (R g ), the part of the free energy corresponding to the ideal chain reads 15,19,43 , where α = R g /R g.id is the expansion factor, with R g.id being the gyration radius of the ideal chain, R 2 g.id = Na 2 /6, where a is the inter-monomer distance.
The second part of the system free energy, F en (R g ), accounting for the entropy of the uncondensed counterions is proportional to the logarithm of the volume available for the counterions.It may be shown that 19,34 , where ρ = ρ c.in /ρ g with ρ c.in being the number density of counterions within the globule of gyration volume V g = 4πR 3 g /3 and ρ g = N c /V g = N/(ZV g ) is the maximal counterion number density corresponding to complete condensation.The value of R 0 corresponds to the simulation box size L in our MD simulations † .
The third part of the system free energy, F el (R g ), accounting † If the solution has the total volume Ω and comprises N PE chains, the volume per one chain reads Ω/N .It defines the characteristic size R 0 of a "Wigner-Seitz" cell 44 , that contains one PE chain, with its counterions, that is, (4π/3)R 3 0 = Ω/N .We consider dilute solutions, when R 0 R g .
for the electrostatic interactions among all charges is given by the counterion fluctuation theory 19 as: (8) The first term in the right hand side of Eq. ( 8) gives the mean-field result for the electrostatic interactions in the system.The second term describes the contribution to the free energy from the correlated fluctuations of the charge density and is beyond the Poisson-Boltzmann approximation 19 .Both Eqs.(7) and ( 8) are valid for dilute salt-free PE solutions, R 0 R g and long chains, N 1 19 ; in this limit, the terms proportional to R g /R 0 are negligibly small.Also, note that to obtain Eqs. ( 7) and ( 8), a simplified approximation of two constant counterion densities has been used -one for the condensed counterions and the other for the free ones 19 .Within this approach one can analyze all conformation regimes of a PE chain, from an extended to collapsed one, and to predict that the electrostatic collapse takes place, for N → ∞, as a first-order phase transition 19 .For a collapsed state, addressed here, a more accurate description of the counterion density, that accounts for its non-uniformity on the mean-field level, is performed in Appendix B. Eventually, this leads to the same results as obtained with the use of the more simple model, Eqs. ( 7) and (8).
Finally, the free energy accounting for the volume interactions (LJ), F vol (R g ), may be written as where F vol m.m. , F vol c.m. , and F vol c.c. are the free energy terms for the volume interactions among monomers, counterions and monomers, and among counterions respectively.
When the packing fraction of monomers is small, the free energy is well-approximated by the truncated virial expansion 15 .Keeping up to the third virial term, we obtain for F vol m.m. : where B 2 and B 3 are the second and third virial coefficients for monomer-monomer interactions, ρ m = N/V g is the average density of monomers inside the gyration volume and b = 2πa 3 /9 √ 6 .For a poor solvent, in general, we expect B 2 to be negative, making it necessary to keep the third virial term, the first positive term in the virial expansion.The free energy of the volume interactions involving counterions F vol c.c. and F vol c.m. will have a similar form as in Eq. (10).Combining these expressions for F vol m.m. , F vol c.m. and F vol c.c. we obtain the free energy F vol for the collapsed PE to be: where B2 and B3 are the renormalized second and third virial coefficients respectively that account for all volume interactions.
These coefficients read: where C k is the k-th virial coefficient for the counterioncounterion interactions and D l,k−l are the virial coefficient for monomer-counterion volume interactions of l-th order with respect to the counterions and (k − l)-th order with respect to monomers (see the Appendix A for the detail).The value of the bare virial coefficients B k , C k and D l,k−l are determined by the relative strength of the LJ interactions and the thermal energy k B T .Due to the dominance of repulsive forces in the monomercounterion and counterion-counterion volume interactions considered here, all coefficients C k and D l,k−l are positive, C k > 0 and D l,k−l > 0 for k ≥ 2 and 1 ≤ l ≤ k.We also assume that the renormalized third virial coefficient B3 is positive as well.
The sign of the renormalized second virial coefficient B2 is difficult to determine apriori.Though the "bare" value B 2 is negative for poor solvents, the renormalized coefficient B2 also depends on the counterion valency and the virial coefficients C k and D l,k−l [see Eq. ( 12)].If these positive coefficients are large and the valency Z is small, B2 could potentially change sign and become positive.The values of the virial coefficients C k and D l,k−l are determined by the size of counterions: the larger the counterions, larger are the virial coefficients.Hence, it is expected that large counterions will result in positive renormalized coefficient, B2 > 0 (see the Appendix A for the detail).These predictions will be checked in our MD simulations discussed in Sec.3.3.
The virial expansion truncated at the third virial term [see Eq. ( 11)] is valid for a poor solvent for low packing fractions.For a good solvent, we showed earlier that it is necessary to include the third virial term to explain the observed scaling of R g with B in MD simulations.At the same time, we showed that the packing fractions remain relatively low for a wide range of B values and inclusion of only the third virial term is sufficient 34 .However, for a poor solvent, due to the attractive volume interactions, we anticipate that the packing fractions will be larger than that of a good solvent at the same B .When the packing fraction is high, the truncated expansion Eq. ( 10) loses its accuracy.In this case, either all the terms in the virial expansion must be retained or an effective equation of state (EOS) for dense fluids like Flory-Huggins or van der Waals EOS or EOS for Lennard-Jones mixtures, e.g. 45 is to be used.In the EOS approach, one has to use model parameters, specific for a particular system, which may be obtained only by fitting to available experimental or numerical data.However, we choose the first approach, developed for dense gases 46 , where the virial coefficients are explicitly expressed in terms of the interaction potential.Including the additional terms in the virial expansion we obtain the general form, where the virial coefficients are, where C k l = k!/l!(k − l)! are the combinatorial coefficients (see the Appendix A for details).Note that any non-singular EOS may be expressed in the form of Eq. ( 14), where the virial coefficient depends on the particular EOS and its (model) parameters.
Combining the different contributions [Eqs.( 6), ( 7), ( 8) and ( 14)] the free energy in Eq. ( 5) attains the form, . In the collapsed regime, it is well known that, regardless of the solvent quality, almost all counterions are located within the gyration volume of the PE chain.Hence, the average counterion density inside the gyration volume will be ρ c.in ρ g = N c /V g resulting in ρ ≈ 1 in Eqs. ( 7) and ( 8).This allowing us to ignore the second and third terms in the right hand side of Eq. ( 16) ‡ .These two terms correspond to entropy of uncondensed counterions § and to the meanfield contribution to the electrostatic energy.Moreover, since R g ∼ N 1/3 in a collapsed regime, the expansion factor α ∼ N −1/6 1. Thus, the term proportional to α 2 in F id.ch (α) can also be neglected.It is straightforward to see that, in this limit, the term proportional to α −2 in F id.ch (α) is small compared to the volume terms and may be dropped as well.With these approximations, Eq. ( 16) for the free energy of a PE in a collapsed state, in terms of R g , reduces to the following expression, regardless of the solvent quality:

Dependence of R g and energies of a PE chain on B
The dependence of R g on B (R g ∼ −γ B ) may be obtained by minimizing the free energy, Eq. ( 17), with respect R g : The relative importance of different terms in Eq. ( 18) depends on the packing fraction Na 3 /R 3 g and the virial coefficients Bk .We consider three cases below: 1.When the packing fraction is small and the sign of B2 is posi- ‡ One can also use a more rigorous approach as in the original counterion fluctuation theory 19 and minimize the total free energy with respect to ρ.This will give the value of ρ very close to unity.Moreover, for N → ∞ the collapse occurs as a firstorder phase transition and ρ = 1 exactly in this limit.§ The entropy of the condensed counterions, is given in the Appendix B; it is much smaller for the addressed systems than other terms and may be neglected.
tive, the right hand side of Eq. ( 18), which arises from virial expansion, can be truncated at B2 , yielding the following dependence of R g on B , which implies γ = 1/2.

2.
When the packing fraction is small and the sign of B2 is negative, the most dominant positive term in the right hand side of Eq. ( 18) is the term proportional to B3 .Truncating the series at B3 , we obtain or γ = 1/5.Here we also assume that the electrostatic term strongly dominates over the term with the negative virial coefficient B2 .
3. When the packing fraction is large, truncation of virial expansion fails as discussed in Sec.3.1.In this case, all the terms in right hand side of Eq. ( 18) have to be retained.Qualitatively, we expect R g to decrease with B with a continuously varying "local " exponent γ.Moreover, we expect that γ < 1/5 and that γ will decrease with B , so that the decay of R g with B will be slower than . Indeed, as it follows from the above analysis, the value of γ is determined by the dominating virial coefficient.If for some interval of B the k-th virial coefficient is dominating, similar considerations lead to the conclusion, that R g ∼ 1/(3k−4) B for this interval.Since k > 3, one finds that γ < 1/5.With increasing B the number k that refers to the dominating coefficient is expected to increase.
We also derive the dependence of the internal energies, associated with the electrostatic and volume LJ interactions, on the gyration radius R g .Using the thermodynamic relation for the internal energy ) where E el and E LJ are the electrostatic and LJ components of the internal energy and the constants B k may be expressed in terms of the derivatives of the renormalized virial coefficients Bk with respect to temperature.
We showed earlier using MD simulations 34 19) and (20).We now check, using MD simulations, whether the presence of such sub-regimes is affected when the quality of solvent is changed to poor.The variation of R g with B and that of the internal energies with R g can also be measured in MD simulations.

MD results
We simulate a flexible PE in poor solvents of varying strengths, to obtain the dependence of R g on B as well as strength of the solvent.In B , with γ = 1/2 and 1/5 exist.For trivalent counterions, we do not observe the regime with γ = 1/2 for ε LJ = 1.0, while for the case ε LJ = 1.5 [see Fig. 1(f)] , both the regimes with γ = 1/2 and 1/5 are absent.Also, in the case of divalent and trivalent counterions, for larger B , an additional sub-regime appears where R g decreases as a convex function of B with a continuously decreasing slope (in the log-log plot).Qualitatively, this corresponds to an effective local power-law exponent γ, smaller than 1/5.Following the dis-cussion in Sec.3.2, the regimes with smaller values of γ that are observed in our simulations may potentially correspond to larger packing fraction of the collapsed PE, where the truncation of virial expansion at the third term fails.The MD data presented in Fig. 1 are consistent with the theoretical predictions of Sec.3.2.
The appearance of γ = 1/2 in Fig. 1, which corresponds to a positive second virial coefficient, is surprising since in poor solvent conditions, the bare second virial coefficient B 2 , when restricted to monomer-monomer interactions, is expected to be negative.Indeed when the charge of monomers is zero, the neutral PE always adopts a collapsed state, corresponding to the negative value of B 2 15 .Then from Eqs. ( 18) and ( 20), it can be seen that the largest possible value for the exponent γ should be 1/5, corresponding to the (positive) third virial term.Hence, we envisage that the presence of the counterions inside a collapsed globule leads to the change of the sign of B 2 for the range of LJ parameters and valency of counterions considered here.This agrees with the theoretical analysis in Sec.3.1 , where we argued that for large counterions with a low valency the effective second virial coefficient B2 could be positive, although the bare coefficient B 2 is negative, yielding the regime with γ = 1/2.Physically, the negative sign of B2 follows from the dominance of attractive volume interactions.Therefore, large counterions with a low valency (which implies the larger counterion density inside a globule) keep the chain monomers apart and reduce the effect of attractive volume interactions between them; this results in the alteration of the sign of B2 and the regime with γ = 1/2 is observed.At the same time, small counterions of high valency (which implies lower density of these inside the globule) can not effectively keep the monomers apart, so that their effective attractive volume interactions can yield a negative B2 .For these systems the regime with γ = 1/2 would be absent.In other words, if the regime with γ = 1/2 is observed for some system ( B2 > 0), the decrease of a counterion size should entail the alteration of the sign of B2 and hence disappearance of this regime, as for B2 < 0. On the other hand, if the regime with γ = 1/2 is absent ( B2 < 0), the increase of the counterion size should lead to the change of the sign of B2 and appearance of the regime with γ = 1/2, as for B2 > 0.
To confirm these predictions, we perform additional simulations varying the size of counterions, σ c−c : (1) for the case of monovalent counterions we decrease σ c−c and test whether the regime with γ = 1/2 disappears and (2) for the case of trivalent counterions, we increase σ c−c and check whether the regime with γ = 1/2 emerges.The MD data for these two simulations are shown in Fig. 2. As predicted, in the case of monovalent counterions, as the size of the counterions is reduced, we find that the regime with γ = 1/2 vanishes [see Fig. 2(a)], that is, B2 becomes negative.On the other hand, in the case of trivalent counterions the regime with γ = 1/2, absent for σ c−c = 1 [see Fig. 1(e)], appears when the counterion size increases up to σ c−c = 2 [see Fig. 2(b)].These data confirm that the presence of counterions inside the condensed phase can strongly modulate the effective attractive interactions between monomers of PE and can change the sign of the second virial coefficient B2 .This is a surprising result when compared to the neutral polymer case, in which the second virial coefficient is always negative in poor solvent.
Fig. 3 The dependence of (a) the electrostatic energy, E el , and (b) the monomer-counterion component of the LJ energy, E mc LJ , on the gyration radius R g for different counterion valencies.The data are for ε LJ = 1.0 and r C = 2.5 for monomer-monomer interactions and chain length N = 204.
The counterion fluctuation theory, as developed in Sec 3.1, may be further substantiated by computing the energies E el and E LJ from MD simulations.As can be seen from Eq. ( 21), the scaling of E el / B is independent of B across a wide range of B values and solvent quality and scales as inverse of R g .However, the dependence of E LJ on R g is more complicated and exhibits crossover from one regime to another similar to the dependence of R g on B .The results for dependence of E el on R g from our MD simulations, shown in Fig. 3(a), capture the linear dependence of E el on R −1 g for all valencies.The total LJ energy E LJ has contribution from attractive monomer-monomer interactions and repulsive interactions between other pairs.This makes it difficult to analyse E LJ in terms of second and third virial coefficients.Instead, we measure E mc LJ , the monomer-counterion component of the LJ energy as only repulsive interactions contribute to it.E mc LJ , shown in Fig. 3(b), captures the dependence of E LJ on R 3 g and R 6 g well, validating the free energy expression, Eq (17).
For E LJ , we also note that as the valency of the counterions is increased, powers of R g , corresponding to larger virial terms, appear.A possible physical explanation for this is that as the electrostatic interactions in the system become stronger and the packing fraction of all species inside the globule increases, more terms to account for the volume interactions in the collapsed state are needed.

Discussion and conclusions
In this paper, we studied theoretically and by means of MD simulations the nature of a collapsed state of a polyelectrolyte (PE) in a poor solvent for different strengths of electrostatic and volume interactions.Our main results can be summarized as follows.
We find that the collapsed phase of a PE chain consists of several sub-regimes, characterized by different power law scaling of R g with B (R g ∼ −γ B ).The most dominant sub-regimes have origins in the second (γ = 1/2) and third (γ = 1/5) virial coeffecients, when the packing fraction of PE monomers in collapsed phase is small.As the packing fraction increases, additional sub-regimes with an effective power law exponent γ < 1/5 start to appear.We observe that the effective γ decreases continuously with B , which demonstrates the contribution of additional terms in the virial expansion.
The present findings for collapsed PE in a poor solvent are different from our earlier results for a PE in a good solvent 34 .In that work, we found only two sub-regimes in the collapsed phase for PE over a similar range of B values studied here.It is to be noted that unlike the poor solvent case, the volume interactions between the monomers of a PE in good solvent were all repulsive in nature.This may result in a lower packing fraction of PE monomers throughout the range of B values considered, justifying the truncated virial expansion formalism and results in appearance of only two sub-regimes corresponding to second In the current study, through extensive MD simulations, we demonstrate that the condensation of counterions on a PE chain leads to an effective renormalization of the volume virial coefficients.The renormalized virial coefficients strongly depend on the valency of counterions and the strength of the poor solvent, which is characterized by the value of ε LJ -the energy parameter of the LJ potential.Surprisingly, the MD results for a particular set of parameters for the volume interaction potential show, via the presence of sub-regime with the exponent γ = 1/2, that the renormalized second virial coefficient B2 is positive.This is contrary to the expectation for B2 to be negative, as the simulations were performed for a poor solvent.We note that, while in MD simulations the solvent quality can be controlled by the interaction potential between monomers, in a theory this property is characterized by the sign of the second virial coefficient: B 2 is positive (B 2 > 0) for a good solvent and negative (B 2 < 0) for a poor one.It is not clear, however, whether this definition of solvent quality, based only on monomer-monomer interactions, remains meaningful for charged polymers in the presence of counterions.Within our generalized theory, we expect that the sign of the renormalized B2 will be manifested in MD simulations through the presence of the sub-regime, R g ∼ −γ B , with the exponent γ = 1/2 for B2 > 0 and absence of this regime for B2 < 0.
To understand the role of condensed counterions on the sign of renormalized B2 , we performed a theoretical analysis as well as additional simulations in which we varied the size of counterions and demonstrated that the appearance and disappearance of subregime with γ = 1/2 crucially depends on the size of the counterions [see Fig. 2].This dependence of the sign of the renormalized B2 on counterion size and valency occurs only for a poor solvent, as the condensed counterions can modulate the effective attractive interactions between monomers resulting in the alteration of the sign of B2 .For a good solvent, with repulsive interactions between the monomers, the counterion size and valency play no role.In Fig. 4, we present the results from MD simulations for three different counterion sizes (σ c−c ) which clearly show the presence of the sub-regime with γ = 1/2 for all the σ c−c .The results in Fig. 4 combined with those in Fig. 2, show that while the sign of B2 is unambiguous in a good solvent, the same is not true in the case of a poor solvent and depends on several system parameters such as strength of the solvent, valency and size of the counterions.This suggests that for charged polymers with attractive monomer-monomer volume interactions and in the presence of counterions, the sign of the second virial coefficient cannot be assumed to be always negative.This is in a striking contrast with collapsed neutral polymers, which have the same attractive monomer-monomer interactions, where the sign of B 2 is always negative.
In the current study, we also develop a generalized theory for a collapsed regime of a PE in good and poor solvents based on counterion-fluctuation theory 19 , by explicitly considering the interactions between monomer-monomers, monomers-counterions and counterions-counterions.The good solvent case considered earlier 34 is a special case of this generalized counterionfluctuation theory.We further validate the predictions of the generalized counterion-fluctuation theory through the comparison of the theoretical and MD simulation results for dependence of different parts of the internal energy of the system with R g .The original counterion-fluctuation theory 19 predicts that the electrostatic internal energy of the system, E el , scales with the gyration radius as E el ∼ R −1 g .At the same time, the dependence of internal energy associated with the volume (LJ) interactions is expected to be different for different sub-regimes, similar to the dependence of R g on B described above.The MD simulation results are in complete agreement with these predictions [see Fig. 3] when the packing fraction is small.We also note that the values of R g at which the crossover from one sub-regime of E LJ (R g ) to another take place [see Fig. 3] coincide with the values of R g where the crossover between regimes with different exponents γ are detected [see Fig. 1].
The different regimes of electrostatic collapse observed in the MD simulations reported here may be experimentally realized as well.In organic solvents, like mixtures of water and alcohol, the dielectric permittivity may be reduced (≤ 18) leading to a large Bjerrum length ( B ∼ 10) and collapse of a biological PE like DNA has been observed experimentally 20 under such conditions.We note that if all monomers of DNA were charged, one obtains B 29 and the range of these experimentally realizable B values are similar to the ones considered in this study.Another system with a comparably high linear charge densities is solution of charged worm-like micelles [47][48][49] .We hope that experimentalists find it interesting to check the existence of different collapse regimes for different counterion valencies and sizes.
Based on our findings, we conclude that the effective attractive electrostatic interactions in systems of like-charged polymers in the presence of counterions is described well in terms of correlated fluctuations of counterions, as has been proposed in the counterion fluctuation theory 19 .The electrostatic term of the free energy of the system, based on the counterion fluctuation, is independent of the solvent quality.Also, unlike the case of neutral polymers, the sign of the second virial coefficient of charged polymers such as PEs studied here cannot be determined apriori to be negative, and that the sign of the second virial coefficient can be changed by changing the size and valency of the counterions.We note that none of the other existing theories of effective electrostatic interactions of a PE 26,[28][29][30][31] can explain the sequence of collapsed sub-regimes or the scaling of the electrostatic energy with B , as seen in our MD simulations.

Acknowledgments.
The simulations were carried out on the supercomputing machines Annapurna, Nandadevi and Satpura at the Institute of Mathematical Sciences.

Appendix A
The free energy for the volume interactions among monomers may be written in the form of virial expansion as, where B k is the k-th virial coefficient for monomer-monomer interactions and ρ m = N/V g is the average density of monomers inside the gyration volume.Similar to Eq. ( 23), the free energy of the volume interactions of counterions reads, where ρ c.in N c /V g = N/ZV g is the average counterion density inside the gyration volume and we approximate it by the corresponding density, when all counterions are condensed.C k are the virial coefficients for the counterion-counterion interactions.
In the above equation we neglect the term that accounts for the volume interactions of the free counterion F f.c.vol c.c. , since it is negligibly small as compared to F vol c.c. : Furthermore, the monomer-counterion volume interactions are described by the term, where C k l = k!/l!(k − l)! are the combinatorial coefficients and D l,k−l is the k-th virial coefficient for monomer-counterion volume interactions which refers to l counterions and k − l monomers.
Using Eqs. ( 23), ( 24) and ( 26) one can write the part of the free energy responsible for the volume interactions in the system in the following compact form, where the renormalized virial coefficients Bk , that account for all volume interactions, are defined as Here we consider a collapsed state of a PE chain and the main difference of a globular state, as compared to that of a coiled state is that a "globule contains a number of uncorrelated parts of the chain" so that a concept of quasimonomers 50 and the corresponding virial expansion for pressure or free energy is valid 43 .In our study we use the LJ potential Eq. ( 3) to model all volume interactions.For the monomer-monomer interactions we use the cutoff distance r c = 2.5σ , that is, these volume interactions include both attractive and repulsive forces.For the monomer-counterion and counterion-counterion interactions we use the cutoff r c = σ , which corresponds to purely repulsive forces.The virial coefficient B 2 reads ¶ (29) ¶ For simplicity we ignore the contribution to the second virial coefficient B 2 from the third virial coefficient B 3 that appears due to the chain connectivity 43 .It may be shown that for the addressed range of parameters this contribution is positive and not large to yield a qualitative difference.Similar expressions apply for the coefficients C 2 and D 1,1 , with corresponding change of the LJ potential.Using these expressions for the virial coefficients B 2 , C 2 and D 1,1 one can compute, using Eq. ( 12) the renormalized second virial coefficient B2 .The results are presented in Fig. 5, where we demonstrate the dependence of sign of B2 on the size of the counterions σ c and their valency Z.As it may be seen from the Fig. 5, the renormalized coefficient B2 changes its sign from negative to positive with increasing size of counterions.This effect corresponds to an effective change of the solvent quality due to counterion condensation.

Appendix B
Qualitative analysis.In the counterion fluctuation theory, 19 a simplified approximation of two constant counterion densities, one for the condensed counterions inside the gyration volume, r < R g (r is the distance from the PE center of mass 19 ), and the other one for free counterions, r > R g , has been used.That is, the following model for the reduced counterion density ρc (r) = ρ c (r)/(n/Z) has been exploited: where ρ is a constant, n = N/V g (V g = 4πR 3 g /3) is the concentration of monomers inside the globule, which is assumed to be constant, n/Z is the maximal concentration of counterions within the globule, θ (x) is the unit step-function and η = R g /R 0 1, with R 0 being the radius of the "Wigner-Seitz" cell, which quantifies the volume per one PE chain 44 .Using the above counterion density, Eq. ( 30), one obtains for N 1 and R 0 R g the according expressions (7) for the entropic part of the free energy of the system and ( 8) for the part attributed to the electrostatic interactions on the mean-field level (i.e.without the account of charge fluctuations).The approach based on the model distribution (30)  allows to analyze (at least on qualitative level) all conformational states of a PE chain, from a completely stretched to a collapsed one 19 .
Below we consider a more accurate description of the counterion density, which is not constant, but determined by the electric field.We will restrict ourself to the mean-field, that is, Poisson-Boltzmann (PB) analysis. Let be the charge density, where n m (r) = nθ (R g − r) is the monomer density and ρ c (r) = (n/Z) ρc (r) is the counterion density.Then the Poisson equation, ∇ 2 ψ(r) = −(4π/ε)q(r) for the reduced electrostatic potential u(r) = eψ(r)/k B T reads, which coincides with the according equation for a uniformly charged sphere, penetrable for the surrounding microions (coions and counteirons) 44,51,52 .Using the dimensional distances, x = r/R g we recast the above equation into the form: where we introduce κ2 = 3N(l B /R g ) which has the meaning of the dimensionless (in units of 1/R g ) inverse Debye length, built up on the monomer concentration.In the mean-field approach the Boltzmann distribution for the counterions is used ρc (r) = Ce Zu(r) .When κ is small, which corresponds to weak electrostatic interactions, the reduced potential u is also small; this allows the linearization ρc (r) ≈ C(1+Zu(r)), which further yields the linearized PB equation.In this case one can solve Eq. ( 33) analytically and compute the entropic and mean-field electrostatic parts of the counterion free energy.The resulting expressions are qualitatively similar to that obtained with the simplified counterion distribution (30) used in the main text.Small values of κ generally correspond to a PE chain being in an expanded (or coiled for B → 0) state, which is beyond the scope of the present study.
In contrast, for a collapsed state, R g ∼ N 1/3 , we observe that for N 1 the value of κ2 is very large, κ2 1; for the studied systems it varies in the range of κ2 ≈ 150 − 1700.From Eq. ( 33) then follows, Hence, for κ−2 1 the approximation ρ ≈ 1, used in the main text for the collapsed state becomes valid [see Eq. ( 30)].In other words, the above equation means that the deviation of the couinteiron distribution from the distribution ρ c = (n/Z)θ (R g − r), coinciding (up to a factor 1/Z) with the distribution of monomers is small and is not negligible only for a narrow shell of thickness δ around r = R g .To find δ we write u = χ((x − 1)/δ ), where χ = O(1), then ∇ 2 χ ∼ δ −2 χ.Substituting this into Eq.( 34), we conclude that the thickness δ scales as δ ∼ κ−1 1 and obtain for the counterion density: where χ (y) = d 2 χ/dy 2 is a very narrow function around r = R g with the width of R g / κ R g .For κ 1, the leading term for ρc (r) in Eq. ( 35) coincides with that of the simplified density distribution (30) with ρ = 1.Hence, keeping only the leading term, which dominates at κ 1, corresponds to the approximation ρ ≈ 1 used in the main text, that yields Eq. ( 17).From the above Eq.( 35) also follows that the contribution from the subleading term, ∼ χ , to the system free energy in a collapsed state will be small, decaying as 1/ κ for large κ 1.These conclusions agree with the results of the studies 44,51,52 , where very similar systems have been explored.
Quantitative estimates.We start from Eq. ( 33) and exploit the Boltzmann approximation for the counterion density, ρc (r) = Ce Zu(r) , which implies u = (1/Z) log( ρc /C).For simplicity and notation brevity we consider the case of Z = 1, then Eq. ( 33) may be written as we assume that ρc may be written as where δ ρ is small everywhere except the narrow region around x ≈ 1.We consider separately solutions of Eq. ( 36) for x < 1 and x > 1.
• For x < 1 we have ρc (x) = 1 + δ ρ 1 , so that Eq. ( 36) reads, where we neglect δ ρ 1 , as compared to 1 and the quadratic term (δ ρ 1 ) 2 as compared to other terms.The solution of the above equation, finite at x = 0, reads • For x > 1 we have ρc (x) = δ ρ 2 , and Eq. ( 36) may be written as, where we neglect the cubic term κ(δ ρ 2 ) 3 as compared to other terms.Solving the above equation we obtain, The integration constants A, B and b may be found from the conditions of continuity of ρc and of its derivatives at x = 1 and using the normalization condition, Hence the counterion density reads, where n = N/V g , and δ ρ 1/2 (x) are given by Eq. ( 39).Using the above counterion density one can compute the entropic part of the counterion free energy: where Λ c is the thermal wavelength of counterions, φ = na 3 = Na 3 /V g is the packing fraction of monomers inside the PE globule, η = R g /R 0 1 and we skip unimportant constant terms.In the last line of Eq. ( 41) the first term gives the entropy of condensed counterions inside globule and may be neglected, as compared to other terms, for B 1 and 1/α 1; the second term scales as 1/ κ and also may be neglected for κ 1, as it has been shown by the above qualitative analysis.To derive the above result we keep, for κ 1, only the leading terms in the exponential integral function expansion, Ei( κ) = e κ ( κ−1 + κ−2 + 2 κ−3 + . ..).
The mean-field electrostatic part of the counterion free energy reads: |r 1/2 |<R 0 dr 1 dr 2 q(r 1 )q(r 2 ) where q(r) is the charge density.Using Eq. ( 31) for q(r) and Eq.(40) for the counterion density one obtains from Eq. ( 42), where we again use the above expansion for Ei(x) and keep only the leading terms.Finally we find for κ 1: For B 1 and collapsed state, 1/α 1, these terms may be safely neglected, as compared to other terms of the system free energy, corresponding to the correlated fluctuations of the charge density [the last term in Eq. ( 8)] and the volume interactions term [Eq.( 14)].Obviously, the approximation F en + F el m.f.≈ 0, valid for B 1 and 1/α 1 is equivalent to the approximation ρ ≈ 1,

a
The Institute of Mathematical Sciences, C.I.T. Campus,Taramani, Chennai 600113, India b Homi Bhabha National Institute, Training School Complex, Anushakti Nagar, Mumbai-400094, India c Department of Mathematics, University of Leicester, Leicester LE1 7RH, United Kingdom that in the case of good solvent, the collapsed phase of PE exhibits two subregimes, characterized by different scaling relation of the gyration radius R g with B namely, R g ∼ −1/2 B and R g ∼ −1/5 B as predicted by Eqs. (

3 Fig. 1
Fig. 1 Variation of the gyration radius R g with B for a PE chain with counterion valency Z = 1 (a,d); Z = 2 (b,e) and Z = 3 (c,f).The data are for poor solvent conditions (ε LJ = 1.0 and ε LJ = 1.5 and r c = 2.5 for monomers).The solid straight lines correspond to power laws with exponents γ = 1/2 and γ = 1/5 as predicted by the theory.

Fig. 4
Fig.4 The dependence of the gyration radius R g for different radius of counterions for good solvent (ε LJ = 1 and r c = 1).The chain length is N = 204 and the valency of counterions Z = 1

Fig. 5
Fig. 5 The dependence of the renormalized second virial coefficient B2 on the ratio of the counterion and monomer sizes, σ c /σ m and the counterions valency Z, as it follows from Eq. (12).The bare coefficients B 2 , C 2 and D 1,1 are computed using Eq.(29) with the corresponding Lennard-Jones potentials: ε LJ = 1.5 and r c = 2.5 for the monomer-monomer interactions, ε LJ = 1 and r c = 1 for the counterion-counterion interactions and ε LJ = √ 1 • 1.5 and r c = (2.5+1)/2= 1.75 for the monomer-counteiron interactions.Depending on the valency of the counterions the renormalized coefficient B2 changes its sign with the increasing size of the counterions. 4π