Sequence disorder-induced first order phase transition in confined polyelectrolytes

V. Stepanyan Yerevan State University, Yerevan, Armenia A. Badasyan University of Nova Gorica, Nova Gorica, Slovenia V. Morozov Yerevan State University, Yerevan, Armenia Y. Mamasakhlisov [email protected] Yerevan State University, Yerevan, Armenia Institute of Applied Problems of Physics, Yerevan, Armenia R. Podgornik [email protected] School of Physical Sciences, University of Chinese Academy of Sciences, Beijing, 100049, China CAS Key Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing, 100190, China Kavli Institute for Theoretical Sciences, University of Chinese Academy of Sciences, Beijing, 100049, China Wenzhou Institute of the University of Chinese Academy of Sciences, Wenzhou, Zhejiang 325011, China Department of Physics, Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19, 1000 Ljubljana, Slovenia

Abstract

We consider a statistical mechanical model of a generic flexible polyelectrolyte, comprised of identically charged monomers with long range electrostatic interactions, and short-range interactions quantified by a disorder field along the polymer contour sequence, which is randomly quenched. The free energy and the monomer density profile of the system for no electrolyte screening are calculated in the case of a system composed of two infinite planar bounding surfaces with an intervening oppositely charged polyelectrolyte chain. We show that the effect of the contour sequence disorder, mediated by short-range interactions, leads to an enhanced localization of the polyelectrolyte chain and a first order phase transition at a critical value of the inter-surface spacing. This phase transition results in an abrupt change of the pressure from negative to positive values, effectively eliminating polyelectrolyte mediated bridging attraction.

I Introduction

Interactions involving biologically relevant heteropolymers, such as nucleic acids, polysaccharides and polypeptides Leckband and Israelachvili (2001), are generally of two types Rubinstein and Colby (2003). The relatively generic and longer-ranged type interactions originate in the electrostatic charges, being at the origin of the polyelectrolyte phenomenology Dobrynin and Rubinstein (2005), and is often invoked as the primary cause of their solution behavior Muthukumar (2023). The more specific, shorter ranged interactions French et al. (2010), related to non-universal chemical identity of the various monomer units Elias (2012), cannot be quantified by a single parameter, analogous to the electrostatic charge, and are standardly invoked when the observed polyelectrolyte phenomenology cannot be understood solely based on electrostatic interactions Smiatek (2020).

The origin of the specific short-range interactions can be understood within different theoretical frameworks as stemming from the differences in polarizability of the monomer subunits in the context of van der Waals interactions between polymers Lu et al. (2015), from the differences in the dissociation properties of ionizable moieties within a generalized charge-regulated electrostatics Borkovec et al. (2001), or within a detailed description of water structuring models accounting for the solvent mediated interactions Blossey and Podgornik (2022), to list just a few. While each of these various theoretical frameworks of accounting for specific interactions within the context of polymer solution theory leads to differing predictions on the details of these interactions, they all agree in the sense that these interactions are important to understand the behavior of polymers in aqueous solutions.

The diversity of these specific interactions, originating from different types of monomer units, such as 4 canonical nucleotides in nucleic acids Bloomfield et al. (2000) or 20 canonical amino acids in polypeptides Finkelstein and Ptitsyn (2016), is therefore much more pronounced in terms of short range non-electrostatic interactions than long range electrostatic effects French et al. (2010) and is closely related to the sequence of monomers along the chain, leading to interesting effects for polymer-polymer interaction as well as polymer-substrate interaction. Gauging the relative importance of the generic long-range Coulomb interaction vs. the sequence specific short-range interaction is particularly important when trying to assess their role in the complicated multi-level problem of RNA folding Fallmann et al. (2017), RNA substrate interaction Poblete et al. (2021), when formulating realistic models of RNA virus co-assembly Twarock et al. (2018), or manipulating electrostatic repulsion in the background of sequence-specific Watson-Crick base pair interaction of ssRNAs Adhikari et al. (2020). In all these cases an interplay of both types of interaction needs to be properly accounted for.

In what follows we will investigate a model hetero-polyelectrolyte chain, characterized by both generic long-range Coulomb interactions between its equally charged monomers, with superimposed specific short-range interactions depending on the type of the monomer along its sequence. We will assume that the monomer sequence can be modeled as a quenched, random sequence, represented by a single random variable characterized by a Gaussian distribution with zero mean. We will write the free energy of the system within the self-consistent field theory of a confined, single chain within a replica symmetric Ansatz, and investigate the interactions between two planar surfaces, charged oppositely from the chain monomers, as a function of the separation between them. We will show that the model confined polymer chain exhibits an additional localization induced by sequence disorder, akin to the Anderson localization, and will assess the relative importance of the disorder effects superimposed on an electrostatic background.

II Model and Methods

II.1 Model

We consider a generic flexible hetero-polyelectrolyte, e.g. ssRNA, see Fig. 1, confined between two planar surfaces, comprised of equally charged monomers interacting via long-range Coulomb interactions, as well as short range "chemical" interactions dependent on the monomer sequence. While the charge of the monomers is assumed to be the same for all the monomers, the short range interactions are sequence specific, with this specificity being encoded in the sequence disorder. This implies that instead of considering a detailed sequence of monomers characterized by different chemical identities, we characterize the strength of the short-range interactions by a random, Gaussian distributed, variable. Similar systems without the sequence disorder affecting the short-range interactions Borukhov et al. (1999); Shafir et al. (2003); Siber and Podgornik (2008); Podgornik (1991); Andelman (1995); Budkov and Kalikin (2023) and the problem of a polymer chain with excluded volume interactions in quenched random media Muthukumar (1989); Hribar-Lee et al. (2011); Bratko and Chakraborty (1995) have been studied before. Our main goal here is to address the interplay between the disordered short-range interactions and long-range Coulomb interactions. In order to simplify the calculations we consider a salt free system with the surfaces oppositely charged from the polymer chain.

Refer to caption — Figure 1: The negatively charged polyelectrolyte chain is confined between two positively charged infinite plates with surface charge density $\sigma$ . The charge per monomer is $pe$ , where $e$ is the unit charge, and the short-range interaction between two parts of the chain is $v_{\tau\tau^{\prime}}({\bf x})=v_{0}\xi_{\tau}\xi_{\tau^{\prime}}\delta({\bf x})$ , where $\xi_{\tau}$ is the random sequence disorder variable and $\tau$ is the position along the contour of the chain.

The charge per monomer is assumed to be $pe$ , where $e$ is the electron charge and $0<p\leq 1$ . The value $p=1$ corresponds to the fully charged polyelectrolyte, while $p<1$ describes the weak polyelectrolyte, where ion binding/condensation partially neutralises the charges. The position of a monomer is characterized by a continuous Cartesian coordinate ${\bf r}(\tau)$ , where $\tau\in[0,N]$ marks the position along the contour of the chain, and $N$ is the dimensionless contour length of the chain. The random variable characterizing the strength of the short-range interaction is $\xi_{\tau}$ . $\{\xi\}$ s are considered to be independent, quenched random variables with the probability distribution

{\cal P}\{\xi\}=\prod_{\tau}p(\xi_{\tau}),

(1)

where the sequence disorder is supposed to be Gaussian in order to make the formalism tractable, therefore

p(\xi_{\tau})=\frac{\mathit{e}^{-\frac{(\xi_{\tau}-\bar{\xi})^{2}}{2\xi^{2}}}}{\sqrt{2\pi\xi^{2}}}.

(2)

For the sake of simplicity we restrict ourselves to the vanishing average value, $\bar{\xi}=0$ , case.

The Hamiltonian of the system in this case has a general Edwardsian form with

	$\displaystyle\beta{\cal H}[{\bf r}(\tau)]=\frac{3}{2\ell^{2}}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf r}(\tau))^{2}+$
	$\displaystyle+\frac{\beta}{2}\int_{0}^{N}\,d\tau\int_{0}^{N}\,d\tau^{\prime}v_{\tau\tau^{\prime}}({\bf r}(\tau)-{\bf r}(\tau^{\prime}))+\beta V[\{{\bf r(\tau)}\}],$		(3)

where $\ell$ is the Kuhn length, $V[\{{\bf r}(\tau)\}]$ describes all the long-range interactions in the system, i.e., the Coulomb interactions between monomers as well as the Coulomb interactions between monomers and the surface, while

v_{\tau\tau^{\prime}}({\bf x})=v_{0}\xi_{\tau}\xi_{\tau^{\prime}}\delta({\bf x}),

(4)

with a positive interaction constant $v_{0}>0$ . Thus, the different monomers will attract, but the similar ones will repel.

II.2 Partition function and free energy

We formulate the partition function of this three-dimensional polymer chain by consistently including the contribution of the monomer sequence disorder and analyzing the combined effect of randomly quenched short-range and homogeneous long-range interactions on the confinement free energy. The partition function of the model polyelectrolyte chain for every fixed sequence $\{\xi\}$ is given by

Z\{\xi\}=\int\,{\cal D}[{\bf r}(\tau)]~{}e^{-\beta{\cal H}[{\bf r}(\tau)]}.

(5)

In the limit of a very long chain, $N\rightarrow\infty$ , consistent with the ground-state dominance Ansatz, the disorder part of the free energy obeys the principle of self-averaging yielding Binder and Young (1986); Parisi (2023):

{\cal F}=-k_{B}T\langle{\ln Z\{\xi\}}\rangle_{{\cal P}},

(6)

where $\langle{...}\rangle_{{\cal P}}$ stands for the average with the distribution function (Eq. 1). Self-averaging physically means that the distribution of the free energy has a vary narrow peak in the vicinity of its maximum, corresponding to the mean value of the free energy (Eq. 6).

The quenched free energy (Eq. 6) can now be estimated using the replica trick Binder and Young (1986); Parisi (2023)

-\beta{\cal F}=\lim_{n\rightarrow 0}\frac{\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}-1}{n},

(7)

where $\beta=(k_{B}T)^{-1}$ . In order to calculate the $n$ -replica partition function

	$\displaystyle\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}=$
	$\displaystyle\qquad=\int\,{\cal D}{\bf r}e^{-\frac{3}{2\ell^{2}}\sum_{a=1}^{n}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf r}^{a}(\tau))^{2}-\beta\sum_{a=1}^{n}V\{{\bf r}^{a}\}}\times$
	$\displaystyle\quad\int\,{\cal D}\xi{\cal P}\{\xi\}e^{-\frac{\beta v_{0}}{2}\int_{0}^{N}\,d\tau\int_{0}^{N}\,d\tau^{\prime}\xi_{\tau}\xi_{\tau^{\prime}}\sum_{a=1}^{n}\delta({\bf r}^{a}(\tau)-{\bf r}^{a}(\tau^{\prime}))},$		(8)

and thus to estimate the free energy, several order parameters need to be introduced, namely, the inter-replica overlap $q_{ab}({\bf x},{\bf x}^{\prime})$ , with $a<b$ , the density of $a$ -th replica monomers $\rho_{a}({\bf x})$ , and the density of the random variable $\xi_{\tau}$ of $a$ -th replica monomers $m_{a}({\bf x})$ as

$\displaystyle q_{ab}({\bf x},{\bf x}^{\prime})$	$\displaystyle=$	$\displaystyle\int_{0}^{N}\,d\tau\delta({\bf x}-{\bf r}^{a}(\tau))\delta({\bf x}^{\prime}-{\bf r}^{b}(\tau))$
$\displaystyle\rho_{a}({\bf x})$	$\displaystyle=$	$\displaystyle\int_{0}^{N}\,d\tau\delta({\bf x}-{\bf r}^{a}(\tau))$
$\displaystyle m_{a}({\bf x})$	$\displaystyle=$	$\displaystyle\int_{0}^{N}\,d\tau\xi_{\tau}\delta({\bf x}-{\bf r}^{a}(\tau)).$	(9)

Introducing next the Lagrange multipliers ${\hat{q}}_{ab}({\bf x},{\bf x}^{\prime})$ and ${\hat{\rho}}_{a}({\bf x})$ , conjugated to the first two order parameters in (Eq. II.2), allows us to rewrite the $n$ -replica partition function as Garel et al. (1997)

	$\displaystyle\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}$		$\displaystyle\propto e^{-\frac{nV}{2\tilde{v}}\ln(2\pi\beta v_{0})}\times$		(10)
			$\displaystyle\times\int\,{\cal D}{\rho}\,{\cal D}{{\hat{\rho}}}\,{\cal D}{q}\,{\cal D}{\hat{q}}~{}e^{g(\rho,{\hat{\rho}},q,{\hat{q}})+\ln\zeta({\hat{\rho}},{\hat{q}})},$		(10)

where $\tilde{v}$ is the monomer volume. The explicit expression for the functional $g(\rho,{\hat{\rho}},q,{\hat{q}})+\ln\zeta({\hat{\rho}},{\hat{q}})$ , as given by (Eq. A), allows us to obtain the final form of the free energy functional as :

	$\displaystyle-\beta{\mathcal{F}}(\rho,{\hat{\rho}},q,{\hat{q}})=\lim_{n\rightarrow 0}\frac{1}{n}\bigg{[}-\frac{nV}{2\tilde{v}}\ln(\beta v_{0})-$
	$\displaystyle\frac{1}{2}{\rm Tr}\ln\hat{A}-\frac{1}{2}{\rm Tr}\ln\bigl{\{}I+\xi^{2}\hat{A}^{-1}\hat{B}\bigr{\}}-$
	$\displaystyle\beta\sum_{a}W_{el}^{a}+i\sum_{a}\langle\rho_{a}\|\hat{\rho}_{a}\rangle+i\sum_{a<b}\langle q_{ab}\|\hat{q}_{ab}\rangle+\ln\zeta\bigg{]}$		(11)

where the operators $\hat{A}$ and $\hat{B}$ have been defined as $n\times n$ operator matrices $\langle{\bf x}|{\hat{A}}_{ab}|{\bf x}^{\prime}\rangle=\delta_{ab}\delta({\bf x}-{\bf x}^{\prime})[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})]$ and $\langle{\bf x}|{\hat{B}}_{ab}|{\bf x^{\prime}}\rangle=(1-\delta_{ab})q_{ab}({\bf x},{\bf x}^{\prime})$ . To estimate the free energy of the system in the self-consistent field approximation we need to minimize the functional (Eq. II.2) over the fields $\rho$ , ${\hat{\rho}}$ , $q$ , and ${\hat{q}}$ . For relevant details and derivations see the Appendix A.

The term $\ln\zeta({\hat{\rho}},{\hat{q}})$ describes the entropy of the polymeric chain in terms of a “quantum-like" particle confined in the restricted volume. The energetic spectrum of such a system is expected to have a gap in the energy spectrum between the ground state and the rest of the spectrum. Thus, the free energy of the system can be calculated using the ground state dominance Ansatz, where the partition function (Eq. 10) is dominated by the ground state "wave function" $\psi({\bf r})$ with energy $\mathcal{E}_{0}$ Garel et al. (1997) which is equivalent to the one-replica density field in the equation (Eq. II.3) (see below)

\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}\propto e^{-N\mathcal{E}_{0}}\int\,d^{3}R\psi({\bf 0})\psi({\bf R}).

(12)

The ground state dominance Ansatz is known to work in the case when the polymer length is long enough compared to the separation between the bounding surfaces Li et al. (2018).

II.3 Replica - symmetric solution

In order to estimate the free energy of the system we maximize the functional $g(\rho,{\hat{\rho}},q,{\hat{q}})+\ln\zeta({\hat{\rho}},{\hat{q}})$ from (Eq. 10). We will restrict ourselves to the replica - symmetric order parameters: $q_{ab}({\bf x},{\bf x}^{\prime})=q({\bf x},{\bf x}^{\prime})$ , $\hat{q}_{ab}({\bf x},{\bf x}^{\prime})=\hat{q}({\bf x},{\bf x}^{\prime})$ , $\rho_{a}({\bf x})=\rho({\bf x})$ and $\hat{\rho}_{a}({\bf x})=\hat{\rho}({\bf x})$ . See Appendix Eq. B for details.

The electrostatic part of the free energy within the mean-field Poisson-Boltzmann electrostatics can be obtained as Borukhov et al. (1999); Shafir et al. (2003); Siber and Podgornik (2008)

	$\displaystyle W_{el}(\rho,\varphi,c^{\pm})=\int\,d^{3}{\bf x}\bigg{\{}-\frac{\epsilon\epsilon_{0}}{2}(\nabla\varphi({\bf x}))^{2}+$
	$\displaystyle\varphi({\bf x})\bigg{[}ec^{+}({\bf x})-ec^{-}({\bf x})-pe\rho({\bf x})+\rho_{surf}({\bf x})\bigg{]}+$
	$\displaystyle\sum_{i=\pm}\bigg{[}k_{B}T(c^{i}({\bf x})\ln c^{i}({\bf x})-$
	$\displaystyle c^{i}({\bf x})-(c_{0}^{i}\ln c_{0}^{i}-c_{0}^{i}))-\mu^{i}(c^{i}({\bf x})-c_{0}^{i})\bigg{]}\bigg{\}},$		(13)

where $\varphi({\bf r})$ is electrostatic potential, $c^{\pm}$ are the concentrations of $\pm$ monovalent electrolyte ions in the bathing solution, with $c_{0}^{\pm}$ being their bulk concentrations, $\mu^{\pm}$ are their chemical potentials, $\varepsilon\varepsilon_{0}$ is the permittivity of water, and $\rho_{surf}$ is the charge density over the bounding surfaces, confining polyelectrolyte. The interconnection between the electrostatic field and polymer density is provided by term $-pe\varphi({\bf x})\rho({\bf x})$ in (Eq. II.3). The ground state dominance approximation suppresses polymer density fluctuations and the Poisson-Boltzmann approximation is expected to be reasonable for a monovalent electrolyte. Finite size chain would lead to a decrease of the energy gap between the ground state and first excited state, thus invalidating the ground state dominance approximation Li et al. (2018), while the presence of multivalent electrolyte would invalidate the saddle-point approximation and thus the Poisson-Boltzmann approximation Naji et al. (2013).

Solving the system of the saddle - point equations (Eqs. B.12,B.16,B.17) we get two separate equations for the electrostatic and configurational degrees of freedom, corresponding to the polymer Poisson-Boltzmann equation Markovich et al. (2021) for the saddle-point electrostatic potential $\varphi({\bf r})$

\varepsilon\varepsilon_{0}\nabla^{2}\varphi({\bf r})=2ec_{0}\sinh(\beta e\varphi({\bf r}))+peN\psi({\bf r})^{2}-\rho_{s}({\bf r})

(14)

as well as the Edwards equation for the saddle-point polymer one-replica density field $\psi({\bf r})$

	$\displaystyle\mathcal{E}_{0}\psi({\bf r})=\bigg{[}-\frac{\ell^{2}}{6}\nabla^{2}-\beta pe\varphi({\bf r})+\frac{1}{2}\frac{1}{\mu+N\psi({\bf r})^{2}}$
	$\displaystyle\bigg{(}\frac{1}{\tilde{v}}+\frac{\kappa}{(1-\kappa)^{2}}\frac{2\mu\psi^{2}({\bf r})+N\psi({\bf r})^{4}}{\mu+N\psi({\bf r})^{2}}\bigg{)}\bigg{]}\psi({\bf r}),$		(15)

corresponding to the monomer density $\rho({\bf r})=N\psi({\bf r})^{2}$ and we introduced te parameter

\kappa=\int d^{3}x\frac{N\psi({\bf x})^{4}}{\mu+N\psi({\bf x})^{2}}~{}~{}~{}~{}{\rm and}~{}~{}~{}~{}\mu=(\beta v_{0}\xi^{2})^{-1}.

Upon substitution of the saddle - point equations (see Appendix, Eq. B.12 - B.17) into (Eq. II.2), we remain with

\begin{split}\mathcal{F}&=\frac{Nl^{2}}{6\beta}\int(\nabla\psi)^{2}d^{3}{\bf x}-\frac{V}{2\beta\tilde{v}}\ln{\mu}+\frac{1}{2\beta}\bigg{(}\frac{\kappa}{1-\kappa}+\ln{(1-\kappa)}\bigg{)}+\\ &+\frac{1}{2\beta\tilde{v}}\int d^{3}{\bf x}\ln{(\mu+N\psi^{2}({\bf x}))}+\frac{\varepsilon\varepsilon_{0}}{2}\int d^{3}{\bf x}(\nabla\varphi({\bf x}))^{2},\end{split}

(16)

which is the form of the free energy that we will use in what follows.

II.4 No-salt regime and planar confinement geometry

Let us now consider the model polyelectrolyte chain confined between two infinite planar surfaces of area $S$ separated by a separation $2al$ , so that $2a$ is the inter – plate separation in the units of the Kuhn segment length $l$ , assuming the limiting salt-free case, corresponding to zero added electrolyte concentration $c_{0}=0$ . This is the simples possible model of a confined polyelectrolyte Podgornik (1991). We will further assume that the homogeneous surface charge density is described by

\rho_{s}({\bf r})=\sigma(\delta(zl-al)+\delta(zl+al)),

(17)

and that

	$\displaystyle\varphi(z)$	$\displaystyle=$	$\displaystyle\varphi({\bf r})$		(18)
	$\displaystyle\Phi(z)$	$\displaystyle=$	$\displaystyle\psi({\bf r})\sqrt{Sl},$		(19)

where due to electroneutrality, $\sigma=\frac{\tau N}{2S}$ , and $\tau=pe$ . We also introduce the dimensionless coordinates ${\bf r}=(xl,yl,zl)$ and stipulate the polymer density normalization $\int_{-a}^{a}dz~{}\Phi^{2}(z)=1$ . With these definitions the Poisson-Boltzmann equation (Eq. (14)) becomes simply the Poisson equation in the form

\frac{d^{2}}{dz^{2}}\varphi(z)=\frac{\sigma l}{\varepsilon\varepsilon_{0}}\bigg{(}2\Phi^{2}(z)-\delta(z-a)-\delta(z+a)\bigg{)},

(20)

and thus has an obvious solution

\varphi(z)=\frac{\sigma l}{\varepsilon\varepsilon_{0}}\bigg{[}\int_{-a}^{a}\!\!\!|z-z^{\prime}|\Phi^{2}(z^{\prime})dz^{\prime}-m(z,a)\bigg{]}.

(21)

with $m(z,a)\equiv\max\big{(}z-a,-a-z,0\big{)}$ . In the same limit the Edwards equation (Eq. (II.3)) can be cast into the following simplified form

\begin{split}\mathcal{E}_{0}\Phi(z)&=\bigg{[}-\frac{1}{6}\frac{d^{2}}{dz^{2}}-\beta\tau\varphi(z)+\frac{1}{2\tilde{v}}\frac{1}{\mu+\frac{2\sigma}{\tau l}\Phi^{2}(z)}+\\ &+\frac{\Phi^{2}(z)}{Sl}\frac{\kappa}{2(1-\kappa)^{2}}\frac{2\mu+\frac{2\sigma}{\tau l}\Phi^{2}(z)}{\big{(}\mu+\frac{2\sigma}{\tau l}\Phi^{2}(z)\big{)}^{2}}\bigg{]}\Phi(z),\end{split}

(22)

where

\kappa=\int_{-a}^{a}dz\frac{2\sigma\Phi^{4}(z)}{\mu\tau l+2\sigma\Phi^{2}(z)}.

(23)

As we consider the bounding plates to be of infinite extension, $S\to\infty$ , the Edwards equation is further simplified and assumes the form

\frac{1}{6}\Phi^{\prime\prime}(z)=-\mathcal{E}_{0}\Phi(z)-\beta\tau\varphi(z)\Phi(z)+\frac{1}{2\tilde{v}}\frac{\Phi(z)}{\mu+\frac{2\sigma}{\tau l}\Phi^{2}(z)}.

(24)

Following the linearization procedure for (Eq. (21)), suggested in Ref. Podgornik (1991), we get for $z\in[0,a]$ an even simpler version of the Edwards equation that can be written as

\frac{1}{6}\Phi^{\prime\prime}(z)=-\bigg{(}\mathcal{E}_{0}+\frac{\beta\sigma\tau l}{\varepsilon\varepsilon_{0}}z\bigg{)}\Phi(z)+\frac{1}{2\tilde{v}}\frac{\Phi(z)}{\mu+\frac{2\sigma}{\tau l}\Phi^{2}(z)},

(25)

coupled with normalization

\begin{split}2\int_{0}^{a}dz~{}\Phi^{2}(z)&=1\end{split}

(26)

and the boundary conditions for the ground state "wave function"

\begin{split}\Phi^{\prime}(0)=0\;\;\;\Phi(a)&=0.\end{split}

(27)

and the first excited states

\begin{split}\Phi(0)=0\;\;\;\Phi(a)&=0.\end{split}

(28)

Furthermore, we introduce the following notations: $\lambda_{B}={\beta e^{2}}/{4\pi\varepsilon_{0}\varepsilon},\nu={N}/{(Sl)}$ and $\alpha=12\pi p^{2}l^{2}\lambda_{B}$ , where $\lambda_{B}$ is the Bjerrum length and $\nu$ is the monomer concentration. The final form of the Edwards equation is then given in the form of a non-linear Schrödinger equation

\Phi^{\prime\prime}(z)=-\bigg{(}6\mathcal{E}_{0}+\alpha\nu~{}z\bigg{)}\Phi(z)+\frac{3}{\tilde{v}}\frac{\Phi(z)}{\mu+\nu\Phi^{2}(z)},

(29)

which is the equation that we will solve numerically in what follows.

II.5 Reduced Free Energy

Using (Eq. (16)) and the results from the previous section, we then get the surface density of the free energy for the salt-free case in the following form

\begin{split}A&=\frac{\mathcal{F}}{S}=\frac{2\sigma}{\tau\beta}\mathcal{E}_{0}-\frac{\sigma}{\tilde{v}\beta\tau}\int_{-a}^{a}dz\frac{\Phi^{2}(z)}{\mu+\frac{2\sigma}{\tau l}\Phi^{2}(z)}+\\ &\frac{1}{2S\beta}\bigg{(}\frac{\kappa}{1-\kappa}+\ln{(1-\kappa)}\bigg{)}+\\ &+\frac{l}{2\tilde{v}\beta}\int_{-a}^{a}dz\ln{\bigg{(}1+\frac{2\sigma}{\tau l\mu}\Phi^{2}(z)\bigg{)}}+\\ &\sigma\bigg{(}\varphi(a)+\int_{-a}^{a}\Phi^{2}(z)\varphi(z)dz\bigg{)}.\end{split}

(30)

After taking the $S\to\infty$ limit and using the notations introduced in (Eq. 29) together with $\tilde{A}=A/k_{B}T\nu l$ , we obtain:

\begin{split}\tilde{A}&=\mathcal{E}_{0}+\frac{\alpha\nu}{12}a+\frac{\alpha\nu}{12}\int_{-a}^{a}dz\int_{-a}^{a}dz^{\prime}|z-z^{\prime}|\Phi^{2}(z^{\prime})\Phi^{2}(z)-\\ &-\int_{-a}^{a}dz\frac{\Phi^{2}(z)}{2\tilde{v}\mu+2\tilde{v}\nu\Phi^{2}(z)}+\frac{1}{2\tilde{v}\nu}\int_{-a}^{a}dz\ln{\big{(}1+\frac{\nu}{\mu}\Phi^{2}(z)\big{)}}.\end{split}

(31)

No additional analytic approximations are feasible and/or appropriate at this point. The integrals indicated above can only be preformed numerically with the numerical solutions of the non-linear Schrödinger equation (Eq. 29).

III Results

First, we compare our numerical results with the previously reported analytical results in Ref. Podgornik (1991), where the confined polyelectrolyte without short-range disordered interactions was considered and the limits allowing the analytical solution were analyzed. The reduced free energy obtained in Ref. Podgornik (1991) can be shown to correspond exactly to (Eq. (31)) with $\xi=0$ , while additional linearization of the third term in terms of $a$ results in a simple analytical solution Podgornik (1991).

We can see from Fig.2 that these additional approximations are not important and the results are effectively the same. The pressure between the two bounding surfaces, i.e., the force per unit area is calculated from

P=-\frac{1}{l}\frac{\partial\tilde{A}}{\partial a}.

(32)

Using this equation we calculate numerically the pressure of the system that is presented in Fig.2, and compared to the analytical results of Ref. Podgornik (1991).

The configurational behavior of the confined polyelectrolyte chain with short–range quenched disorder interactions is governed by the interplay between the inter-plate separation $2a$ and the standard deviation of disorder $\xi$ . The distribution of the dimensionless polyelectrolyte monomers $\Phi^{2}(z)$ is presented in Figs. 3, 4, and 5. For the given value of $\xi$ the behavior of the monomer density distribution on increase of the inter–plate separation is qualitatively the same both for the pure case, corresponding to $\xi=0$ , as well as the disordered cases, $\xi\neq 0$ , up to the critical inter-plate separation $a_{c}$ : upon further increase of the separation the chain goes through a phase transition from the density profile having two maxima separated by the region of lower non-zero density to the two-maxima density profile, with complete depletion of density at the midpoint between the two bounding surfaces. The critical separation $a_{c}$ depends on the magnitude of the disorder dispersion $\xi$ as well as on the charge parameters of the system. The main effect of the quenched sequence disorder at the inter-plate separations corresponding to $a<a_{c}$ is an overall depletion of the monomer density close to the bounding surfaces and accumulation at the relative maximum located at the midpoint between the two bounding surfaces ( Fig. 3) or at the points of the two maxima in case of two-maxima density distribution, leading to an increase of the monomer density in the region of the density maximum ( Fig. 4), i.e., there is an additional localization of the polyelectrolyte chain on increase of the disorder parameter $\xi^{2}$ . However, for separations corresponding to $a>a_{c}$ , there is also a complete depletion of the monomer at the midpoint between the two bounding surfaces, with preferential accumulation actually at the positions of the relative maxima of the monomer density in the left and right halves of the slit ( Fig. 5).

In the presence of disorder the reduced free energy (Eq. 31) is defined by the eigenvalue $\mathcal{E}_{0}$ , corresponding to the "wave function" $\Phi(z)$ of the non-linear Schrödinger equation (Eq. 29). The free energy behavior corresponding to the ground as well as the first excited states is given in Fig. 6, where the singularity of the first derivative of the free energy at $a=a_{c}$ is clearly shown. Thus, the confined polyelectrolyte with short–range sequence disorder exhibits a phase transition, accompanied by a loss of the original ground state at the critical separation $a_{c}$ . In fact, at the critical separation $a_{c}$ the ground state free energy actually coincides with that of the first excited state.

The sequence disorder phase transition at a critical inter-surface separation $a=a_{c}$ is of a first order type, and is accompanied by an abrupt change of the inter-plate pressure (Fig. 7), as well as the vanishing of the monomer density at the midpoint between the surfaces (Fig. 5). The pressure behavior for the pure ( $\xi=0$ ) and disordered ( $\xi=10$ ) cases (Fig. 7) are qualitatively similar below the critical separation $a<a_{c}$ , but diverge sharply for $a>a_{c}$ , when the inter-surface pressure behavior changes drastically from the negative (attractive) pressure in the pure case to the positive (repulsive) pressure in the disordered case. The negative (attractive) pressure for the pure ( $\xi=0$ ) case is the result of the polyelectrolyte bridging Podgornik and Ličer (2006), where part of the chain is partially adsorbed to one plate and part to the other plate, pulling them together because of the chain connectivity. Thus, the bimodal polyelectrolyte chain density profile between the similarly charged plates (see Fig. 4) leads to chain-mediated attractions between them. These bridging interactions stem from the fact that the chain is partially adsorbed to both surfaces and the remaining part in between provides an entropic elastic force between them. This situation is very similar to the case of uncharged polymers, except that the polymer-surface interaction is long-range electrostatics for a polyectrolyte chain, and short-range local for an uncharged polymer chain de Gennes (1987). This furthermore implies that above the critical separation $a_{c}$ the repulsion between the two positively charged plates is no longer compensated by the polyelectrolyte bridging attraction mediated by the negatively charged polymer chain with the short-range disordered interaction. The reason for this qualitative difference can be traced to a complete charge separation of the polyelectrolyte density between the left and right halves of the inter-surface region, promoted by an increased localization of the monomer density as clearly shown in Fig. 5.

The critical separation $a_{c}$ and the very existence of the phase transition strictly depends on the partial charge $pe$ and, consequently on the surface charge density $\sigma$ as summarized in Fig. 8. A decrease in the partial charge $pe$ for the given value of the strength of disorder $\xi$ leads to an increase in the critical separation $a_{c}$ (see in Fig. 8). Without electrostatics ( $p=0$ ), the inter-plate pressure is always positive (repulsive) without the first-order phase transition. On the other hand, the effect of the strength of disorder $\xi$ on the phase behavior of the system is not unique. We have no phase transition for small enough values of $\xi$ , while for higher values of $\xi$ , as well as dependent on the partial charge $pe$ , the system exhibits a first-order phase transition. In general, the increase in $\xi$ , for a given value of partial charge $pe$ , leads to an increase in the critical separation $a_{c}$ (see in Fig. 8). The values of the critical separation $a_{c}$ presented in Fig. 8 for the different values of $\xi$ and $pe$ parameters are summarized in Table 1.

Table 1: The dependence of the critical separation

a_{c}=a_{c}(p,\xi)

on parameters

\xi

and

p

$p$	$\xi$	$a_{c}$
1	10	1.410
1	16	1.670
0.8	10	1.525
0.6	10	1.970
0.6	13	1.975

IV Discussion

We developed a self-consistent field (SCF) formalism describing behavior of the polyelectrolyte chain confined between two charged impenetrable bounding surfaces. There is no added salt present in the system and the charges on the polyelectrolyte chain are compensated solely by the charges on the surfaces. The identity of the monomers along the polyelectrolyte chain is described by a random variable, $\xi_{\tau}$ , characterized by a Gaussian distribution of zero mean and variance $\xi^{2}$ . The quenched distribution of monomers along the contour only affects the short range interactions, while the long-range interactions are Coulombic, corresponding to a fixed charge per monomer.

To calculate the free energy and density profile of the monomers, we have used the replica approach and obtained a phase transition of the first order at a critical separation $a_{c}$ for the replica-symmetric solution. We have considered the polyelectrolyte chain in no-added electrolyte regime, confined between two oppositely charged, impenetrable planar surfaces with the specific aim to asses the contribution of the quenched sequence disorder and the short-range interactions between the monomers. The obtained results can be interpreted in terms of a “quantum-like" Hamiltonian ${\cal{\hat{H}}}_{n}=-\frac{\ell^{2}}{6}\sum_{a}\nabla_{a}^{2}+\imath\sum_{a}{\hat{\rho}}_{a}({\bf x}_{a})+\imath\sum_{a<b}{\hat{q}}_{ab}({\bf x}_{a},{\bf x}_{b})$ , describing an effective particle in the random external field $\imath\sum_{a}{\hat{\rho}}_{a}({\bf x}_{a})+\imath\sum_{a<b}{\hat{q}}_{ab}({\bf x}_{a},{\bf x}_{b})$ . The confinement of the polyelectrolyte chain inside the slit with charged walls would thus correspond to the quantum-mechanical particle localized in a potential well. The behavior of the polymer density $\Phi(z)^{2}$ presented in Figs.(3, Eq. 4) exhibits the additional localization induced by sequence disorder and consequently by the random effective field in the “quantum-like" Hamiltonian ${\cal{\hat{H}}}_{n}$ , akin to the Anderson localization of the polymer chain in a random medium, as was e.g. demonstrated in Ref. Shiferaw and Goldschmidt (2001). The essence of Anderson localization phenomenon Anderson (1958); Crisanti et al. (2012) consists of the fact that under certain conditions the wave function of an electron in a random potential does not spread over the whole space, but is localized in a finite region. In the self-consistent field approximation the confined polyelectrolyte is already localized, but the disorder in the sequence creates an effective random potential that leads to an additional, disorder driven localization of the "wave function" $\Phi(z)$ .

Besides the additional, disorder driven localization, the sequence disorder also results in a phase transition of a first order type, accompanied by an abrupt change of the inter-surface pressure and monomer density (Fig. 5), on increase of the inter-surface separation beyond the critical separation $a_{c}$ . The pressure behavior for the pure and disordered cases are actually qualitatively similar below the critical separations, $a<a_{c}$ , but at separations $a>a_{c}$ , the inter-surface pressure behavior changes drastically, flipping from a negative pressure in the pure case to a positive one in the disordered case. Above the critical separation $a_{c}$ the repulsion between the two positively charged plates ceases to be mediated by the polyelectrolyte chain bridging mechanism, due to the complete separation of the polyelectrolyte density between the left and right halves of the slit as a consequence of the disorder-driven localization of the monomer density.

The work presented above introduces an interesting and important new feature in the phenomenology of polyelectrolyte-driven interactions in charged systems, where the polymer chains exhibit a "chemical" sequence disorder coupled to short-range interactions.

Acknowledgements.

The SCS of Armenia supported this work, grant No. 22AA-1C023 (V.S., Y.M.), grant No. 21AG-1C038 (V.S.), grant No. 21T-1F307 (Y.M.). A.B. and Y.M. acknowledge the partial financial support from Erasmus+ Project No. (2023-1-SI01-KA171-HED-000138423). R.P. acknowledges funding from the Key Project No. 12034019 of the National Natural Science Foundation of China.

author declarations

Conflict of Interest

The authors have no conflicts to disclose.

Author Contribution

V. Stepanyan: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Soft- ware (equal); Validation (equal); Visualization (equal).

A. Badasyan: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Validation (equal); Writing – original draft (equal).

V. Morozov: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (equal).

Y. Mamasakhlisov: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – original draft (equal)

R. Podgornik: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (equal); Validation (equal); Writing – original draft (equal).

Data Availability Statement

The data that support the findings of this study are available within the article.

References

Leckband and Israelachvili (2001) D. Leckband and J. Israelachvili, Quarterly Reviews of Biophysics 34, 105 (2001).
Rubinstein and Colby (2003) M. Rubinstein and R. H. Colby, Polymer Physics (Oxford University Press, New York, 2003).
Dobrynin and Rubinstein (2005) A. V. Dobrynin and M. Rubinstein, Progress in Polymer Science 30, 1049 (2005).
Muthukumar (2023) M. Muthukumar, Physics of Charged Macromolecules: Synthetic and Biological Systems (Cambridge University Press, 2023).
French et al. (2010) R. H. French, V. A. Parsegian, R. Podgornik, R. F. Rajter, A. Jagota, J. Luo, D. Asthagiri, M. K. Chaudhury, Y.-m. Chiang, S. Granick, S. Kalinin, M. Kardar, R. Kjellander, D. C. Langreth, J. Lewis, S. Lustig, D. Wesolowski, J. S. Wettlaufer, W.-Y. Ching, M. Finnis, F. Houlihan, O. A. von Lilienfeld, C. J. van Oss, and T. Zemb, Rev. Mod. Phys. 82, 1887 (2010).
Elias (2012) H.-G. Elias, Macromolecules· 1: Volume 1: Structure and Properties (Springer Science & Business Media, 2012).
Smiatek (2020) J. Smiatek, Molecules 25 (2020), 10.3390/molecules25071661.
Lu et al. (2015) B.-S. Lu, A. Naji, and R. Podgornik, The Journal of Chemical Physics 142, 214904 (2015), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4921892/13254056/214904_1_online.pdf .
Borkovec et al. (2001) M. Borkovec, B. Jönsson, and G. J. Koper, Ionization processes and proton binding in polyprotic systems: Small molecules, proteins, interfaces, and polyelectrolytes (Springer, 2001).
Blossey and Podgornik (2022) R. Blossey and R. Podgornik, Europhysics Letters 139, 27002 (2022).
Bloomfield et al. (2000) V. A. Bloomfield, D. M. Crothers, and I. Tinoco, Nucleic acids: structures, properties, and functions (University Science Books, 2000).
Finkelstein and Ptitsyn (2016) A. V. Finkelstein and O. Ptitsyn, Protein physics: a course of lectures (Elsevier, 2016).
Fallmann et al. (2017) J. Fallmann, S. Will, J. Engelhardt, B. Grüning, R. Backofen, and P. F. Stadler, Journal of biotechnology 261, 97 (2017).
Poblete et al. (2021) S. Poblete, B. Anže, M. Kanduč, R. Podgornik, and H. V. Guzman, ACS Omega 6, 32823 (2021).
Twarock et al. (2018) R. Twarock, R. J. Bingham, E. C. Dykeman, and P. G. Stockley, Current opinion in virology 31, 74 (2018).
Adhikari et al. (2020) P. Adhikari, N. Li, M. Shin, N. F. Steinmetz, R. Twarock, R. Podgornik, and W.-Y. Ching, Physical Chemistry Chemical Physics 22, 18272 (2020).
Borukhov et al. (1999) I. Borukhov, D. Andelman, and H. Orland, J. Phys. Chem. B 103, 5042 (1999).
Shafir et al. (2003) A. Shafir, D. Andelman, and R. R. Netz, J. Chem. Phys. 119, 2355 (2003).
Siber and Podgornik (2008) A. Siber and R. Podgornik, Phys. Rev. E 78, 051915 (2008).
Podgornik (1991) R. Podgornik, J. Phys. Chem. 95, 5249 (1991).
Andelman (1995) D. Andelman, Structure and Dynamics of Membranes Chap. 12, edited by E. S. R. Lipowsky, Vol. 1B (Elsevier, Amsterdam, Handbook of Biological Physics, 1995).
Budkov and Kalikin (2023) Y. A. Budkov and N. N. Kalikin, Phys. Rev. E 107, 024503 (2023).
Muthukumar (1989) M. Muthukumar, The Journal of Chemical Physics 90, 4594 (1989).
Hribar-Lee et al. (2011) B. Hribar-Lee, M. Lukšič, and V. Vlachy, Annu. Rep. Prog. Chem., Sect. C: Phys. Chem. 107, 14 (2011).
Bratko and Chakraborty (1995) D. Bratko and A. Chakraborty, Physical Review E 51, 5805 (1995).
Binder and Young (1986) K. Binder and A. P. Young, Rev. Mod. Phys. 58, 801 (1986).
Parisi (2023) G. Parisi, Rev. Mod. Phys. 95, 030501 (2023).
Garel et al. (1997) T. Garel, H. Orland, and E. Pitard, Spin glasses and random fields, edited by A. Young (World Scientific, 1997).
Li et al. (2018) S. Li, H. Orland, and R. Zandi, Journal of Physics: Condensed Matter 30, 144002 (2018).
Naji et al. (2013) A. Naji, M. Kanduč, J. Forsman, and R. Podgornik, The Journal of Chemical Physics 139, 150901 (2013).
Markovich et al. (2021) T. Markovich, D. Andelman, and R. Podgornik, in Handbook of lipid membranes (CRC Press, 2021) pp. 99–128.
Podgornik and Ličer (2006) R. Podgornik and M. Ličer, Current Opinion in Colloid & Interface Science 11, 273 (2006).
de Gennes (1987) P. de Gennes, Advances in Colloid and Interface Science 27, 189 (1987).
Shiferaw and Goldschmidt (2001) Y. Shiferaw and Y. Y. Goldschmidt, Phys. Rev. E 63, 051803 (2001).
Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
Crisanti et al. (2012) A. Crisanti, G. Paladin, and A. Vulpiani, Products of Random Matrices: in Statistical Physics, Springer Series in Solid-State Sciences (Springer Berlin Heidelberg, 2012).

Appendices

Appendix A The free energy functional

The quenched free energy (Eq. 6) can be estimated using the replica trick Binder and Young (1986)

-\beta{\cal F}=\lim_{n\rightarrow 0}\frac{\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}-1}{n},

(A.1)

where $\beta=(k_{B}T)^{-1}$ . We therefore need to calculate the $n$ -replica partition function as

	$\displaystyle\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}=$
	$\displaystyle\int\,{\cal D}{\bf r}e^{-\frac{3}{2\ell^{2}}\sum_{a=1}^{n}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf r}^{a}(\tau))^{2}-\beta\sum_{a=1}^{n}V\{{\bf r}^{a}\}}\times$
	$\displaystyle\int\,{\cal D}\xi{\cal P}\{\xi\}e^{-\frac{\beta v_{0}}{2}\int_{0}^{N}\,d\tau\int_{0}^{N}\,d\tau^{\prime}\xi_{\tau}\xi_{\tau^{\prime}}\sum_{a=1}^{n}\delta({\bf r}^{a}(\tau)-{\bf r}^{a}(\tau^{\prime}))}.~{}~{}~{}~{}~{}$		(A.2)

The average over the distribution function (Eq. 1) in the equation (Eq. A) is transformed as

	$\displaystyle\int\,{\cal D}\xi{\cal P}\{\xi\}e^{-\frac{\beta v_{0}}{2}\int_{0}^{N}\,d\tau\int_{0}^{N}\,d\tau^{\prime}\xi_{\tau}\xi_{\tau^{\prime}}\sum_{a=1}^{n}\delta({\bf r}^{a}(\tau)-{\bf r}^{a}(\tau^{\prime}))}=$
	$\displaystyle\int\,{\cal D}\xi e^{-\frac{1}{2\xi^{2}}\int_{0}^{N}\,d\tau\xi_{\tau}^{2}}e^{-\frac{\beta v_{0}}{2}\int\,d^{3}{\bf x}\sum_{a=1}^{n}m_{a}({\bf x})^{2}},~{}~{}~{}~{}$		(A.3)

where the field $m_{a}({\bf x})=\int_{0}^{N}\,d\tau\xi_{\tau}\delta({\bf x}-{\bf r}^{a}(\tau))$ , and $v_{0}$ is the interaction strength parameter. The r.h.s. of the equation (Eq. A) is linearized using the Hubbard-Stratonovich transformation as

	$\displaystyle e^{-\frac{\beta v_{0}}{2}\int\,d^{3}{\bf x}\sum_{a=1}^{n}m_{a}({\bf x})^{2}}=e^{-\frac{nV}{2}\ln(2\pi\beta v_{0})}$
	$\displaystyle\int\,{\cal D}{\Pi}e^{-\frac{1}{2\beta v_{0}}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})^{2}+\imath\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})m_{a}({\bf x})}.$		(A.4)

Insertion of (Eqs. A,A) into the (Eq. A) upon averaging over the $\{\xi\}$ variables yields the $n$ -replica partition function as

	$\displaystyle\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}\propto e^{-\frac{nV}{2\tilde{v}}\ln(2\pi\beta v_{0})}\int\,{\cal D}{\Pi}e^{-\frac{1}{2\beta v_{0}}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})^{2}}\times$
	$\displaystyle\int\,{\cal D}{\bf r}e^{-\frac{3}{2\ell^{2}}\sum_{a=1}^{n}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf r}^{a}(\tau))^{2}-\beta\sum_{a=1}^{n}V\{{\bf r}^{a}\}}\times$
	$\displaystyle e^{-\frac{\xi^{2}}{2}\sum_{a,b}\int_{0}^{N}\,d\tau\Pi_{a}({\bf r}^{a}(\tau))\Pi_{b}({\bf r}^{b}(\tau))}.$		(A.5)

For further consideration it is useful to introduce two order parameters, viz., the inter-replica overlap $q_{ab}({\bf x},{\bf x}^{\prime})$ , with $a<b$ , and the density of $a$ -th replica monomers $\rho_{a}({\bf x})$ as

	$\displaystyle q_{ab}({\bf x},{\bf x}^{\prime})=\int_{0}^{N}\,d\tau\delta({\bf x}-{\bf r}^{a}(\tau))\delta({\bf x}^{\prime}-{\bf r}^{b}(\tau))$
	$\displaystyle\rho_{a}({\bf x})=\int_{0}^{N}\,d\tau\delta({\bf x}-{\bf r}^{a}(\tau)).$		(A.6)

Using these order parameters (Eq. A), the $n$ -replica partition function (Eq. A) is transformed as

	$\displaystyle\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}\propto\int\!\!{\cal D}{\varphi}\!\!\int\!\!{\cal D}{c^{\pm}}\!\!\int\!\!{\cal D}{\Pi}~{}e^{-\frac{1}{2\beta v_{0}}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})^{2}}\times$
	$\displaystyle\int\,{\cal D}{\bf r}e^{-\frac{3}{2\ell^{2}}\sum_{a=1}^{n}\int_{0}^{N}\,d\tau(\partial_{\tau}{\bf r}^{a}(\tau))^{2}}\times$
	$\displaystyle e^{-\beta\sum_{a=1}^{n}W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})-\frac{\xi^{2}}{2}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})^{2}\rho_{a}({\bf x})}\times$
	$\displaystyle e^{-\xi^{2}\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf x}^{\prime}\Pi_{a}({\bf x})\Pi_{b}({\bf x}^{\prime})q_{ab}({\bf x},{\bf x}^{\prime})},~{}~{}~{}~{}~{}$		(A.7)

where the electrostatic part of the free energy $W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})$ is written as a sum of electrostatic energies of polyelectrolyte segments, fixed charges residing on the surface, confining polyelectrolyte and charges of salt ions, as well as the entropies of the salt ions. As shown in Borukhov et al. (1999); Shafir et al. (2003); Siber and Podgornik (2008),

	$\displaystyle W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})=\int\,d^{3}{\bf x}\bigg{\{}-\frac{\epsilon\epsilon_{0}}{2}(\nabla\varphi_{a}({\bf x}))^{2}+$
	$\displaystyle\varphi_{a}({\bf x})\bigg{[}ec_{a}^{+}({\bf x})-ec_{a}^{-}({\bf x})-pe\rho_{a}({\bf x})+\rho_{surf}({\bf x})\bigg{]}+$
	$\displaystyle\sum_{i=\pm}\bigg{[}k_{B}T(c_{a}^{i}({\bf x})\ln c_{a}^{i}({\bf x})-$
	$\displaystyle c_{a}^{i}({\bf x})-(c_{0}^{i}\ln c_{0}^{i}-c_{0}^{i}))-\mu^{i}(c_{a}^{i}({\bf x})-c_{0}^{i})\bigg{]}\bigg{\}}.$		(A.8)

Here $c^{\pm}$ are the concentrations of $\pm$ monovalent salt ions, with $c_{0}^{\pm}$ being their bulk concentrations, and $\mu^{\pm}$ their chemical potentials, $\varepsilon\varepsilon_{0}$ is the permittivity of water, and $\rho_{surf}$ is the charge distribution over the bounding surfaces, confining polyelectrolyte.

Introducing the Lagrange multipliers ${\hat{q}}_{ab}({\bf x},{\bf x}^{\prime})$ and the ${\hat{\rho}}_{a}({\bf x})$ conjugated with the order parameters (A), the $n$ -replica partition function reads Garel et al. (1997)

	$\displaystyle\langle{Z\{\xi\}^{n}}\rangle_{{\cal P}}\propto e^{-\frac{nV}{2\tilde{v}}\ln(2\pi\beta v_{0})}\int\,{\cal D}{\Pi}\,{\cal D}{\rho}\,{\cal D}{{\hat{\rho}}}\,{\cal D}{q}\,{\cal D}{\hat{q}}$
	$\displaystyle e^{G(\Pi,\rho,{\hat{\rho}},q,{\hat{q}})+\ln\zeta({\hat{\rho}},{\hat{q}})},$		(A.9)

where

	$\displaystyle G(\Pi,\rho,{\hat{\rho}},q,{\hat{q}},\varphi,c^{\pm})=-\frac{1}{2\beta v_{0}}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})^{2}-$
	$\displaystyle-\beta\sum_{a=1}^{n}W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})-\frac{\xi^{2}}{2}\sum_{a}\int\,d^{3}{\bf x}\Pi_{a}({\bf x})^{2}\rho_{a}({\bf x})-$
	$\displaystyle-\xi^{2}\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf x}^{\prime}\Pi_{a}({\bf x})\Pi_{b}({\bf x}^{\prime})q_{ab}({\bf x},{\bf x}^{\prime})+$
	$\displaystyle+\imath\sum_{a}\int\,d^{3}{\bf x}\rho_{a}({\bf x}){\hat{\rho}}_{a}({\bf x})+\imath\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf x}^{\prime}q_{ab}({\bf x},{\bf x}^{\prime}){\hat{q}}_{ab}({\bf x},{\bf x}^{\prime})$		(A.10)

and

	$\displaystyle\ln\zeta({\hat{\rho}},{\hat{q}})=-N\min_{\psi}\bigg{\{}\int\,d^{3n}{\bf x}\psi({\bf x}_{1},...,{\bf x}_{n})\times$
	$\displaystyle\bigg{(}-\frac{\ell^{2}}{6}\sum_{a}\nabla_{a}^{2}+\imath\sum_{a}{\hat{\rho}}_{a}({\bf x}_{a})+\imath\sum_{a<b}{\hat{q}}_{ab}({\bf x}_{a},{\bf x}_{b})\bigg{)}\times$
	$\displaystyle\psi({\bf x}_{1},...,{\bf x}_{n})-\mathcal{E}_{0}\bigg{(}\int\,d^{3n}{\bf x}\psi({\bf x}_{1},...,{\bf x}_{n})^{2}-1\bigg{)}\bigg{\}}.$		(A.11)

Here $\mathcal{E}_{0}$ is the ground state energy and $\psi({\bf x}_{1},...,{\bf x}_{n})$ is the corresponding eigenfunction of the “quantum-like" Hamiltonian ${\cal{\hat{H}}}_{n}=-\frac{\ell^{2}}{6}\sum_{a}\nabla_{a}^{2}+\imath\sum_{a}{\hat{\rho}}_{a}({\bf x}_{a})+\imath\sum_{a<b}{\hat{q}}_{ab}({\bf x}_{a},{\bf x}_{b})$ . Since $v_{0}>0$ , we have no trend of segregation between different kinds of monomers and we expect large fluctuations of the fields $\{\Pi\}$ . Thus, we can consider the saddle-point approximation in the equation (Eq. A) over the fields $\{\rho,{\hat{\rho}},q,{\hat{q}}\}$ only and to integrate over the fields $\{\Pi\}$ . Let us introduce $g(\rho,{\hat{\rho}},q,{\hat{q}},\varphi,c^{\pm})=\ln\int\,{\cal D}{\Pi}~{}e^{G(\Pi,\rho,{\hat{\rho}},q,{\hat{q}})}$ , where

	$\displaystyle\int\,{\cal D}{\Pi}~{}e^{G(\Pi,\rho,{\hat{\rho}},q,{\hat{q}})}=e^{-\beta\sum_{a=1}^{n}W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})}\times$
	$\displaystyle e^{\imath\sum_{a}\int\,d^{3}{\bf x}\rho_{a}({\bf x}){\hat{\rho}}_{a}({\bf x})+\imath\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf x}^{\prime}q_{ab}({\bf x},{\bf x}^{\prime}){\hat{q}}_{ab}({\bf x},{\bf x}^{\prime})}\times$
	$\displaystyle\int\,{\cal D}{\Pi}e^{-\frac{1}{2}\sum_{a,b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf x}^{\prime}P_{ab}({\bf x},{\bf x}^{\prime})\Pi_{a}({\bf x})\Pi_{b}({\bf x}^{\prime})},$		(A.12)

and $P_{ab}({\bf x},{\bf x}^{\prime})=\delta_{ab}\delta({\bf x}-{\bf x}^{\prime})[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})]+\xi^{2}(1-\delta_{ab})q_{ab}({\bf x},{\bf x}^{\prime})$ . Thus, $g(\rho,{\hat{\rho}},q,{\hat{q}})$ can be obtained as

	$\displaystyle g(\rho,{\hat{\rho}},q,{\hat{q}},\varphi,c^{\pm})=-\frac{1}{2}\ln\det\frac{{\hat{P}}}{2\pi}-$
	$\displaystyle\beta\sum_{a=1}^{n}W_{el}(\rho_{a},\varphi_{a},c_{a}^{\pm})+\imath\sum_{a}\int\,d^{3}{\bf x}\rho_{a}({\bf x}){\hat{\rho}}_{a}({\bf x})+$
	$\displaystyle\imath\sum_{a<b}\int\,d^{3}{\bf x}\int\,d^{3}{\bf x}^{\prime}q_{ab}({\bf x},{\bf x}^{\prime}){\hat{q}}_{ab}({\bf x},{\bf x}^{\prime}).$		(A.13)

Let us consider the term $\ln\det$ in the equation (Eq. A). These calculations may be carried out by employing properties of block-matrices and the operator algebra defined over the Hilbert space $\{|{\bf x}\rangle\}$ . We shall use the compact notation by defining the operators $\langle{\bf x}|{\hat{A}}_{ab}|{\bf x}^{\prime}\rangle=\delta_{ab}\delta({\bf x}-{\bf x}^{\prime})[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})]$ as the $n\times n$ operator matrix and $\langle{\bf x}|{\hat{P}}_{ab}|{\bf x^{\prime}}\rangle=P_{ab}({\bf x},{\bf x^{\prime}})$ as an element of the $n\times n$ operator matrix

{\hat{P}}={\hat{A}}+\xi^{2}{\hat{B}},

(A.14)

where $\langle{\bf x}|{\hat{B}}_{ab}|{\bf x^{\prime}}\rangle=(1-\delta_{ab})q_{ab}({\bf x},{\bf x}^{\prime})$ is the $n\times n$ operator matrix. Therefore

	$\displaystyle\ln\det{\hat{P}}={\rm Tr}\ln{\hat{P}}={\rm Tr}\ln\bigg{\{}{\hat{A}}+\xi^{2}{\hat{B}}\bigg{\}}={\rm Tr}\ln{\hat{A}}+$
	$\displaystyle{\rm Tr}\ln\bigg{\{}{\hat{1}}_{n}\otimes{\hat{1}}+\xi^{2}{\hat{A}}^{-1}{\hat{B}}\bigg{\}},$		(A.15)

where $({\hat{1}}_{n})_{ab}=\delta_{ab}$ and $\langle{\bf x}|{\hat{1}}|{\bf x}^{\prime}\rangle=\delta({\bf x}-{\bf x}^{\prime})$ . Calculation of the first term on the r.h.s. of the equation (Eq. A) is straightforward

	$\displaystyle{\rm Tr}\ln{\hat{A}}=\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k}\sum_{a}\int\,d^{3}{\bf x}\langle{\bf x}\|({\hat{A}}-{\hat{1}})^{k}_{aa}\|{\bf x}\rangle=$
	$\displaystyle=\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k}\sum_{a}\int\,d^{3}{\bf x}\delta({\bf x}-{\bf x})[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})-1]^{k}=$
	$\displaystyle=\sum_{a}\int\,d^{3}{\bf x}\delta({\bf x}-{\bf x})\ln[\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})]\approx$
	$\displaystyle\approx\frac{1}{\tilde{v}}\sum_{a}\int\,d^{3}{\bf x}\ln\bigg{(}\frac{1}{\beta v_{0}}+\xi^{2}\rho_{a}({\bf x})\bigg{)},$		(A.16)

using $\delta(0)\approx\frac{1}{\tilde{v}}$ , where $\tilde{v}$ is volume of the monomer.

Appendix B Replica symmetric solution

In the replica-symmetric case we will consider only the Hartree-like replica-symmetric wave functions

\psi({\bf x}_{1},...,{\bf x}_{n})=\prod_{a=1}^{n}\psi({\bf x}_{a}),

(B.1)

where $\psi({\bf x}_{a})$ is the one-replica wave function. In the replica symmetric regime the expression (Eq. II.2) may be written (using (Eq. A)) as

	$\displaystyle-\beta{\mathcal{F}}=\lim_{n\rightarrow 0}\frac{1}{n}\bigg{[}-\frac{n}{2\tilde{v}}\int\,d^{3}{\bf x}\ln\bigg{(}1+\frac{\rho_{a}({\bf x})}{\mu}\bigg{)}-\frac{1}{2}{\rm Tr}\ln\bigl{\{}I+\hat{C}\bigr{\}}-$
	$\displaystyle-\beta nW_{el}+in\langle\rho\|\hat{\rho}\rangle+i\frac{n(n-1)}{2}\langle q\|\hat{q}\rangle+\ln\zeta\bigg{]},$		(B.2)

where $\hat{C}=\xi^{2}\hat{A}^{-1}\hat{B}$ and $\mu=(\beta v_{0}\xi^{2})^{-1}$ . With the replica symmetry we have $\hat{C}=\hat{T}\otimes\hat{C^{\prime}}$ where $\hat{T}_{ab}=1-\delta_{ab}$ and $\langle x|\hat{C^{\prime}}|x^{\prime}\rangle=\frac{q(x,x^{\prime})}{\mu+\rho(x)}$ . Taking the limit we get the reduced free energy

	$\displaystyle f=-\beta{\mathcal{F}}=-\frac{1}{2\tilde{v}}\int\,d^{3}{\bf x}\ln\bigg{(}1+\frac{\rho_{a}({\bf x})}{\mu}\bigg{)}-\frac{\Xi}{2}-$
	$\displaystyle-\beta W_{el}+i\langle\rho\|\hat{\rho}\rangle-\frac{i}{2}\langle q\|\hat{q}\rangle+\ln\zeta,$		(B.3)

where we defined

\Xi\equiv\lim_{n\to 0}\frac{1}{n}{\rm Tr}\ln\bigl{\{}I+\hat{C}\bigr{\}}=\lim_{n\to 0}\frac{1}{n}\sum_{k=1}^{\infty}\frac{(-1)^{k+1}}{k}{\rm Tr}(\hat{C}^{k}).

(B.4)

Thus, using ${\rm Tr}(\hat{T}\otimes\hat{C^{\prime}})={\rm Tr}(\hat{T}){\rm Tr}(\hat{C^{\prime}})$ , the terms of the series are written

	$\displaystyle{\rm Tr}(\hat{C})^{k}=T_{k}\int d^{3k}x\frac{q(x^{(1)},x^{(2)})q(x^{(2)},x^{(3)})...q(x^{(k)},x^{(1)})}{(\mu+\rho(x^{(1)}))(\mu+\rho(x^{(2)}))..(\mu+\rho(x^{(k)}))}$
	$\displaystyle T_{k}=\sum_{a_{1},..,a_{k}}\bigg{[}\prod_{i=2}^{k}\biggl{\{}(1-\delta_{a_{i},a_{i-1}})\biggr{\}}(1-\delta_{a_{k},a_{1}})\bigg{]}$		(B.5)

Where $T_{k}={\rm Tr}(\hat{T}^{k})$ . To calculate $T_{k}$ first notice the following recursive equation

T_{k}=T_{k-1}(n-2)+T_{k-2}(n-1)

(B.6)

with the initial conditions $T_{2}=n(n-1)$ and $T_{3}=(n-2)T_{2}$ . Let us define a generating function $P(x)=\sum_{i=0}^{\infty}x^{i}T_{i+2}$ . Using

	$\displaystyle x(n-2)P(x)=\sum_{i=1}^{\infty}x^{i}T_{i+1}(n-2)$		(B.7)
	$\displaystyle x^{2}(n-1)P(x)=\sum_{i=2}^{\infty}x^{i}T_{i}(n-1)$

we can furthermore write

P(x)\bigg{(}1-(n-2)x-(n-1)x^{2}\bigg{)}=T_{2}=n(n-1)

(B.8)

and

P(x)=\sum_{i=0}^{\infty}x^{i}\bigg{(}(n-1)^{i+2}-(-1)^{i+1}(n-1)\bigg{)}=\sum_{i=0}^{\infty}x^{i}T_{i+2}.

(B.9)

Thus, we can take the limit

	$\displaystyle\lim_{n\to 0}\frac{T_{k}}{n}=\lim_{n\to 0}\frac{(n-1)^{k}-(-1)^{k-1}(n-1)}{n}=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
	$\displaystyle=\lim_{n\to 0}\bigg{(}k(n-1)^{k-1}-(-1)^{k-1}\bigg{)}=(-1)^{k-1}(k-1)~{}~{}~{}~{}~{}~{}$		(B.10)

Substituting (Eq. B) into $\Xi$ defined by (Eq. B.4) we get

\Xi=\sum_{k=2}^{\infty}\frac{k-1}{k}\int d^{3k}x\frac{q({\bf x}^{(1)},{\bf x}^{(2)})q({\bf x}^{(2)},{\bf x}^{(3)})...q({\bf x}^{(k)},{\bf x}^{(1)})}{(\mu+\rho({\bf x}^{(1)}))(\mu+\rho({\bf x}^{(2)}))..(\mu+\rho({\bf x}^{(k)}))}

(B.11)

The saddle - point equations can be written in the following form

0=\frac{\delta f}{\delta\hat{\rho}({\bf r})}\Rightarrow\rho({\bf r})=N\psi({\bf r})^{2}

(B.12)

0=\frac{\delta f}{\delta\hat{q}({\bf r},{\bf r}^{\prime})}\Rightarrow q({\bf r},{\bf r}^{\prime})=N\psi({\bf r})^{2}\psi({\bf r}^{\prime})^{2}

(B.13)

0=\frac{\delta f}{\delta\rho({\bf r})}\Rightarrow\mathrm{i}\hat{\rho}({\bf r})=-\beta pe\varphi({\bf r})+\frac{1}{2\tilde{v}}\frac{1}{\mu+\rho(\bf r)}+\frac{1}{2}\frac{\delta\Xi}{\delta\rho}

(B.14)

0=\frac{\delta f}{\delta q({\bf r},{\bf r}^{\prime})}\Rightarrow\mathrm{i}\hat{q}({\bf r},{\bf r}^{\prime})=-\frac{\delta\Xi}{\delta q}

(B.15)

0=\frac{\delta f}{\delta\varphi({\bf r})}\Rightarrow\varepsilon\varepsilon_{0}\nabla^{2}\varphi({\bf r})=-(ec^{+}-ec^{-}-pe\rho({\bf r})+\rho_{s}({\bf r}))

(B.16)

0=\frac{\delta f}{\delta c^{\pm}({\bf r})}\Rightarrow c^{\pm}({\bf r})=c_{0}\exp({\mp\beta e\varphi({\bf r})})

(B.17)

	$\displaystyle 0=\frac{\delta f}{\delta\psi({\bf r})}=-N\frac{\delta}{\delta\psi({\bf r})}\bigg{\{}\frac{\ell^{2}}{6}\langle\nabla\psi\|\nabla\psi\rangle+$
	$\displaystyle\mathrm{i}\langle\psi^{2}\|\hat{\rho}\rangle--\frac{\mathrm{i}}{2}\langle\psi^{2}\|\hat{q}\|\psi^{2}\rangle-\mathcal{E}_{0}(\langle\psi\|\psi\rangle-1)\bigg{\}}=$
	$\displaystyle-N\bigg{\{}-\frac{\ell^{2}}{3}\nabla^{2}\psi({\bf r})+2\mathrm{i}\psi({\bf r})\hat{\rho}({\bf r})-\mathrm{i}\psi({\bf r})\int\,d^{3}{\bf x}\hat{q}({\bf r},{\bf x})\psi({\bf x})^{2}-$
	$\displaystyle\mathrm{i}\psi({\bf r})\int\,d^{3}{\bf x}\hat{q}({\bf x},{\bf r})\psi({\bf x})^{2}-2\mathcal{E}_{0}\psi({\bf r})\bigg{\}}$		(B.18)

where

\begin{split}&\frac{\delta\Xi}{\delta\rho({\bf x})}=-\sum_{k=2}^{\infty}\frac{k-1}{(\mu+\rho({\bf x}))}\\ &\int d^{3k-3}{\bf x}\frac{q({\bf x},{\bf x}^{(1)})q({\bf x}^{(1)},{\bf x}^{(2)})...q({\bf x}^{(k-1)},{\bf x})}{(\mu+\rho({\bf x}))(\mu+\rho({\bf x}^{(1)}))..(\mu+\rho({\bf x}^{(k-1)}))}\end{split}

(B.19)

\begin{split}&\frac{\delta\Xi}{\delta q({\bf x},{\bf x}^{\prime})}=\sum_{k=2}^{\infty}\frac{k-1}{(\mu+\rho({\bf x}))}\\ &\int d^{3k-6}{\bf x}\frac{q({\bf x}^{\prime},{\bf x}^{(1)})q({\bf x}^{(1)},{\bf x}^{(2)})...q({\bf x}^{(k-2)},{\bf x})}{(\mu+\rho({\bf x}^{\prime}))(\mu+\rho({\bf x}^{(1)}))..(\mu+\rho({\bf x}^{(k-2)}))}\end{split}

(B.20)

By further defining

\kappa\equiv\int d^{3}x\frac{N\psi({\bf x})^{4}}{\mu+N\psi({\bf x})^{2}}

and using (Eq. B.12, B.13) we get

\frac{\delta\Xi}{\delta\rho({\bf r})}=-\frac{\kappa}{(1-\kappa)^{2}}\frac{N\psi({\bf r})^{4}}{(\mu+N\psi({\bf r})^{2})^{2}}

(B.21)

\frac{\delta\Xi}{\delta q({\bf r},{\bf r}^{\prime})}=\frac{1}{(1-\kappa)^{2}}\frac{N\psi({\bf r})^{2}\psi({\bf r}^{\prime})^{2}}{(\mu+N\psi({\bf r})^{2})(\mu+N\psi({\bf r}^{\prime})^{2})}

(B.22)

Substituting equations (Eqs. B.14,B.15) into the equation (Eq. B) we finally obtain

	$\displaystyle\frac{\ell^{2}}{6}\nabla^{2}\psi({\bf r})=\psi({\bf r})\bigg{[}-\beta pe\varphi({\bf r})+\frac{1}{2\tilde{v}}\frac{1}{\mu+\rho(\bf r)}+\frac{1}{2}\frac{\delta\Xi}{\delta\rho}\bigg{]}+$
	$\displaystyle+\psi({\bf r})\int\,d^{3}{\bf x}\frac{1}{2}\frac{\delta\Xi}{\delta q({\bf r},{\bf x})}\psi({\bf x})^{2}+$
	$\displaystyle\psi({\bf r})\int\,d^{3}{\bf x}\frac{1}{2}\frac{\delta\Xi}{\delta q({\bf x},{\bf r})}\psi({\bf x})^{2}--\mathcal{E}_{0}\psi({\bf r})$		(B.23)

Together with equations (Eqs. B.21,B.22) the latter equation is transformed into

	$\displaystyle\frac{\ell^{2}}{6}\nabla^{2}\psi({\bf r})=\psi({\bf r})\bigg{[}-\beta pe\varphi({\bf r})-\mathcal{E}_{0}+\frac{1}{2\tilde{v}}\frac{1}{\mu+N\psi({\bf r})^{2}}-$
	$\displaystyle-\frac{\kappa}{(1-\kappa)^{2}}\frac{N\psi({\bf r})^{4}}{2(\mu+N\psi({\bf r})^{2})^{2}}+$
	$\displaystyle\int\,d^{3}{\bf x}\frac{1}{(1-\kappa)^{2}}\frac{N\psi({\bf r})^{2}\psi({\bf x})^{2}}{(\mu+N\psi({\bf r})^{2})(\mu+N\psi({\bf x})^{2}}\psi({\bf x})^{2})\bigg{]}.$		(B.24)