text:
# Neural Propagation of Beliefs

## 1 Introduction

### 1.1 Fundamental Hypothesis

This is the second of two papers elaborating upon the proposition \[Anderson, 1994\] that neural populations encode and process information about analog variables in the form of probability density functions (PDFs). As demonstrated in the first paper \[Barber et al., 2001\], explicit representation of probabilistic descriptors of the state of knowledge of physically relevant variables subserves a powerful strategy for modeling neural circuits. By exploiting mathematical tools developed within the theory of Bayesian inference, we can establish general procedures for building and understanding models of cortical circuits that carry out well-posed information-processing tasks.

Importantly, if neural populations encode the joint probability distribution over the variables of interest, then such neural networks are able to answer any probabilistic question about those variables. Motivated by this fact, we shall now extend the PDF hypothesis by devising methods for embedding joint probabilities into neural networks. This will enable us to construct neural circuit models that pool multiple sources of evidence, such as sensory inputs and any evolutionarily determined priors on the joint distribution. More specifically, we will focus on developing neural networks that use “bottom-up” sensory inputs to build an internal model of the data, which in turn uses “top-down” signals to impose global regularities on the sensory data. In the resulting neural networks, there will naturally arise distinct feedforward, feedback, and lateral connections, analogous to the neural pathways observed in the anatomy of the cerebral cortex \[Van Essen et al., 1992\].

Our treatment is based on Bayesian belief networks \[Pearl, 1988, Smyth et al., 1997\], graphical representations of probabilistic models that provide an efficient means for organizing the relations between the random variables of a given model. The resulting neural networks will have several properties of Bayesian belief networks, as well as more typical neural-network properties, so we will call them *neural belief networks*. Bayesian belief network have been previously utilized in the genesis of certain neural networks (\[Neal, 1992, Zemel, 1999\]), although with different methods for generating the neural network architecture and dynamics.

### 1.2 Three Levels of Representation

In the context of the PDF hypothesis \[Anderson, 1996\], we assert that a physical variable $`x`$ is described by a neural population at time $`t`$ in terms of a PDF $`\rho (x;t)`$, rather than as a single-valued estimate $`x(t)`$. In general, we consider a PDF described at time $`t`$ in terms of a set of $`D`$ parameters $`\left\{A_\mu \right\}`$. Experimentally observed linear decoding rules \[Georgopoulos et al., 1986, Schwartz, 1993\] suggest that the PDFs may be represented as

$$\rho (x;t)=\underset{\mu =1}{\overset{D}{}}A_\mu (t)\mathrm{\Phi }_\mu (x)$$
(1)
The basis functions $`\mathrm{\Phi }_\mu (x)`$ are orthonormal functions which serve to define the PDFs that the neural circuit can represent, but the amplitudes $`A_\mu (t)`$ cannot be interpreted as neuronal firing rates due to their arbitrarily high precision and their ability to take on negative values. Therefore, we introduce an additional representation of the PDF in terms of firing rates $`a_i(t)`$ and decoding functions $`\varphi _i(x)`$ assigned to $`N`$ neurons, so that

$$\rho (x;t)=\underset{i=1}{\overset{N}{}}a_i(t)\varphi _i(x)$$
(2)
Unlike the basis functions $`\mathrm{\Phi }_\mu (x)`$, the decoding functions $`\varphi _i(x)`$ form a highly redundant, overcomplete representation ($`ND`$) suitable for use with neuronal units having biologically realistic precision (some 2–3 bits; see Rieke et al., 1997).

The abstract representation defined by equation 1 will underlie the representation in the neuron space defined by equation 2. This allows us to deal with the issue of how PDFs can be precisely implemented in populations of neurons by focusing on the mapping between the minimal space and the space of neurons \[Barber et al., 2001\]. Thus, neural belief networks can be developed in the theoretically convenient abstract representation, and then be implemented in more realistic networks of low-precision model neurons. Adopting the terminology of Zemel et al. (1998), we denote the set of physical variables as the *implicit space* and the measurable quantities as the *explicit space*. Extending their nomenclature, we denote the abstract space of equation 1 as the *minimal space*. The explicit space of neurons constitutes a biological implementation of the desired computations in the implicit space, while the minimal space, whose properties are more conducive to formal analysis, provides a valuable bridge between the two other spaces.

We have developed rules to transform between representations in the three spaces \[Barber et al., 2001\]. Equations 1 and 2 respectively define how to transform representations in the minimal and explicit spaces into a representation (PDF) in the implicit space. With the remaining rules (summarized below), we can readily switch between the three representations and approach any task in the most appropriate space.

For the orthonormal basis $`\left\{\mathrm{\Phi }_\mu (x)\right\}`$, the minimal space coefficients $`A_\mu (t)`$ are found using the encoding rule

$$A_\mu (t)=\mathrm{\Phi }_\mu (x)\rho (x;t)𝑑x$$
(3)
The coefficients in the explicit space, i.e. the neural firing rates $`a_i(t)`$, cannot be found in this direct fashion. We utilize a set of encoding functions $`\widehat{\varphi }_i(x)`$ and an encoding rule of the form

$$a_i(t)=f\left(\widehat{\varphi }_i(x)\rho (x;t)𝑑x\right)$$
(4)
This explicit-space encoding rule is patterned after the neural responses associated with the population vector of Georgopoulos et al. (1986).

To relate the explicit and minimal spaces, we express the encoding and decoding functions in terms of the orthonormal basis, so that

$`\varphi _i(x)={\displaystyle \underset{\nu =1}{\overset{D}{}}}\kappa _{\nu i}\mathrm{\Phi }_\nu (x)`$ (5)
$`\widehat{\varphi }_i(x)={\displaystyle \underset{\nu =1}{\overset{D}{}}}\widehat{\kappa }_{i\nu }\mathrm{\Phi }_\nu (x)`$ (6)
The transformation coefficients $`\kappa _{\nu i}`$ and $`\widehat{\kappa }_{i\nu }`$ also relate the $`A_\mu (t)`$ and $`a_i(t)`$, such that

$`A_\mu (t)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\kappa _{\nu i}a_i(t)`$ (7)
$`a_i(t)`$ $`=`$ $`f\left({\displaystyle \underset{\nu =1}{\overset{D}{}}}\widehat{\kappa }_{i\nu }A_\nu (t)\right)`$ (8)
Methods for determining the encoding functions $`\widehat{\varphi }_{i(x)}`$, decoding functions $`\varphi _i(x)`$, and transformation coefficients $`\kappa _{\nu i}`$ and $`\widehat{\kappa }_{i\nu }`$ have been addressed in detail in the earlier work \[Barber et al., 2001\].

## 2 Neural Belief Networks

### 2.1 Bayesian Belief Networks

In section 1.2, we have summarized methods for encoding and decoding probability density functions into and from the firing rates of populations of neurons. These methods are oriented towards encoding a single random variable (or vector), but we do not wish to restrict ourselves to only the simplest implicit spaces. In this work, we will explore ways in which we can apply the methods so far developed \[Barber et al., 2001\] to more complicated implicit spaces. In particular, we will use Bayesian belief networks to efficiently organize the implicit random variables, and then use these Bayesian belief networks to generate neural networks.

Bayesian belief networks are directed acyclic graphs that represent probabilistic models (Figure 1).

Each node represents a random variable, and the arcs (or directed line segments) signify the presence of direct causal influences between the linked variables. The strengths of these influences are defined using conditional probabilities. The directionality of a specified link indicates the direction of causality (or, more simply, relevance); an arc points from direct cause to effect.

Bayesian belief networks have two properties that we will find very useful, both of which stem from the independencies shown by the graph structure. First, the value of a node $`X`$ is not dependent upon all of the other graph nodes. Rather, it depends only on a subset of the nodes, called a Markov blanket of $`X`$, that separates node $`X`$ from all the other nodes in the graph. The Markov blanket of interest to us can be readily determined from the graph structure. It is comprised of the union of the direct parents of $`X`$, the direct successors of $`X`$, and all direct parents of the direct successors of $`X`$. Second, the joint probability over the random variables is decomposable as

$$P(x_1,x_2,\mathrm{},x_n)=\underset{i=1}{\overset{n}{}}P(x_i\mathrm{Pa}(x_i))$$
(9)
where $`\mathrm{Pa}(x_i)`$ refers to the (possibly empty) set of direct parent nodes of $`X_i`$. This decomposition comes about from repeated application of Bayes’ rule and from the structure of the graph.

### 2.2 Probabilistic Inference Performed by Neural Networks

Before exploring arbitrary Bayesian belief networks, it is enlightening to consider a BBN with a simple graph consisting of two connected nodes $`XY`$. This graph represents any probabilistic model where a single random variable is inferred from one source of evidence. For convenience, we will work in the minimal space.

Our objective is to find the best marginal PDFs $`\rho (y;t)`$ and $`\rho (x;t)`$ to describe the system. We represent the PDFs using equation 1 and

$$\rho (y;t)=\underset{\nu }{}B_\nu (t)\mathrm{\Psi }_\nu (y)$$
(10)
In this network, we find values of the output firing rates $`\left\{B_\nu (t)\right\}`$ only, with input firing rates $`\left\{A_\mu (t)\right\}`$ assumed to be fully determined by an encoding process.

We define a cost function by

$`E_y`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \left(\rho (y;t)\rho (yx)\rho (x;t)𝑑x\right)^2𝑑y}`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \left(\underset{\nu }{}B_\nu (t)\mathrm{\Psi }_\nu (y)\underset{\mu }{}A_\mu (t)\rho (yx)\mathrm{\Phi }_\mu (x)𝑑x\right)^2𝑑y}`$
The cost function has a minimum corresponding to taking a weighted average of the conditional probability $`\rho (yx)`$

$$\rho (y;t)=\rho (yx)\rho (x;t)𝑑x$$
(12)
We assume in equation 12 that the relationship between $`x`$ and $`y`$ is independent of the values of the underlying parameters $`A_\mu (t)`$ of the minimal space \[Barber et al., 2001\].

To minimize $`E_y`$, we calculate the derivatives $`E_y/B_\nu `$. Using gradient descent, we obtain the update rule

$`{\displaystyle \frac{dB_\nu (t)}{dt}}`$ $`=`$ $`\eta {\displaystyle \frac{E_y}{B_\nu }}`$ (13)
$`=`$ $`\eta \left(B_\nu (t){\displaystyle \underset{\mu }{}}\mathrm{\Omega }_{\nu \mu }A_\mu (t)\right)`$
where $`\eta `$ is a rate constant and we identify weights

$$\mathrm{\Omega }_{\nu \mu }=\mathrm{\Psi }_\nu (y)\rho (yx)\mathrm{\Phi }_\mu (x)𝑑x𝑑y$$
(14)
Equations 13 and 14 define a network in the minimal space which can be converted to a neural network in the explicit space using equation 8.

Equation 13 has a fixed point at

$$B_\nu (t)=\underset{\mu }{}A_\mu (t)\mathrm{\Psi }_\nu (y)\rho (yx)\mathrm{\Phi }_\mu (x)𝑑x𝑑y$$
(15)
This is identical to the result obtained by encoding the inference relation in equation 12 into the minimal space using equation 3.

### 2.3 Predictive and Retrospective Support

Neural belief networks can feature two distinct types of information propagation that provide support for the PDFs represented at each graph node. Predictive support, also called causal support, is probabilistic information that propagates, along the directions of the graph links, from cause to effect \[Pearl, 1988\]. The network considered in section 2.2 involves only predictive support.

The second type of support is retrospective, or diagnostic, support. In this case, information propagates against the directions of the graph links, from effect to cause, or, equivalently, from evidence to hypothesis \[Pearl, 1988\]. By specifying $`\rho (y;t)`$ instead of $`\rho (x;t)`$, the $`XY`$ inference network features only retrospective support.

Using the same minimal-space representation (equations 1 and 10) that we used for the predictive network, we determine the update rule for the retrospective network. We find

$$\frac{dA_\mu (t)}{dt}=\eta \left(\underset{\beta }{}B_\beta (t)\mathrm{\Omega }_{\beta \mu }\underset{\alpha }{}A_\alpha (t)\mathrm{{\rm Y}}_{\alpha \mu }\right)$$
(16)
with feedback weights

$$\mathrm{\Omega }_{\beta \mu }=\mathrm{\Psi }_\beta (y)\rho (yx)\mathrm{\Phi }_\mu (x)𝑑x𝑑y$$
(17)
and lateral weights

$$\mathrm{{\rm Y}}_{\alpha \mu }=\left(\rho (yx)\mathrm{\Phi }_\alpha (x)𝑑x\right)\left(\rho (yx)\mathrm{\Phi }_\mu (x)𝑑x\right)𝑑y$$
(18)
While the feedback weights $`\mathrm{\Omega }_{\beta \mu }`$ of the retrospective network are identical to the feedforward weights of the predictive network (equation 13) in the minimal space, note that the resulting neural networks weights in the explicit space need not be identical. The lateral weights $`\mathrm{{\rm Y}}_{\alpha \mu }`$ provide a measure of the correlation between what the different basis functions in the parent node $`X`$ predict about the child node $`Y`$; these lateral connections act to ensure consistency between evidence and hypothesis.

Although the networks driven by retrospective support are closely related to the networks driven by predictive support, their function is quite different. For example, consider a network with $`\rho (yx)=\delta \left(y\left|x\right|\right)`$ as the underlying computation. In a predictive network, wherein we specify $`\rho (x;t)`$ and infer $`\rho (y;t)`$, the absolute-value relationship is approximated in a straightforward fashion (Figure 2a). Conversely, a retrospective network based upon the same conditional PDF is called upon to “invert” the absolute value, a non-invertible function. The inferred $`\rho (x;t)`$ decoded from the retrospective network (Figure 2b) captures both of the possible solutions—with positive and negative values—in response to a unimodal PDF $`\rho (y;t)`$.

### 2.4 Encoding Bayesian Belief Networks Into Neural Networks

Following a strategy similar to that presented in section 2.2, we can develop neural-network update rules from arbitrary Bayesian belief networks. We assume that the Bayesian belief network consists of $`R`$ nodes, symbolizing random variables $`X_1`$, $`X_2`$, $`X_3`$, …,$`X_R`$. Introducing representations

$$\rho (x_i;t)=\underset{\mu }{}A_\mu ^{(i)}(t)\mathrm{\Phi }_\mu ^{(i)}(x_i)$$
(19)
for the marginal distributions in the minimal space, we define auxiliary cost functions, analogous to that in equation LABEL:eq:xyinferencecostfunction, with the forms

$`E_i`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \left(\rho (x_i;t)\rho (x_i\mathrm{Pa}(x_i))\underset{j}{}\rho (x_j;t)dx_j\right)^2𝑑x_i}`$
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \left(\genfrac{}{}{0pt}{}{{\displaystyle \underset{\mu }{}}A_\mu ^{(i)}(t)\mathrm{\Phi }_\mu ^{(i)}(x_i)}{{\displaystyle \rho (x_i\mathrm{Pa}(x_i))\underset{j}{}\underset{\nu _j}{}A_{\nu _j}^{(j)}(t)\mathrm{\Phi }_{\nu _j}^{(j)}(x_j)dx_j}}\right)^2𝑑x_i}`$
where the index $`j`$ for the product runs over the direct parents of $`X_i`$.

We further define an aggregate cost function

$$E=\underset{i=1}{\overset{R}{}}K_iE_i$$
(21)
The parameters $`K_i`$ can be used to emphasize particular portions of the network; except where otherwise noted, we assume that $`K_i=1`$ for all $`i`$. Employing gradient descent to minimize $`E`$, we find

$$\frac{dA_\sigma ^{(k)}}{dt}=\eta \underset{i=1}{\overset{R}{}}K_i\frac{E_i}{A_\sigma ^{(k)}}$$
(22)
Since $`A_\sigma ^{(k)}`$ does not appear in all of the cost functions, $`E_i/A_\sigma ^{(k)}`$ is nonzero only for $`i=k`$ and when the graph node for $`X_i`$ is one of the direct children of the graph node for $`X_k`$. Thus,

$$\frac{1}{\eta }\frac{dA_\sigma ^{(k)}}{dt}=K_k\frac{E_k}{A_\sigma ^{(k)}}\underset{i\{i:X_i\mathrm{Ch}(X_k)\}}{}K_i\frac{E_i}{A_\sigma ^{(k)}}$$
(23)
As was the case for the two-node inference network, the input PDF or PDFs will be specified by an encoding process rather than an update rule of this sort.

The derivatives in equation 23 are straightforward but lengthy to evaluate. We omit the details of the calculation and simply state the results

$`{\displaystyle \frac{1}{\eta }}{\displaystyle \frac{dA_\sigma ^{(k)}}{dt}}=`$ (24)
$`K_k\left(A_\sigma ^{(k)}{\displaystyle \underset{j}{}}{\displaystyle \underset{\nu _j}{}}A_{\nu _j}^{(j)}\mathrm{\Omega }_{\sigma \nu _1\nu _2\mathrm{}\nu _{j_{max}}}^{(k)}\right)`$
$`+{\displaystyle \underset{i}{}}K_i{\displaystyle \underset{\mu }{}}A_\mu ^{(i)}{\displaystyle \underset{jk}{}}{\displaystyle \underset{\nu _j}{}}A_{\nu _j}^{(j)}\mathrm{\Omega }_{\mu \nu _1\nu _2\mathrm{}\nu _{j_{max}}\sigma }^{(ik)}`$
$`{\displaystyle \underset{i}{}}K_i{\displaystyle \underset{j}{}}{\displaystyle \underset{j^{}k}{}}{\displaystyle \underset{\nu _j}{}}A_{\nu _j}^{(j)}{\displaystyle \underset{\mu _j^{}}{}}A_{\mu _j^{}}^{(j^{})}\mathrm{{\rm Y}}_{\nu _1\nu _2\mathrm{}\nu _{j_{max}}\mu _1\mu _2\mathrm{}\mu _{j_{max}^{}}\sigma }^{(jj^{}k)}`$
where we have defined

$`\mathrm{\Omega }_{\sigma \nu _1\nu _2\mathrm{}\nu _{j_{max}}}^{(k)}=`$ (25)
$`{\displaystyle \mathrm{\Phi }_\sigma ^{(k)}(x_k)\rho (x_k\mathrm{Pa}(x_k))\underset{j}{}\mathrm{\Phi }_{\nu _j}^{(j)}(x_j)dx_jdx_k}`$
$`\mathrm{\Omega }_{\mu \nu _1\nu _2\mathrm{}\nu _{j_{max}}\sigma }^{(ik)}=`$ (26)
$`{\displaystyle \mathrm{\Phi }_\mu ^{(i)}(x_i)\rho (x_i\mathrm{Pa}(x_i))\mathrm{\Phi }_\sigma ^{(k)}(x_k)\underset{jk}{}\mathrm{\Phi }_{\nu _j}^{(j)}(x_j)dx_jdx_idx_k}`$
and

$`\mathrm{{\rm Y}}_{\nu _1\nu _2\mathrm{}\nu _{j_{max}}\mu _1\mu _2\mathrm{}\mu _{j_{max}^{}}\sigma }^{(jj^{}k)}=`$ (27)
$`{\displaystyle }({\displaystyle }\mathrm{\Phi }_\sigma ^{(k)}(x_k)\rho (x_i\mathrm{Pa}(x_i)){\displaystyle \underset{j^{}k}{}}\mathrm{\Phi }_{\mu _j^{}}^{(j^{})}(x_j^{})dx_j^{}dx_k`$
$`{\displaystyle }\rho (x_i\mathrm{Pa}(x_i)){\displaystyle \underset{j}{}}\mathrm{\Phi }_{\nu _j}^{(j)}(x_j)dx_j)dx_i`$
The sums over the index $`i`$ in equation 24 run over the children of $`X_k`$. The products are over the parents of either node $`X_k`$ (equation 25, first term in 24) or node $`X_i`$ (equations 26 and 27, second and third term in 24), possibly excluding node $`X_k`$ itself. From these equations, it can be seen that the PDF represented at a node $`X_i`$ is updated based only on its direct parents, its direct descendents, and the direct parents of its direct descendents, so that $`X_i`$ is separated from all other nodes in the neural belief networks by an appropriate Markov blanket.

Equations 24 through 27 can be awkward to use directly. It can be more convenient to write out the cost functions for a specific probabilistic model using equation LABEL:eq:generalenergies and directly evaluate the necessary derivatives (specified by equation 23).

## 3 Applications of Neural Belief Networks

### 3.1 Bidirectional Propagation

So far, we have restricted our attention to developing neural networks that encode simple probabilistic models involving only a single source of evidence. Given suitable representations, we can design neural networks that capture a wide variety of probabilistic relations \[Barber et al., 2001\]. Further, by regarding representations of probabilistic dependence models as Bayesian belief networks, we have seen that two types of information propagation come into play: predictive and retrospective.

In this section, we will examine several applications of neural belief networks. Unlike those considered before, the probabilistic models will now feature multiple sources of evidence, with bidirectional propagation of both predictive and retrospective support. The corresponding neural networks thus have neurons with both feedforward and feedback connections, as well as lateral connections.

### 3.2 Neural Propagation of Evidence in Trees

To facilitate investigation of information propagation in neural belief networks, we first focus on networks whose implicit spaces are specified by tree-structured Bayesian belief networks (Figure 3a).

These implicit networks are general enough to illustrate the concepts, while yielding neural networks that are readily understood. We will examine binary trees, where each node has at most two children, but the results extend simply to more general tree structures.

We may assume that evidence is only available in the root node and the leaf nodes (although all such nodes need not provide evidence). The root node is the single node which has no father and is located at the top of the tree, while the leaf nodes are all the nodes which have no children. If another node $`X`$ were externally specified, the subtree rooted at $`X`$ could be broken off and treated separately. Conversely, the father node of $`X`$ is unaffected by the descendants of $`X`$, so they could be deleted from the original tree, leaving $`X`$ as a leaf node.

Clearly, an unspecified root node can only receive retrospective support, while an unspecified leaf node can only receive predictive support. All other unspecified nodes in the tree will receive both retrospective and predictive support. We will thus need to consider separately these three types of nodes when determining the update rules for the neural network.

An unspecified root node $`X`$ with children $`Y`$ and $`Z`$ receives feedback inputs (retrospective support) from both of them (Figure 3). We introduce the representations

$`\rho (x;t)`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}A_\alpha (t)\mathrm{\Phi }_\alpha (x)`$ (28)
$`\rho (y;t)`$ $`=`$ $`{\displaystyle \underset{\beta }{}}B_\beta (t)\mathrm{\Psi }_\beta (y)`$ (29)
$`\rho (z;t)`$ $`=`$ $`{\displaystyle \underset{\gamma }{}}C_\gamma (t)\mathrm{\Theta }_\gamma (z)`$ (30)
Following the procedures described in section 2.4, the update rule for the root is

$$\frac{1}{\eta }\frac{dA_\mu }{dt}=K_y\frac{E_y}{A_\mu }K_z\frac{E_z}{A_\mu }$$
(31)
with

$`{\displaystyle \frac{E_y}{A_\mu }}`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}A_\alpha (t){\displaystyle \left(\rho (yx)\mathrm{\Phi }_\alpha (x)𝑑x\right)\left(\rho (yx)\mathrm{\Phi }_\mu (x)𝑑x\right)𝑑y}`$ (32)
$`{\displaystyle \underset{\beta }{}}B_\beta (t){\displaystyle \mathrm{\Psi }_\beta (y)\rho (yx)\mathrm{\Phi }_\mu (x)𝑑x𝑑y}`$
and

$`{\displaystyle \frac{E_z}{A_\mu }}`$ $`=`$ $`{\displaystyle \underset{\alpha }{}}A_\alpha (t){\displaystyle \left(\rho (zx)\mathrm{\Phi }_\alpha (x)𝑑x\right)\left(\rho (zx)\mathrm{\Phi }_\mu (x)𝑑x\right)𝑑z}`$ (33)
$`{\displaystyle \underset{\gamma }{}}C_\gamma (t){\displaystyle \mathrm{\Theta }_\gamma (z)\rho (zx)\mathrm{\Phi }_\mu (x)𝑑x𝑑z}`$
Thus, the firing rates for the root node are driven by a sum of feedback inputs that individually are identical to the input produced by a single source of retrospective support (section 2.3). The parameters $`K_y`$ and $`K_z`$ need not be the same; different values may be used to give greater significance to one of the inputs.

The firing rates for an unspecified leaf node are also updated by a familiar rule. Consider a leaf node $`X`$ with father $`U`$. Using equation 28 and

$$\rho (u;t)=\underset{\delta }{}D_\delta (t)\mathrm{\Lambda }_\delta (u)$$
(34)
we obtain the update rule

$$\frac{1}{\eta }\frac{dA_\mu }{dt}=K_x\frac{E_x}{A_\mu }$$
(35)
where

$$\frac{E_x}{A_\mu }=A_\mu (t)\underset{\delta }{}D_\delta (t)\mathrm{\Lambda }_\delta (u)\rho (xu)\mathrm{\Phi }_\mu (x)𝑑u𝑑x$$
(36)
This update rule is of course identical to the update rule for the $`XY`$ network with predictive support (section 2.2).

All nonroot, nonleaf nodes have similar update rules. For a node $`X`$ with father $`U`$ and sons $`Y`$ and $`Z`$, we impose the representations given above in equations 28 through 30 and equation 34. Applying the procedure of section 2.4 yields an update rule with form

$$\frac{1}{\eta }\frac{dA_\mu }{dt}=K_y\frac{E_y}{A_\mu }K_z\frac{E_z}{A_\mu }K_x\frac{E_x}{A_\mu }$$
(37)
The partial derivatives of the cost functions are identical to those evaluated for the root and leaf nodes. The descendants and ancestors of $`X`$ thus communicate, in the neural network, only through the intermediary of $`X`$ itself. This is consistent with the tree structure of the underlying Bayesian belief network; given $`X`$, the descendants of $`X`$ are conditionally independent of the ancestors of $`X`$ \[Pearl, 1988\].

With these update rules, evidence provided at the root node or at leaf nodes will propagate throughout the network. The manner in which evidence is specified will depend on how the probabilistic model is posed. Therefore, the same graph could have different nodes specified for different purposes. For instance, the small tree in Figure 3b could have any of its nodes represent sensory inputs.

If we specify the PDF for the root node, the leaf nodes $`Y`$ and $`Z`$ will receive predictive support from $`X`$. By selecting appropriate representations for the PDFs and for the conditional probabilities associated with the links, the resulting neural network could subdivide a complex, highly general sensory input into simpler, more specialized components. For example, a visual input at a particular retinal location might be separated into contrast and color.

Conversely, if we specify PDFs for the leaf nodes, the root node $`X`$ will receive retrospective support from $`Y`$ and $`Z`$. Amongst other possibilities, this provides a simple way to model redundant sensory inputs. If we take the conditional probabilities to be of narrow Gaussian form, $`\rho (yx)=N(y;x,\sigma _y^2)`$ and $`\rho (zx)=N(z;x,\sigma _z^2)`$, then the firing rates representing $`\rho (x;t)`$ will be updated so as to pool the diagnostic information from both the sensory inputs, $`Y`$ and $`Z`$.

Similar arguments apply to larger trees. Additionally, larger trees may well have the root node and leaf nodes specified simultaneously (which is not of interest for the small tree discussed above). These nodes may correspond to sensory inputs, or can represent priors that are built into the neural network.

It is important to recognize that the neural update rules developed above apply only to the binary trees. If the arrows in the binary tree graphs are reversed, rather different neural update rules are produced. For example, the directions of the small tree shown in Figure 3b can be reversed, as shown in Figure 3c. Specifying $`Y`$ and $`Z`$ yields the multiplicative update rule

$`{\displaystyle \frac{1}{\eta K_x}}{\displaystyle \frac{dA_\alpha (t)}{dt}}=A_\alpha (t)`$
$`+{\displaystyle \underset{\beta ,\gamma }{}}B_\beta (t)C_\gamma (t){\displaystyle \mathrm{\Phi }_\alpha (x)\rho (xy,z)\mathrm{\Psi }_\beta (y)\mathrm{\Theta }_\gamma (z)𝑑x𝑑y𝑑z}`$
while specifying $`X`$ and $`Z`$ yields

$`{\displaystyle \frac{1}{\eta K_x}}{\displaystyle \frac{dB_\nu (t)}{dt}}=`$
$`{\displaystyle \underset{\alpha ,\gamma }{}}A_\alpha (t)C_\gamma (t){\displaystyle \mathrm{\Phi }_\alpha (x)\rho (xy,z)\mathrm{\Psi }_\nu (y)\mathrm{\Theta }_\gamma (z)𝑑x𝑑y𝑑z}`$
$`{\displaystyle \underset{\beta ,\gamma _1,\gamma _2}{}}B_\beta (t)C_{\gamma _1}(t)C_{\gamma _2}(t){\displaystyle }\{\left({\displaystyle }{\displaystyle }\rho (xy,z)\mathrm{\Psi }_\beta (y)\mathrm{\Theta }_{\gamma _1}(z)dydz\right)`$
$`\left({\displaystyle }{\displaystyle }\rho (xy,z)\mathrm{\Psi }_\nu (y)\mathrm{\Theta }_{\gamma _2}(z)dydz\right)\}dx`$
This latter update rule features a feedforward term that is multilinear in the firing rates $`\left\{A_\alpha (t)\right\}`$ and $`\left\{C_\gamma (t)\right\}`$. It also features a nonlinear lateral combination of the firing rates $`\left\{B_\beta (t)\right\}`$ and $`\left\{C_\gamma (t)\right\}`$ which serves to ensure that the two parent nodes $`Y`$ and $`Z`$ are mutually consistent with the PDF of the child node $`X`$. If only $`X`$ is specified, there will be an additional update rule for the firing rates $`\left\{C_\gamma (t)\right\}`$ that is similar in form to equation LABEL:eq:reversetreeretro.

Although the update rules are more complicated with the directions reversed in this manner, the probabilistic model may demand it. For instance, the Bayesian belief network in Figure 3c is appropriate for implementing the arithmetic operations (add, subtract, multiply, and divide).

An interesting application of these two types of neural belief networks (tree and reversed-tree) is the estimation of the velocity of a moving object. A small tree (Figure 3b) can generate two copies $`Y`$ and $`Z`$ of the input position $`X`$ by taking the conditional probabilities to be Dirac delta functions $`\delta (yx)`$ and $`\delta (zx)`$. We can set the parameters $`K_i`$ so that the values of the copies will be held for different lengths of time. In particular, we can establish the relations

$`\rho (y;t)`$ $`=`$ $`\rho (x;t\tau )`$ (40)
$`\rho (z;t)`$ $`=`$ $`\rho (x;t2\tau )`$ (41)
These copies of the position can then be used as the inputs to a Bayesian belief network of the type shown in Figure 3c. By setting $`\rho (xy,z)=\delta (x(yz)/\tau )`$, the velocity (given here by $`x`$) at time $`t\tau `$ can be estimated (Figure 4).

### 3.3 Top-Down Feedback From a High-Level Model

In section 3.2, we examined the neural update rules for Bayesian belief networks having binary-tree structure. By doing so, we examined the means by which information propagates throughout the network from one or more sources. In particular, we saw that conditional independence in the Bayesian belief network was preserved in the neural network connectivity.

We now turn from the general rules by which probabilistic information propagates in the neural network and investigate in more detail the effects that multiple sources of evidence have on the encoded PDFs. We consider PDFs encoded by chains of nodes with evidence provided at one or both ends. (Specifying the PDF for any other node breaks the chain into two chains that can be treated independently.) A chain is a special case of a tree, so the neural update rules can be obtained as in section 3.2.

For the chosen demonstration, we employ representations of higher dimensionality than those generally used in the preceding developments. For each graph node, the interval $`[1,1]`$ is covered with ten Gaussians (Figure 5a), from which an orthonormal basis is generated. This makes it possible to represent a variety of PDFs, including multimodal distributions.

The neural firing rate patterns in the explicit space are also assumed to be of Gaussian form. To eliminate any possibility of negative firing rates, the nonlinear activation function is taken to be rectification.

We have already studied at some depth the behavior of chains consisting of two nodes (sections 2.2 and 2.3). These chains are able to transmit probabilistic information from one node of the graph to another, with the accuracy limited by the representations adopted. For the predictive network, a feedforward input $`X`$ drives the output $`Y`$.

To keep the focus on the interaction of multiple sources of evidence, we take the conditional probability to be a Dirac delta function, $`\rho (yx)=\delta (yx)`$. Of course, this defines a version of the communication channel that we have extensively considered, but the goal is now to duplicate the PDF rather than a particular numerical value. As seen in Figure 5b, the input PDF $`\rho (x;t)`$ is copied with great fidelity by the output PDF $`\rho (y;t)`$.

Every chain with three or more nodes will admit both predictive and retrospective support. The behavior of all such chains is well characterized by a chain with three nodes: the update rule for any node depends only on the neighboring nodes, so a chain with three nodes covers all possibilities (only predictive support, only retrospective support, and both predictive and retrospective support). Therefore, we consider the effect of adding a third node to the chain (Figure 3d). We set $`\rho (zy)=\delta (zy)`$ along with $`\rho (yx)=\delta (yx)`$ and use identical parameters $`K_y`$ and $`K_z`$.

The first possibility is just to add the third node without adding any additional evidence to the network. We decode $`\rho (y;t)`$ and $`\rho (z;t)`$ from neural firing rates determined using update rules derived previously for more general trees. For the same input PDF as above, $`\rho (y;t)`$ is unchanged by the addition of the retrospective support from the third node, and the structure of $`\rho (z;t)`$ is identical to that of $`\rho (y;t)`$.

Since the same evidence was presented to the network, the introduction of a third node did not change the behavior of the original nodes. This is entirely appropriate, given the probabilistic foundations of the neural belief networks. However, it is feasible that, by adjusting the $`K_i`$ parameters or the PDFs representable by the minimal spaces, there can be a change in the dynamics of neural networks extended in this fashion. For example, appropriate parameter choice could produce a slowly varying PDF in $`Z`$ that is relatively stable against noise, thus stabilizing a more rapidly varying PDF in $`Y`$.

A more proactive role that the third node $`Z`$ can play is as an additional source of evidence. We directly specify $`\rho (z;t)`$ and encode it into the neural network. The neural firing rates for $`Y`$ are driven by predictive support from $`X`$ and retrospective support from $`Z`$, and then decoded to find $`\rho (y;t)`$. One possible use of this second source of evidence is to resolve an ambiguous input; the inferior mode of a bimodal predictive input can be de-emphasized using more specific retrospective evidence (Figure 5c).

Although this use of the retrospective evidence resolves the ambiguous predictive input, it does not explain how the retrospective evidence comes about. Ideally, we would like to build up a high-level model of the predictive input, and use this model to generate a top-down signal that imposes global regularity on the bottom-up predictive input. The uniform priors that we have assumed throughout this work are somewhat of a hindrance for this purpose. With these priors, probabilities account well for multiple possibilities, but here we want to choose only one of those possibilities.

In principle, we could restrict the allowed PDFs for $`\rho (z;t)`$ by imposing a prior on the PDFs and rederiving our neural update rules. However, we will adopt a more direct approach that illustrates the first steps towards implementing decision theory \[Raiffa, 1968\] in neural networks and take the neurons representing $`\rho (z;t)`$ to be winner-take-all units. Only one of these neurons will be active at a time, based on a simple utility function: the neurons compete to be active, and the single neuron that is most strongly driven will be the one that is activated. The manner of implementation of the winner-take-all units is immaterial, so we directly choose the most strongly driven neuron in the computer simulations. (It is possible to implement a set of winner-take-all units in a real network through a suitable choice of nonlinear activation function and lateral weights. One may let each neuron inhibit the others and have a self-excitatory connection. See Hertz et al., 1991.)

Since only the most strongly driven winner-take-all neuron is activated, this strategy provides us with a way to generate a high-level model that selects the dominant mode of a multimodal input distribution. We allow $`Z`$ to be driven by the predictive support from $`Y`$, and the winner-take-all nature of the $`Z`$ neurons permits only narrow Gaussian PDFs to be represented. Thus, $`\rho (z;t)`$ serves as our high-level model, and can resolve ambiguous inputs, as demonstrated in Figure 5d. This type of neural network could be used as the starting point for coherent theoretical accounts of attentive effects in the primate visual system, of the electrosensory system of weakly electric fish, and of other neural systems where an internal model is built up to impose global constraints on neural representations of information.

## 4 Conclusions

We have extended the hypothesis that neural networks represent information as probability density functions. These PDFs are assumed to be obtainable by a linear combination of some implicit decoding functions, with the decoder for each neuron being weighted by its firing rate. The firing rates in turn are obtainable from the PDF using a complementary set of encoding functions.

Success in representing an individual PDF with a population of neurons \[Barber et al., 2001\] has led us to inquire whether more complex probabilistic models can be represented and implemented in neural terms. To this end, we have adopted graphical representations of probabilistic dependence models, in the form of Pearl’s Bayesian belief networks, to organize and simplify the relations between the random variables in such a model.

In summary, we have introduced, analyzed, and applied three ways to represent probabilistic information. These are the implicit model, depicted as a Bayesian belief network; the representation of the probabilistic model in the minimal space; and the representation of the probabilistic model as a neural network in the explicit space. Further, we have devised rules that permit us to convert one type of representation into another type. This formalization of the representation and processing of information in neurobiological computation therefore suggests a general protocol by which probabilistic models can be embedded in neural networks. First, we specify the probabilistic model, using a Bayesian belief network to organize the random variables. We then consider what sorts of functions are exemplars of the PDFs describing the random variables and use these exemplars to define the minimal space. Finally, we utilize the relations between the minimal space and the explicit space to generate the neural network itself.

## Acknowledgements

This research was supported by the National Science Foundation under Grant Numbers IBN-9634314, PHY-9602127, and PHY-9900713. While writing this work, MJB was a member of the Graduiertenkolleg Azentrische Kristalle, GK549 of the DFG; the writing was completed while MJB and JWC were participants in the Research Year on “The Sciences of Complexity: From Mathematics to Technology to a Sustainable World” at the Zentrum für interdisziplinäre Forschung, University of Bielefeld.


format: text

text:
# 1 The viable parameter space of the COBE-normalized hybrid inflation model. The red (dashed) line is the lower bound on 𝑘 from the SUSY-breaking soft terms. The thick black solid line is the upper bound on 𝑘 from the requirement 𝜑_𝑁_𝑒2.0×10⁻⁶. We also show the contours for fixed Λ_𝐼. The lines except the red (dashed) line are cited (with some modification) from Ref.[]. We thank the authors for permission.

## Acknowledgments

We thank Tsutomu Yanagida for his careful reading of this manuscript, precious comments, and continuous encouragements. T.W. thanks the Japan Society for the Promotion of Science for financial support.


text:
# Weak violation of universality for Polyelectrolyte Chains: Variational Theory and Simulations

## I Introduction

Charged polymers play a central role in theory and experiments. Polyelectrolytes are macromolecules containing ionizable groups. They dissociate in water to form charged groups and low-molecular-weight counter ions. It is important to remark that interactions between different polyions in a $`manychainsystem`$ are of crucial importance for any realistic system with a low but finite overlap concentration. The problem of the $`singlechain`$ conformation is the one of the essential questions for polyelectrolyte solutions and has been under investigation in the recent years. This will be the topic of the present paper as well, but in a more general context. As common to most of the recent investigation, we consider the limit of infinite dilution and zero salt concentration for polyelectrolytes. In this limit a polyelectrolyte can be represented as a connected sequence of charged and uncharged monomers in a “dielectric vacuum” that constitutes the solvent. ¿From a scaling point of view, the properties of single polyelectrolyte chains, despite the long time debates are still not completely clear. Indeed, for the neutral polymer chain system, where a clear understanding of the scaling properties has been achieved , the details of the local chain structure are not essential, while in the polyelectrolyte chain problem they could be important . In particular the long range nature of the Coulomb interaction couples via the presence of the solvent, the short and long length scales of the system . Thus a scaling approach to polyelectrolytes results to be more involved than for neutral chains. The interaction between monomers, counterions and solvent molecules leads to the important phenomenon of counter ion condensation  and a realistic description should take into account all the different characters of this scenario. On the other hand, a minimal or bare essential description of the system is required in order to achieve a quantitative description of the single polyelectrolyte chain problem.

In this paper we will present a generalized form of des Cloizeaux Gaussian variational method  to study the scaling behavior and local properties of polyelectrolytes. Variational techniques have been shown to be well suited to describe a self-interacting polymer with long range repulsion . Scaling properties of the end-to-end distance of the chain have been discussed for the long - ranged intra-molecular interaction $`1/r^\lambda `$ at an arbitrary spatial dimensionality $`d`$. A recent numerical implementation of such a variational technique, for the case of Coulombic interaction (when $`\lambda =d2`$), seems to be promising as well . On the other hand, the validity of the variational approach for short range interacting homopolymer chains must be questioned  for natural and well known reasons. The most studied example in the context of polymer physics is the problem of the self avoiding chain. It was proven that the problem shows universality in the limit of infinitely long chains  and can be treated by field theoretic methods common to phase transitions. Indeed the renormalization group (RG) theory manifests these statements and develops systematic methods to compute the scaling exponents, e.g. for the size of the chain . The variational technique provides an alternative approach since most of the approximations can be controlled. Although, it is well known from elementary statistical mechanics that the variational technique cannot provide such powerful tool due to their mean - field character. Therefore this technique is expected to work in the mean field limit, i.e. at large dimensions (at the upper critical dimension) or for long range potentials. Thus it is well investigated and generally believed that the Gaussian variational approach is better suited for long-range interacting chains.

Moreover, it has been suggested that in the study of random heteropolymer  and directed polymer in a random media , a variational approach in replica space can be also successfully applied. We postpone these issues for future investigations. Here we will mainly concentrate on the pure long-range and short-range homopolymer problems within the Gaussian variational technique.

In section II we introduce the Gaussian polymer model and calculate the variational free energy using a general Gaussian trial Hamiltonian in the spirit of des Cloizeaux  and reference . We then critically review the asymptotic analysis given originally in reference  and suggest a new one which leads to the Flory exponents corrected by the presence of logarithmic factors. In section III we describe a numerical technique in order to solve the Euler equations corresponding to the variational principle, by means of a new algorithm specially devised for this problem. An independent algorithm, that simply minimizes the variational free energy is shown to be quantitatively equivalent to the direct solution of the Euler equations. The numerical solution of the Euler equations is then carefully analyzed for polymer chains of increasing size. We show here that our new asymptotic analysis is in quantitative agreement with our numerical findings. We discuss connections with the short range homopolymer problem within the Gaussian Variational Principle and results for the end-to-end distance are obtained. In section IV Monte Carlo simulation results, obtained by means of a pivot algorithm, are presented, showing quantitative agreement with our variational analysis. The case of a self-interacting homopolymer chain, when an attractive two-body interaction between the monomers is considered together with a three body repulsive term, can also be studied within our variational approach. Preliminary calculations shows that the globular formation, expected from simple phenomenological arguments and due to the balancing between the two body and three body interactions, takes place. The instability of the globule state due to a Coulombic repulsion between the monomers is also described by our variational principle and will be presented elsewhere . Conclusions and further prospects are summarized in section V.

## II The Gaussian Variational principle

### A Basic equations

First we start from a well known example for a Gaussian variational approach suitable for the investigations presented below. We introduce a simple problem to compute the free energy of a self-interacting polymer chain with Coulombic repulsion between the monomers, based on a discrete representation of a chain of $`N`$ monomers. The method relies on the well known variational principle, where a Gaussian trial probability is proportional to the exponential of a quadratic form of the monomer coordinates $`\stackrel{}{r_j}`$,

$$P_V(\stackrel{}{r_1},..,\stackrel{}{r_N})=Z_V^1\mathrm{exp}\{H_0(\stackrel{}{r_1},..,\stackrel{}{r_N})\}.$$
(1)
In eq.(1) the trial Hamiltonian $`H_0`$ is given by

$$H_0(\stackrel{}{r_1},..,\stackrel{}{r_N})=\frac{d}{2}\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{N}{}}G_N(ij)(\stackrel{}{r_i}\stackrel{}{r_j})^2$$
(2)
and $`Z_V`$ is a normalization constant determined by the condition

$$P_V(\stackrel{}{r_1},..,\stackrel{}{r_N})d^d\stackrel{}{r_1}\mathrm{}d^d\stackrel{}{r_N}=1.$$
(3)
Up to this point different forms of the Gaussian variational principle could be considered. We decided here to consider the variational form originally proposed by des Cloizeaux  in the context of the self-avoiding homopolymer problem. The main essential features of the Gaussian variational principle can be obtained within this choice of the variational form according to equation (2). More sophisticated forms of the variational kernel including a classical path $`\stackrel{}{r_0}`$ and eventual anisotropy along parallel and perpendicular directions to the polymer backbone, can be also included . By the appropriate choice of a periodic boundary conditions for the chain the quadratic form $`H_0`$ can be diagonalized by introducing cyclic coordinates (Rouse modes) for the monomer positions

$`\stackrel{}{\rho _q}`$ $`=`$ $`N^{\frac{1}{2}}{\displaystyle \underset{j=1}{\overset{N}{}}}\mathrm{exp}[\text{i}2\pi jq/N]\stackrel{}{r_j}`$ (4)
$`\stackrel{}{r_j}`$ $`=`$ $`N^{\frac{1}{2}}{\displaystyle \underset{q=1}{\overset{N}{}}}\mathrm{exp}[\text{i}2\pi jq/N]\stackrel{}{\rho _q}.`$ (5)
Indeed, due to cyclic invariance the Cartesian components of $`\stackrel{}{\rho }_q`$ satisfy

$$\rho _q^{(j)}\rho _q^{}^{(j^{})}=\frac{1}{d}\delta _{jj^{}}\delta _{qq^{}}g_N^1(2\pi q/N),$$
(6)
where $`g_N^1(2\pi q/N)`$ is positive and discontinuous function of $`k=2\pi q/N`$ and is related to the Gaussian propagator $`G`$ of equation (2) via the following expression

$$g_N(2\pi q/N)=\underset{n=1}{\overset{N1}{}}G_N(n)[1\mathrm{cos}(2\pi qn/N)].$$
(7)
The variational free energy for a uniformly charged polymer of polymerization index $`N`$ in a $`d`$ dimensional space can be computed according to the usual Feynman inequality

$$FF_V=HH_0_0+F_0,$$
(8)
where $`F_0`$ is the free energy of the Gaussian model defined by equation (2). In equation (8) $`H`$ represents the full Hamiltonian of a charged polymer chain with a pure long - range form for the monomer-to-monomer interaction,

$$V(|\stackrel{}{r_i}\stackrel{}{r_j}|)=|\stackrel{}{r_i}\stackrel{}{r_j}|^\lambda ,$$
(9)
where $`\lambda <d`$. The expression for the Hamiltonian reads

$$H=\frac{d}{2b^2}\underset{i=1}{\overset{N}{}}(\stackrel{}{r_{i+1}}\stackrel{}{r_i})^2+\beta /2\underset{i=1}{\overset{N}{}}\underset{j=1}{\overset{N}{}}V(|\stackrel{}{r_i}\stackrel{}{r_j}|),$$
(10)
where $`\beta =l_\mathrm{B}^\lambda /A^2`$ depends on the Bjerrum length $`l_\mathrm{B}=(e^2/4\pi ϵk_BT)`$ and where the distance between charged monomers along the chain is given by $`A`$ in units of the Kuhn segment $`b`$. $`A`$ defines thus the charge fraction $`f=1/A`$.

It is straightforward to derive the following expression for the variational free energy $`F_V`$ in the limit of infinitely long chain $`N\mathrm{}`$,

$`F_V[g(k)]/N={\displaystyle \frac{d}{2}}b(1)`$ $`+`$ $`\beta \left({\displaystyle \frac{d}{2}}\right)^{\frac{\lambda }{2}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{d\lambda }{2}\right)}{\mathrm{\Gamma }\left(\frac{d}{2}\right)}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}[b(n)]^{\frac{\lambda }{2}}+`$ (11)
$`{\displaystyle \frac{d}{4\pi }}{\displaystyle _\pi ^{+\pi }}\mathrm{log}g(k)𝑑k`$ $``$ $`{\displaystyle \frac{d}{2}}[1+\mathrm{log}(2\pi /d)],`$ (12)
where

$$b_N(ij)=\frac{1}{N}_\pi ^{+\pi }\frac{(1\mathrm{cos}((ij)k))}{g(k)}𝑑k,$$
(13)
$$b_N(ij)(\stackrel{}{r_i}\stackrel{}{r_j})^2,$$
(14)
and

$$\underset{N\mathrm{}}{lim}g_N(k_N)g(k),g(k)0.$$
(15)
The function $`g(k)`$ satisfies the symmetry requirements

$$g(k)=g(k),g(k+2\pi )=g(k).$$
(16)
Minimization of the variational free energy $`F_V`$ yields

$`g(k)=1\mathrm{cos}(k)`$ $``$ $`\beta \left({\displaystyle \frac{\lambda }{2}}\right)\left({\displaystyle \frac{d}{2}}\right)^{\frac{\lambda }{2}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{d\lambda }{2}\right)}{\mathrm{\Gamma }\left(\frac{d}{2}\right)}}`$ (17)
$`\times {\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1`$ $``$ $`\mathrm{cos}(nk))[b(n)]^{\frac{\lambda +2}{2}}.`$ (18)
Equations (12)-(17) represent the generalization of the Gaussian variational principle for the long-range potential (9) which was introduced first in reference . For the pure Coulombic case $`\lambda =d2`$, comparing the expression for the variational free energy (12) with the correspondent expression for the short range interaction problem ,

$`F_V[g(k)]/N={\displaystyle \frac{d}{2}}b(1)`$ $`+`$ $`w{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}[b(n)]^{\frac{d}{2}}+`$ (19)
$`{\displaystyle \frac{d}{4\pi }}{\displaystyle _\pi ^{+\pi }}\mathrm{log}g(k)𝑑k`$ $``$ $`{\displaystyle \frac{d}{2}}[1+\mathrm{log}(2\pi /d)],`$ (20)
leads to the following conclusion which holds within the Gaussian Variational Principle. All universal exponents for the Coulombic case are related to the corresponding exponents for the short-range interaction case through the dimensional shift

$$dd2.$$
(21)
The result of minimization of equation (20) reads

$$g(k)=1\mathrm{cos}(k)w\underset{n=1}{\overset{\mathrm{}}{}}(1\mathrm{cos}(nk))[b(n)]^{\frac{d+2}{2}},$$
(22)
which is of the same form like equation (17).

Before moving to the next section, where a careful numerical implementation of the Euler equations (13)-(17) for polymer chains of increasing length $`N`$ is presented, we revise the asymptotic analysis of the Euler equations in the limit of small $`k0`$,which was originally given by des Cloizeaux . The behavior of $`g(k)`$ at small $`k`$ determines the asymptotic behavior of $`b(n)`$ for $`n>>1`$ via equation (13). Let us assume first that in this limit $`g(k)`$ has the following form

$$g(k)g|k|^{1+2\nu },\frac{1}{2}<\nu <1,$$
(23)
then the correlation function $`b(n)`$, reads

$$b(n)bn^{2\nu }.$$
(24)
We will now proceed to analyze the Euler equations in the scaling limit of small $`k`$ values within this simple ansatz.

### B The asymptotic analysis - neutral chains

To demonstrate several problems within the Gaussian variational principle we are going to rewise in this section the asymptotic analysis of the Euler equation (22) for neutral chains, which was originally given in the context of the self avoiding chain problem , and later extended to self interacting polymer chains with long-ranged Coulomb like interactions . It has been suggested that due to the long ranged nature of the interactions between the monomers the Gaussian Variational Principle is particularly well suited to describe the asymptotic behavior of polyelectrolyte chains. Special attention has also been devoted  to the case of screened Coulombic interaction, within the Debye Hückel approximation.

Let us consider the interaction term

$$I(k)=w\underset{n=1}{\overset{\mathrm{}}{}}(1\mathrm{cos}(nk))[b(n)]^{\frac{d+2}{2}}.$$
(25)
We can obtain the asymptotic expression for $`I(k)`$ in the $`k0`$ limit, substituting the series with an integral  and assuming that the dominant contribution would correspond to large values of $`n`$. The result can be written in the form

$$I(k)w\frac{\pi b^{\frac{d+2}{2}}}{2\mathrm{\Gamma }(\delta )\mathrm{cos}(\frac{\pi }{2}\delta )}k^{\delta 1}+wb^{\frac{d+2}{2}}\zeta (\delta 2)\frac{k^2}{2},$$
(26)
where $`\delta =\nu (d+2)`$.

An alternative derivation can be obtained in the asymptotic regime $`n>>1`$ evaluating exactly the series

$$S(\delta ,z)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{z^n}{n^\delta },$$
(27)
at $`z=\mathrm{exp}(ik)`$ . The explicite calculations are given in Appendix where it is stressed that we should be distinguish between noninteger and integer $`\delta `$. We apply these results for calculation of eq. (25) where $`b(n)`$ is given by eq.(24).In this case we find

$$I(k)wb^{\frac{d+2}{2}}\left[S(\delta ,z=1)\text{R}eS(\delta ,z)\right],$$
(28)
where $`\delta =\nu (d+2)`$. Considering non-integer values of the exponent $`\delta `$ (this assumption will be checked afterwards) eq. (56) yields

$`I(k)`$ $`=`$ $`wb^{\frac{d+2}{2}}\left[{\displaystyle \frac{\pi }{2\mathrm{\Gamma }(\delta )\mathrm{cos}(\frac{\pi }{2}\delta )}}k^{\delta 1}{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}\zeta (\delta 2p)(1)^p{\displaystyle \frac{(k)^{2p}}{(2p)!}}\right].`$ (29)
This equation generalizes eq. (26) ( see also eq. (B.9) of ref.  ). The Euler equation (22) together with the expansion (29) gives, in the $`k0`$ limit

$`gk^{1+2\nu }={\displaystyle \frac{1}{2}}k^2wb^{\frac{d+2}{2}}\zeta (\delta 2){\displaystyle \frac{k^2}{2}}+w{\displaystyle \frac{\pi b^{\frac{d+2}{2}}}{2\mathrm{\Gamma }(\delta )\mathrm{cos}(\frac{\pi }{2}\delta )}}k^{\delta 1}+O(k^4),`$ (30)
where $`\delta =\nu (d+2)`$.

Let us impose first the balancing conditions in the spirit of des Cloizeaux: (i) the first two terms on the r.h.s. balance each other, (ii) the term on the l.h.s. and the third term on the r.h.s. are also balanced. Then we immediately obtain

$$\nu =\frac{2}{d},$$
(31)
as a valid exponent, but together with the additional conditions

$`w`$ $`b^{\frac{d+2}{2}}\zeta (\delta 2)=1,`$ (32)
$`g`$ $`={\displaystyle \frac{w\pi b^{\frac{d+2}{2}}}{2\mathrm{\Gamma }(\delta )\mathrm{cos}(\frac{\pi }{2}\delta )}},`$ (33)
contained in the above equations, the physical interpretation becomes, however, more complicated. Moreover, we must add to these equations the general relation between $`g`$ and $`b`$, which originates from equation (13)

$$b=\frac{1}{2g\mathrm{\Gamma }(1+2\nu )\mathrm{sin}(\pi \nu )}.$$
(34)
Equation (31) is thus the well known result for the exponent $`\nu `$, whereas equations (32)-(34) set up three equations for the two unknowns $`b`$,$`g`$. Equations (32)-(34) determines not only the amplitudes $`g`$ and $`b`$ but also fixes the value of the interaction parameter $`w`$ itself. This delicate point, (see equation (IV.52) of reference ) has already been discussed by Allegra . This is contradictory in the sense that generally the interaction parameter $`w`$ is an arbitrary quantity. The proper balance condition must consist of a single equation, whereas the conditions (32)-(33), derived in the asymptotic analysis of des Cloizeaux, overconstrains the parameters involved. Besides that the parameter $`\delta =2(d+2)/d`$ is only non-integer for $`d=3`$ and is integer for $`d=1,2,4`$ (where the Flory scaling produces exact results). This indicates that the des Cloizeaux solution can not be valid in this later case since it relies on the expansion (29) which (see Appendix) is only valid for non-integer values of $`\delta `$.

In order to avoid the overconstraints we try to impose the Flory balancing conditions : (i) the elastic term, which is proportional to $`k^2`$, is balanced by the interaction term eq.(25), (ii) the term on the l.h.s. (or the entropic term) is negligible compare to the elastic and interaction terms on the r.h.s. Obviously the Flory balancing conditions have one constraint less. In this case for $`S(\delta ,z)`$ we must use eq.(61) where $`\delta =m=\nu (d+2)=3`$ and $`z=\mathrm{exp}(ik)`$. This immediately leads to the well known Flory exponent $`\nu _\mathrm{F}=3/(d+2)`$ . Then the interaction term (25) becomes

$`I(k)=\underset{k0}{lim}wb^{\frac{d+2}{2}}\left[S(m=3,z=1)ReS(m=3,z)\right]=`$ (35)
$`wb^{\frac{d+2}{2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1\mathrm{cos}(nk)}{n^3}}wb^{\frac{d+2}{2}}\left[\psi (3)\psi (1)\mathrm{log}k\right]{\displaystyle \frac{k^2}{2}}+𝒪(k^4).`$ (36)
The resulting Euler equation takes the form

$`g`$ $`k^{1+2\nu _\mathrm{F}}={\displaystyle \frac{1}{2}}k^2wb^{\frac{d+2}{2}}\left[\psi (3)\psi (1)\mathrm{log}k\right]{\displaystyle \frac{k^2}{2}}+O(k^4).`$ (37)
Eq.(37) can be also obtained directly from eq.(30) if we assume $`\delta =3+ϵ`$, where $`ϵ0`$. Then the poles which have its origin from the second and third terms on the r. h. s. are canceled (see e.g. Appendix and sec. 1.11 in ref.) but the term $`k^2\mathrm{log}k`$ appears. The presence of the $`k^2\mathrm{log}k`$ term in equation (37) leads to the conclusion that the pure power law ansatz $`g(k)k^{1+2\nu }`$ (or $`b(n)n^{2\nu }`$) does not solve in the asymptotic limit $`k0`$ the Euler equation (22) once the exponent $`\delta =(d+2)\nu `$ is assumed to be integer. In fact, the direct numerical solution of the Euler equations (17) and (22) presented in the next section does not show a pure power law behavior in the asymptotic regime $`k<<1`$. We also see that the term $`k^2\mathrm{log}k`$ naturally appears in the previous expressions for the series $`S(\delta ,k)`$, which suggests to consider logarithmic corrections in the initial ansatz for $`b(n)`$. This will allow to stay within a Flory-like solution (corresponding to the integer value $`\delta =\nu (d+2)=3`$) and also to compute explicitly the exponent that characterize the logarithmic correction to the pure power law scaling behavior. Moreover the problem of having three equations for the two parameters $`g`$ and $`b`$, which invalidates the asymptotic analysis of des Cloizeaux, is naturally solved within the new ansatz we propose in the next subsection. The value of $`w`$ is no longer constrained to a certain value as it should be the case in an ideal experiment where it is allowed to tune the Coulombic interaction (or the quality of the solvent for the self avoiding problem).

### C Corrections to scaling

As it was suggested we will stay within the Flory balancing condition but try the next asymptotic ansatz for the correlation function (14)

$$b(n)b_0n^{2\nu }(\mathrm{log}n)^{2\gamma },$$
(38)
where $`\gamma `$ is the exponent that characterize the logarithmic correction to scaling. Substituting the ansatz (38) in the Euler equation (22) gives for the interaction term

$$I(k)=wb_0^{\frac{d+2}{2}}\underset{n=1}{\overset{\mathrm{}}{}}\frac{1\mathrm{cos}(nk)}{n^{\nu (d+2)}[\mathrm{log}n]^{\gamma (d+2)}}.$$
(39)
The largest contribution to this series at small $`k`$ comes from the large $`n`$ limit. We can substitute the sum with an integral (for the power law ansatz, this is explained and justified in the Appendix) to obtain

$$I(k)wb_0^{\frac{d+2}{2}}_{\overline{n}}^{\mathrm{}}𝑑n\frac{1\mathrm{cos}(nk)}{n^{\nu (d+2)}[\mathrm{log}n]^{\gamma (d+2)}},$$
(40)
where $`\overline{n}`$ is a cutoff we introduce to regularize the integral. Equation (40), after substituting $`n=x/k`$ reads

$$I(k)k^{\nu (d+2)1}wb_0^{\frac{d+2}{2}}_{\overline{x}}^{\mathrm{}}𝑑x\frac{1\mathrm{cos}(x)}{x^{\nu (d+2)}[\mathrm{log}x\mathrm{log}k]^{\gamma (d+2)}},$$
(41)
where, because of the Flory balancing conditions ($`I(k)k^2`$), we should impose $`\nu (d+2)=3`$. Let us consider a value $`t`$ such that at $`x<t`$ it is ensured that $`1\mathrm{cos}xx^2/2`$ holds. In this case the integral in equation (41) becomes

$`J(k)`$ $`=`$ $`{\displaystyle _{\overline{x}}^{\mathrm{}}}𝑑x{\displaystyle \frac{1\mathrm{cos}(x)}{x^3[\mathrm{log}x\mathrm{log}k]^{\gamma (d+2)}}}`$ (42)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{k\overline{n}}^t}𝑑x{\displaystyle \frac{1}{x[\mathrm{log}x\mathrm{log}k]^{\gamma (d+2)}}}`$ (43)
$`+`$ $`{\displaystyle _t^{\mathrm{}}}𝑑x{\displaystyle \frac{1\mathrm{cos}(x)}{x^3[\mathrm{log}x\mathrm{log}k]^{\gamma (d+2)}}}=J_1+J_2.`$ (44)
At $`k0`$ $`J_2(k)0`$, while the integral $`J_1(k)`$ defined above can be evaluated as

$$J_1(k)=\frac{1}{2[\gamma (d+2)1]}[[\mathrm{log}(\overline{n})]^{1\gamma (d+2)}],$$
(45)
provided that $`\gamma (d+2)>1`$ and $`\overline{n}>1`$. We conclude that the Euler equation (22) is properly regularized by the logarithmic ansatz (38) and the cutoff parameter $`\overline{n}>1`$. The balance in equation (22), between the elastic energy term and the interaction term, gives

$`\nu `$ $`=`$ $`{\displaystyle \frac{3}{d+2}},`$ (46)
$`\gamma `$ $`>`$ $`{\displaystyle \frac{1}{d+2}}.`$ (47)
### D Correction to scaling - charged chains

For the Coulombic case, considering the dimensional shift rule (21), we obtain

$`\nu `$ $`=`$ $`{\displaystyle \frac{3}{d}},`$ (48)
$`\gamma `$ $`>`$ $`{\displaystyle \frac{1}{d}}.`$ (49)
The inequality we just derived will be verified numerically in the next section. Values for the exponent $`\gamma `$ will be obtained via the numerical solution of the Euler equations (12)-(17) and via Monte Carlo simulations. Within our asymptotic analysis, we can obtain an expression for the exponent $`\gamma `$ that depends on the interaction strength. In particular by balancing the prefactors of the proper terms we have

$$\gamma \frac{1}{d}\left(1+\frac{w}{b_0(w)^{d/2}}\right),$$
(50)
where we also took into account that the amplitude $`b_0`$ can depend from $`w`$. The exponent $`\gamma `$ depends on the interaction strength, being non universal. This will be tested against simulations and numerical minimization of the variational free energy (12).

It was suggested  that above a certain minimum charge fraction $`fN^{3/4}(b/l_\mathrm{B})^{1/2}`$ the electrostatic interactions are relevant and the end-to-end distance can be written in the form

$$RNf^{2/3}(l_\mathrm{B}b^2)(\mathrm{log}N)^\gamma .$$
(51)
In this equation the value of $`\gamma =1/3`$ is confirmed also by a simplified chain under tension variational method  where the trial Hamiltonian is the Hamiltonian of a Gaussian chain subject to an external tension $`\stackrel{}{F}`$. Indeed we also reproduce this result for large enough values of the Coulomb strength as we will discuss in the next section. We observe that the exponent $`\gamma `$ is non universal, i.e. decreasing with the interaction strength, and approaching the limiting value $`\gamma =1/3`$ from above in the region of $`\beta >1`$, in agreement with equation (49). At $`\beta =1`$, the value of $`\gamma 0.38`$ is in good agreement with the calculation of Söderberg and collaborators . It is important to remark that if we assume in the previous asymptotic analysis that the $`\gamma `$ exponent is universal a spurious nested logarithmic form for the end-to-end distance (the way Allegra  discussed long time ago) is obtained. This scenario however is not confirmed by the numerical analysis of the Euler equations as will be presented in section III. The balancing between the elastic energy and the interaction term determines the exponents (46)-(47) for the self avoiding chain (or short-range interaction) problem, whereas for the Coulombic case the result is given in equations(48)- (49).

The term that appears on the l.h.s. of equation (30) is negligible in this case. The proper condition for this to happen reads $`1+2\nu >2`$, which leads to the condition $`d<4`$ for the self-avoiding chain  and $`d<6`$ for the Coulombic chain . Accordingly the upper critical dimensionality is $`d_c=4`$ for the former case and $`d_c=6`$ for the latter one. At $`d>d_c`$ we can neglect the interaction term and the chain becomes Gaussian.

We would like to make a general remark at this point. It is clear that the variational technique cannot cover the results of the RG. From the universality of the self avoiding walk problem it is well known that the choice of the parameter $`w`$ does not change the scaling exponent and below the upper critical dimension no logarithmic corrections are necessary. At the critical dimension, these become obviously real. On the other hand the variational technique produce these corrections naturally and works well in the mean field limit. For the case of polyelectrolytes the situation is more complex. The Coulomb potential changes its nature from long range to short range between three and four dimensions. Therefore the fixed point and scaling behavior is completely different, as it has been discussed already by de Gennes . It is therefore also necessary to check our results by alternative methods.

## III Numerical Solution of the Euler equations

In this section we consider in detail the Euler equations (13) - (17) for a polymer chain of length $`N`$ with cyclic boundary conditions, obtained by minimization of the variational free energy $`F_V`$. Two different ways in order to solve the variational principle minimization for a finite chain of length $`N`$ are possible. On the one side, a multidimensional minimization of the free energy $`F_V`$ in a $`N`$ dimensional space is considered . On the other side, a direct solution of the Euler equations will be presented.

We checked first that our multidimensional algorithm is able to reproduce the simple Gaussian profile for the simple case of the non interacting chain problem. We then implement our algorithm by choosing a very small value of the interaction strength $`\beta `$, where the solution is expected to be very close to the simple Gaussian chain solution, and give $`N+1`$ trial profiles for the function $`g_N(k)`$. The multidimensional algorithm via a sequence of reflections and translations between the $`N+1`$ trial profiles $`g_N(k)`$, relaxes to the true solution that can be resolved within the numerical precision desired. Once the algorithm was properly tuned, we moved on studying systematically the shape of the optimal profile that minimizes the free energy $`F_V`$ for increasing values of the interaction strength $`\beta `$. For any specific value of $`\beta `$ the optimal profile was checked to satisfy the basic symmetry requirements given in equation (16). Once the numerical solution $`g_N(k)`$ that minimize the variational free energy was found, we proceeded to analyze the properties of the solution in the asymptotic regime, in order to check our results for the maximal mean square distance $`b(N/2)`$. Results for the mean square distance $`(\stackrel{}{r}_{i+n}\stackrel{}{r}_i)^2`$ for intermediate values of $`n`$ were also studied in order to confirm the phenomenological predictions of the electrostatic blob picture.

A direct implementation of the equations (13) - (17) does confirm the multidimensional minimization of the free energy $`F_V`$ and also represents a convenient method to solve the variational problem itself. We also consider, for the case of $`\lambda =1`$ a Fast Fourier Algorithm that let us consider chains of length up to $`N=2^{15}`$. The difference between the values of the exponent $`\gamma `$ computed with $`N<1024`$ and $`N<32768`$ are very small and tends to vanish when the interaction strength increases. We consider this as a second, independent test to confirm the non-universal behavior of $`\gamma `$. For this reason we will consider shorter chains only for what concerns the case of $`\lambda =2`$ and for the Self Avoiding Chain Problem (see section IIc). We start with a purely Gaussian profile in order to compute the left hand side of equation (17), via the Fourier transform of (15) that defines the mean square distance $`b(n)=(\stackrel{}{r_{i+n}}\stackrel{}{r_i})^2`$ for a small value of the interaction strength $`\beta `$ and iterate this procedure until convergence is achieved. Once the solution of the Euler equations is found for a specific value of the interaction, one can use it as a starting point for the following iteration procedure, in order to solve the same set of equations for a higher value of $`\beta `$. It is important to note that a simple recursive implementation of the Euler equations above is highly unstable in the sense that the number of iterations required as a function of $`\beta `$ rapidly increases, making a simple iteration procedure quite ineffective. A simple solution to this problem can be found borrowing ideas from the renormalization group literature . Renormalization Group Theory was successfully applied to determine the phase diagram and magnetic properties of disordered spin systems . A functional set of recursion relations for the probability distribution function of bond randomness was studied and a fixed distribution was found for different values of the initial conditions. Instead of a simple iterative procedure, which could be a rather slow converging process, a multidimensional bisection between different initial conditions for the renormalization group trajectory turns out to be a very effective method in order to achieve convergence. In a similar fashion, instead of a simple iteration of the Euler equations above, we prefer to bisect between possible initial profiles of the function $`g(k)`$ given above.

### A Results for the Coulombic chain

We obtained the solution of the Euler equations (13)-(17) via the iterative method discussed above and independently via the direct minimization of the variational free energy $`F_V`$ in equation (12) within the numerical precision . We considered cyclic chains of length $`2^4<N<2^{15}=32768`$, for different values of the interaction strength $`\beta `$ and two different forms of the potential (9), corresponding to the values $`\lambda =1,2`$. We report the results of the end-to-end distance for different values of $`\lambda `$ and for chains of increasing size $`N`$ at values of $`\beta =0.5`$ and $`\beta =1.0`$.

| $`\lambda `$ | $`\beta `$ | N=120 | N=240 | N=320 | N=520 | N=1024 | N=2048 | N=4096 | N=8192 | N=16384 | N=32768 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 1.0 | 0.5 | 60.6 | 128.2 | 174.5 | 293.2 | 601.7 | 1249.4 | 2583.4 | 5372.2 | 11071.2 | 22769.8 |
| 1.0 | 1.0 | 78.0 | 164.4 | 223.6 | 374.2 | 767.59 | 1590.6 | 3283.9 | 6778.8 | 13942.4 | 28627.3 |

| $`\lambda `$ | $`\beta `$ | N=120 | N=240 | N=320 | N=520 | N=1024 |
| --- | --- | --- | --- | --- | --- | --- |
| 2.0 | 0.5 | 26.2 | 48.6 | 63.0 | 98.12 | 140.2 |
| 2.0 | 1.0 | 33.2 | 62.5 | 81.4 | 126.2 | 248.7 |

Tab.1: Results for the end-to-end distance $`R`$ for different values of the chain length $`N`$ and the interaction strength $`\beta `$ and $`\lambda `$. These results are obtained solving explicitly the the Euler equations (13)-(17).

Let us begin to discuss the case of a simple polyelectrolyte chain in $`d=3`$ with pure Coulombic interaction ($`\lambda =1`$). We analyzed the end-to-end distance behavior for increasing values of $`N`$, i.e. $`b(N/2)`$. This is reported in Fig.1. Dependence of the exponent $`\gamma `$ on the interaction strength suggests a weak violation of universality as manifest in Fig.2.

The des Cloizeaux argument, once extended to the polyelectrolyte chain problem in $`d=3`$ implies an overstretching ($`\nu =2/\lambda =2`$), while the ansatz suggested above indicates $`\nu =3/(\lambda +2)=1`$. In order to compute the non-universal exponent $`\gamma `$ we plot $`R/N^{3/(\lambda +2)}`$ as shown in the inset of Fig.1. The non-universal behavior of $`\gamma `$ agrees with the inequality (49) derived in the asymptotic analysis of section II. In Fig.2, the exponent $`\gamma `$ decreases with increasing values of $`\beta `$, approaching values very close to $`\gamma =1/3`$, for large values of the interaction strength. A recent discussion about the logarithmic correction to scaling behavior for polyelectrolyte chains, has been given  and a direct derivation of the value $`\gamma =1/3`$ has been obtained. In the notation of reference  the small temperature limit $`T0`$ corresponds to the large interaction strength limit $`\beta \mathrm{}`$ which is consistent with our findings. The result

$$RN(\mathrm{log}N)^{1/3}$$
(52)
obtained in the large $`\beta \mathrm{}`$ limit, can be interpreted as a stretching of the harmonic bonds of the chain. The value of $`\gamma =1/3`$ was also obtained by a simplified chain under tension variational method  where the trial Hamiltonian is the Hamiltonian of a Gaussian chain subject to an external tension $`\stackrel{}{F}`$. This result, confirmed also by the Monte Carlo simulations presented in the next section, holds in the limit of large interaction strength but should not considered as a limiting case of little practical impact. Ordinary temperatures in the aqueous solutions correspond to surprisingly “high” values of the interaction strength $`\beta `$ .

A rather important problem discussed in the past  is the polyelectrolyte chain with an interaction (9) at $`\lambda =2`$. After the proper rescaling of the interaction strength $`\beta `$ (see equation 12), this problem can be thought as a chain with pure Coulombic interaction for $`d=4`$ . Results for the end-to-end distance versus the chain length $`N`$ are shown in Fig.3.

A similar discussion about the non-universal behavior of the exponent $`\gamma `$, follows the same lines as in the previous case (Fig.1-Fig.2). According to the ansatz we proposed in the asymptotic analysis of section II C we compute the exponent $`\gamma `$ as in Fig.4, where $`\nu =\nu _F`$ is the Flory exponent.

In Fig.5 we compute the non-universal behavior of the exponent $`\gamma `$ assuming a linear dependence of the end-to-end distance of the chain length $`N`$, beside logarithmic terms. This is not the case, for the potential (9) with $`\lambda =2`$, but interesting comparison with previous results follows. The exponent $`\gamma `$, as in Fig.5, approaches values of $`\gamma 0.4`$ at values of $`\beta =1.0`$ where Monte Carlo calculations were performed  and where the computed value of $`\gamma 0.35`$ is in good agreement with our result. A two parameter fit would represent a very conclusive way to discriminate between these two possibilities, but the chain lengths that one can reach in simulations and in a direct minimization of the variational free energy are not long enough to make such analysis effective.

### B Electrostatic blob picture

We will now present a numerical analysis of the electrostatic blob picture. We measured the crossover (shown by the mean square distance as a function of the polymer chain length $`N`$) from a Gaussian regime, observed at small distances, to the stretched regime for long enough distances. This represents an effective way to compute the electrostatic blob size and to check our method. One of the major achievement of the Gaussian variational approach discussed in section I is that it does not simply provides a quantitative description of the large $`n`$ behavior correlation function. A full numerical solution of the Euler equations (13)-(17) provides informations on the intermediate behavior of the correlation functions $`b(n)=(\stackrel{}{r_{i+n}}\stackrel{}{r_i})^2`$. That is why it is interesting to check whether short enough segments of the polymer chains behave as Gaussian chains, as expected from general scaling arguments , whereas the large $`n`$ behavior corresponds to equation (51). Once the crossover is observed, the size of an electrostatic blob $`\xi `$ can be extrapolated. The dependence of the electrostatic blob size $`\xi `$ from the interaction strength $`\beta `$, that we expect from phenomenological arguments, is indeed confirmed within our variational approach (see Fig.6). In particular we expect

$`R`$ $``$ $`(N/g)\xi Na\beta ^{1/3}`$ (53)
$`\xi `$ $``$ $`a\beta ^{1/3}`$ (54)
In Fig. 6 we show the electrostatic blob size as a function of the interaction strength $`\beta `$. The $`\xi \beta ^{1/3}`$ behavior expected from standard scaling arguments is here also found.

### C Results for the Self Avoiding Chain problem

We consider the problem of a polymer chain with a short range interaction between the monomers, in a good solvent. Within the Gaussian Variational Principle, the expression for the free energy and the Euler equation are given by equations (20) and (22). In section II we have discussed the dimensional shift (21). Nevertheless this equivalence is purely formal and no deep insight is hidden behind it. The two problems considered are indeed pretty different and they reduce to the same type of variational equations because of the Gaussian approximation itself. It is important to remark that the Gaussian Variational Principle, specially reliable for the long ranged polyelectrolyte problem, has been originally introduced to study the self avoiding homopolymer chain problem . Even though it is nowadays well established that the Gaussian Variational Principle is not very well suited for the short ranged interacting chain problem, we wish to make clear that within this approximation, apart for spurious logarithmic corrections, the usual Flory-like power law behavior should be recovered.

In the original asymptotic analysis of des Cloizeaux the non integer value of $`\delta `$ is tacitly assumed, so that the result of $`\nu =2/d`$ for the swelling exponent was obtained. Although we believe this result to be incorrect, it is often mentioned in the literature to discuss the asymptotic properties of polyelectrolyte chains . According to our analysis and to the previous analysis of Allegra and coworkers , the problem of a self avoiding chain in good solvent turns out, within the Gaussian Variational Principle, to present the expected Flory-like behavior. This is not surprising since the original derivation of Flory also belong to the class of Gaussian Variational approximations, since a Gaussian probability distribution of the distance between the ends of the chain is assumed. A numerical justification of the des Cloizeaux result for the swelling exponent can also be found in the literature . We believe that this numerical analysis is inconsistent since the constraint on the interaction parameter (assumed by des Cloizeaux and critically reviewed in section II) is enforced in the numerical algorithm. This makes the results obtained highly questionable. In Fig. 7 and Fig. 8 we present the results of the end-to-end distance for the self avoiding homopolymer chain problem studied via the numerical solution of the Euler equations discussed above.

Results for the end-to-end distance as a function of the polymer chain length $`N`$ are shown and a quantitative analysis of the non universal logarithmic part to be considered as a spurious effect due to the short range nature of the problem is presented.

According to the asymptotic analysis given above we expect the exponent $`\gamma `$ to be non universal and to depend indeed on the quality of the solvent.

## IV Numerical Simulations

In order to test the conclusions obtained in the previous section we decided to perform a numerical simulation for a polyelectrolyte chain in an ideal salt free solution, when pure Coulombic interactions are considered between the monomers. A bit of caution is required since we will compare the results of our simulations with the results obtained via the Gaussian variational principle, where a cyclic geometry for the polymer chain is chosen. It was easy however to device a numerical simulation of the Monte Carlo type where close chain configurations can be taken under investigation. As remarked long time ago,  the cyclic geometry chosen within the variational calculation becomes unessential as soon as long enough chains are considered and this turns out to be confirmed by our numerical simulations. We then consider a Monte Carlo simulation, performed in the canonical ensemble via the standard Metropolis algorithm. We consider a pivot move for what concerns the elementary Monte Carlo move, i.e. we assume rigid bonds between the monomers and an elementary pivot move consists in a random rigid rotation of a randomly chosen part of the chain. The pivot algorithm, first described in reference  and lately applied in the context of polyelectrolyte chains with and without salt  is a very effective method to confirm our general conclusions. Several million Monte Carlo moves (the number which was properly increased with the chain length) are required in order to reduce below $`1\%`$ percent the statistical fluctuations for the end-to-end distance. We also combine our Monte Carlo simulation with a Molecular dynamics subroutine that is randomly called, during the Monte Carlo simulation, with a certain fractional occurrence ratio. On the one side, the pivot algorithm is well suited to simulate the polymer chain at long enough scales. The pivot moves are in a sense collective moves and a good description of the chain is certainly obtained at long scales. On the other side the Molecular Dynamic algorithm is well suited to simulate the chain at small distances. We initially perform independent simulations with the Monte Carlo and the Molecular dynamics algorithm in order to tune the simulation parameters involved for a chain of a given length $`N`$. Once the two methods were independently tested, we proceeded to combine them as discussed above.

We now present our results for the end-to-end distance for $`d=3`$ polyelectrolyte problem. In Fig. 9 the end-to-end distance is shown as a function of the polymer length $`N`$ for the specific value of $`\beta =1.0`$. The inset shows the end-to-end distance divided by the pure power law Flory-like dependence $`N^\nu `$ as a function of $`\mathrm{log}N`$. We repeated the same analysis of Fig.9 changing the value of the interaction strength $`\beta `$. The non universal behavior obtained in the asymptotic analysis of section II and in the numerical analysis of the Euler equations (13) - (17), is found within our simulation results (see Fig.10).

The non universal exponent $`\gamma `$ is seen to decrease, as we obtained with the variational method of section II-III. Nevertheless one should treat the results of simulation with a bit caution.

In the numerical simulations monomers are considered to be rigidly connected while in the variational principle approach monomers are connected by springs. In this case one can not expect that the exponent $`\gamma `$ as a function of interaction strength has the same large $`\beta `$ behavior.

## V Conclusions

In the present paper we have reinvestigated the Gaussian variational principle and its applicability to polymer chain with long and short range interactions. The asymptotic analysis originally presented by des Cloizeaux in the context of self avoiding walks has been critically reconsidered. In particularly we have shown that the des Cloizeaux’s analysis overconstains the model’s parameters. It has been found that within the Flory type of balancing conditions (when the elastic term is balanced by the interaction but the entropic term is not relevant) the asymptotic law (38) satisfies the Euler equation (17) . In this asymptotic law the Flory exponent $`\nu =3/(d+2)`$ and logarithmic corrections to scaling is characterized by a non - universal exponent $`\gamma `$. Within the Gaussian variational principle the universal exponents for the case of the Coulombic and short range interactions are related by a dimensional shift (21).

Gaussian variational principle suits much better for the long range problems which was already discussed in reference . The arguments we have presented prove that the scaling form (38) is valid for a general long - range interaction potential (9). The exponents $`\nu `$ and $`\gamma `$ satisfy the equations (48) - (49). According to the condition (49), the exponent $`\gamma `$ is non - universal and depends from the interaction parameter as well as from the chain model (compare e.g. Fig.2 and Fig. 10) . In this sense we may consider the combination $`l_{\mathrm{eff}}=a(\mathrm{log}N)^\gamma `$ as an effective Kuhn segment for chains with unscreened type of interactions.

On the contrary, for the short range interactions the Gaussian variational principle simply overestimates the probability of the segment - to - segment contacts  which leads to more extended configurations compare to the Monte Carlo results . We conclude observing that the models which are characterized by the interplay between long - ranged attraction and repulsion are also useful systems to be treated by the Gaussian variational principle, as presented here. Appropriate systems are polyampholyte . For the randomly charged polyampholytes Kantor and Kardar  have found strongly nonuniversal behavior. For example at $`d<4`$ the spatial configurations depend on the polyampholyte excess charge $`Q`$ so that if $`Q>Q^{}e\sqrt{N}`$ the compact form changes to the stretched one. We believe that the Gaussian variational principle is a good starting point for the investigation of such systems.

## Appendix A

The asymptotic analysis of equations (22) suggests to consider the following series

$$S(\delta ,z)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{z^n}{n^\delta },$$
(55)
In the asymptotic limit $`k0`$ the interaction term of equation (22) is simply related to $`S(\delta ,z)`$. Even the simple power law ansatz $`b(n)n^\nu `$ requires some care. The expression ( 55) is evaluated exactly but two rather different expressions are obtained if the exponent $`\delta =\nu (d+2)`$ is assumed to be integer or non-integer. The case of the exponent $`\delta `$ being non-integer reduces the interaction term to the expression given in equation (29) and leads to the des Cloizeaux result. Assuming an integer value for the exponent $`\delta `$ produces instead a different expression for the interaction term and the proper balancing condition with the elastic energy naturally follows.

Let us consider the series $`S(\delta ,z)`$ in this two distinct cases. In equation (55) $`z`$ and $`\delta `$ can be complex quantities, as soon as $`|z|1`$ and $`Re\delta >1`$ at $`|z|=1`$ whereas $`Re\delta >0`$ if $`|z|<1`$. The calculation of the series (55) ( see sec. 1.11 in ref. ) gives

$$S(\delta ,z)=\mathrm{\Gamma }(1\delta )[\mathrm{log}(\frac{1}{z})]^{\delta 1}+\underset{r=0}{\overset{\mathrm{}}{}}\zeta (\delta r)\frac{(\mathrm{log}z)^r}{r!},$$
(56)
where $`|\mathrm{log}z|<2\pi `$, $`\mathrm{\Gamma }(x)`$ is the gamma function, $`\delta 1,2,..`$ and $`\zeta (x)`$ is the $`\zeta `$-zeta Riemann’s function. We are interested in the case of $`z=\mathrm{exp}(ik)`$ at $`k0`$. Considering the real part of both sides of equation (27) in the asymptotic limit $`k0`$, we obtain equation (29).

We now consider the case of integer values of $`\delta `$. The $`\mathrm{\Gamma }(x)`$ has poles at all negative integer arguments whereas the pole of $`\zeta (x)`$ is placed at $`x=1`$. We write $`\delta =m+ϵ`$ where $`m`$ is a positive integer and $`ϵ0`$. Then in the vicinity of poles the gamma and $`\zeta `$ \- functions can be rewritten as

$`\mathrm{\Gamma }(1mϵ)`$ $`=`$ $`{\displaystyle \frac{(1)^m}{(m1)!}}\left\{{\displaystyle \frac{1}{ϵ}}\psi (m)+𝒪(ϵ)\right\}`$ (57)
$`\zeta (1+ϵ)`$ $`=`$ $`\left\{{\displaystyle \frac{1}{ϵ}}\psi (1)+𝒪(ϵ)\right\}.`$ (58)
We should also take into account that

$`\left[\mathrm{log}\left({\displaystyle \frac{1}{z}}\right)\right]^ϵ=1+ϵ\mathrm{log}\mathrm{log}\left({\displaystyle \frac{1}{z}}\right)+𝒪(ϵ^2).`$ (59)
After that the equation (56), due to the cancellation of poles in gamma and $`\zeta `$ \- functions at small values of $`ϵ`$, becomes

$`S(m,z)`$ $`=`$ $`{\displaystyle \frac{(\mathrm{log}z)^{m1}}{(m1)!}}\left(\psi (m)\psi (1)\mathrm{log}\mathrm{log}(1/z)\right)`$ (60)
$`+`$ $`{\displaystyle \stackrel{}{}_{r=0}^{\mathrm{}}}\zeta (mr){\displaystyle \frac{(\mathrm{log}z)^r}{r!}},`$ (61)
where the prime indicates that the term $`r=m1`$ is to be omitted. In equations (61) (58) $`\psi (x)`$ is the digamma function, defined as the logarithmic derivative of the gamma function $`\psi (\delta )=\frac{d\mathrm{log}\mathrm{\Gamma }(\delta )}{d\delta }`$. In the limit $`k0`$, for $`z=\mathrm{exp}(ik)`$ and $`\delta =3`$, eq. (61) have been considered in eq.(36).

We want to show now that the integral representation of the sum is correct in the asymptotic limit and that finite $`n<l`$ terms do not contribute in this limit. Let us consider explicitly the presence of a finite cutoff in the series (55).

$$S(\delta ,z)=\underset{n=l}{\overset{\mathrm{}}{}}\frac{z^n}{n^\delta },$$
(62)
which reduces to equation (55) at $`l=1`$. For integer values of $`\delta `$ it is easy to derive exactly the following expression :

$`S(m,z,l)`$ $`=`$ $`{\displaystyle \frac{(\mathrm{log}z)^{m1}}{(m1)!}}\left(\psi (m)\psi (l)\mathrm{log}\mathrm{log}(1/z)\right)`$ (63)
$`+`$ $`{\displaystyle \stackrel{}{}_{r=0}^{\mathrm{}}}\zeta (mr){\displaystyle \frac{(\mathrm{log}z)^r}{r!}}.`$ (64)
Consider the series (62) with a finite cutoff $`l`$ in the asymptotic limit $`k0`$. We want to evaluate the contributions to the series from finite $`n<l`$ terms, when $`m=3`$. We obtain

$$\underset{n=l}{\overset{\mathrm{}}{}}\frac{1\mathrm{cos}(nk)}{n^3}\frac{k^2}{2}(\psi (3)\psi (l)\mathrm{log}k)+O(k^4),$$
(65)
so that the contribution from finite $`n`$ takes the following form

$$\underset{n=1}{\overset{l}{}}\frac{1\mathrm{cos}(nk)}{n^3}\frac{k^2}{2}(\psi (l)\psi (1))$$
(66)
Summing up equations (65) and (66) we obtain again the result (36). The expression (65) shows that for small $`k`$ values only large $`n`$ terms contribute to the sum (62). Indeed the contribution coming from $`n>l`$ produces a term $`k^2\mathrm{log}(k)`$. This, whenever $`k<<\mathrm{exp}(\psi (l))`$, will dominate the simple $`k^2`$ behavior. Moreover, for large values of $`l`$ (see ref. ) we have $`\psi (l)\mathrm{log}(l)`$ and the condition above takes the form $`k<<1/l`$. It follows that for $`k<<1/l`$ the term $`k^2\mathrm{log}(k)`$ dominates and the sum does not depends on the low limit value. In this case, the series can be well approximated by an integral, i.e.

$$S(k)=\underset{n=l}{\overset{\mathrm{}}{}}\frac{1\mathrm{cos}(nk)}{n^3}_l^{\mathrm{}}𝑑n\frac{1\mathrm{cos}(nk)}{n^3}=k^2_{kl}^{2\pi }𝑑x\frac{1\mathrm{cos}(x)}{x^3}=k^2I(k)$$
(67)
In order to compute $`I(k)`$, consider the derivative of the previous expression, $`I^{}(kl)=1/2kl`$ so that $`I(kl)\frac{1}{2}\mathrm{log}(kl)+c`$, where $`c`$ is a constant, and $`S(k)=\frac{k^2}{2}(c\mathrm{log}(kl))+O(k^4)`$. Comparison with equation (65) shows that for $`k<<1/l`$ only large $`n`$ contributions are relevant.

## Acknowledgments

We thank Christian Holm, Burkhard Dünweg and Ralf Everaers for continuous and fruitful discussions. In particular we thank Dr. Holm for several comments on the simulation part. One of us (G.M.) wish to express many thanks to Giorgio Parisi and Roland Netz.