text:
# Dirac Decomposition of Wheeler-DeWitt Equation on the Bianchi Class A Models

## 1 Introduction

Arnowitt-Deser-Misner (ADM) reformulated the general relativity using the canonical formalism, which makes a 3+1 split of the space-time metric, and introduced some constraints through the variational principle. According to the ADM approach, the Einstein equation is derived as the equation of dynamical evolution. The canonical formalism by the ADM approach is a system of dynamics including two constraints: one is the momentum constraint and the other is the Hamiltonian constraint. The latter governs the time evolution of the system. According to Dirac  the procedure of canonical quantization replaces the Hamiltonian constraint with the supplementary condition on the wave function which represents the quantum state of space-time. The Hamiltonian constraint thus becomes the Wheeler-DeWitt equation, which is fundamental to the quantum gravity and determines the quantum state of space-time. The Wheeler-DeWitt equations for specific cosmological models like the Friedmann-Robertson-Walker model  or the Bianchi models  have been investigated.

The Wheeler-DeWitt equation for the Bianchi models is in general reduced to the second-order differential equation if the operator-ordering is neglected. The resultant Wheeler-DeWitt equation is similar to the Klein-Gordon type, which has a problem that the probability density is not positive-definite. The Wheeler-DeWitt equation has also often this kind of problem that the probability density becomes negative. The negative values of probability density have been shown by Furusawa  through numerical calculation, who treated the Wheeler-DeWitt equation for the quantum mixmaster model using wave-packet solution.

To avoid the problem, the approach called Dirac-Square-Root formalism has been studied . Since the Hamiltonian for the Bianchi type models is expressed in terms of the quadratic, the first-order differential equation can be derived by applying a procedure similar to the Dirac method. Mallett  applied it to the Friedmann model coupled to the charged scalar field. Kim and Oh  mainly discussed the Bianchi type-IX model. In their papers, the first-order equation included an unknown function dependent on the time-parameter is postulated and the wave function is required to be a solution of two-component column vector. Iterating the first-order equation thus postulated, the differential equation with respect to the unknown function is yielded so as to be consistent with the Wheeler-DeWitt equation for each component. The differential equation with respect to the time parameter is thus obtained. Indeed, the differential equation is reduced to Riccati or Bernoulli equation in case of the Friedmann model with charged scalar field or the Bianchi type-IX model, respectively. However, one sees that it is difficult to obtain the exact solutions in both cases.

The purpose of this paper is to propose a new method, which factorizes the Wheeler-DeWitt equation into the first-order differential equation, using the Pauli matrices. The resultant differential equation thus obtained is similar to the Dirac equation. As a necessary condition for the Dirac factorization, it is postulated that the resultant Hamiltonian is “self-adjoint”. We apply this method to the the Wheeler-DeWitt equation of the Bianchi Class A vacuum models except for type-IX model. The solution of the Dirac type equation is written by two-component column vector. Using the Dirac type equation and the complex conjugate, we can define the probability density. Then, we get the equation of continuity and the conserved current. Here, we restrict ourselves to determine the form of Hamiltonian so as to make the probability density positive-definite. It is to be noted that our formalism cannot be applied to the Bianchi type-IX model, because the structure of anisotropic potential for the Bianchi type-IX model is different from that of the other Bianchi Class A models.

The organization of this paper is as follows. In §2, we review the ADM canonical formulation for the vacuum Bianchi Class A models and introduce the Wheeler-DeWitt equation. In §3, we propose the new method in which the Wheeler-DeWitt equation is factorized into the first-order differential equation in terms of the Pauli matrices and apply the method to the Bianchi type-I, II, VI0, VII0 and VIII models. For each of them, the Dirac type equation is derived and solved except for the Bianchi type-VIII model. Then, we assume a trail function which is expressed in terms of two-component spinor. We show that the probability density becomes positive-definite and there is a conserved current. As the time parameter approaches the limit of the singularity of the universe, the Wheeler-DeWitt equation for the Bianchi type-II, VI0, VII0 and VIII models are reduced to that of the Bianchi type-I model, so that the Dirac type equations must possess the solution for the Bianchi type-I model. In §4, we discuss the behavior for the Dirac type equation near the singularity in detail. It is shown that from the Heisenberg equation the quantum behavior of the anisotropic parameters $`\beta _\pm `$ make the quivering motion and the Hamiltonian system has a ‘spin-like’ degree of freedom. The ‘spin-like’ degree of freedom causes the universe to go necessarily through an early quantum stage of evolution with agitated anisotropy-oscillation. In §5, summary and discussion are made on the results of our formalism and comparisons are made with the Dirac-Square-Root formalism. Henceforward, we take that $`i`$, $`j`$ and $`k`$ runs from 1 to 3 and the unit $`c=\mathrm{}=16\pi G=1`$.

## 2 Canonical Quantization formulation for the Bianchi Class A models

The Einstein-Hilbert action is given by

$$I=\sqrt{^{\left(4\right)}g}^{\left(4\right)}Rd^4x,$$
(1)
where $`^{\left(4\right)}g`$ and $`^{\left(4\right)}R`$ are the determinant of the space-time metric and the curvature, respectively; the indices $`(4)`$ implies 4-geometry quantity. The metric for the Bianchi models is generally given by

$$ds^2=N^2dt^2+e^{2\alpha }e^{2\beta _{ij}}\chi ^i\chi ^j,$$
(2)
where $`N`$, $`\alpha `$ and $`\beta _{ij}`$ are functions of $`t`$ only: $`N`$ is the lapse function; $`e^{2\alpha }`$ plays the role of a scale factor, so that the universe has the initial singularity as $`\alpha \mathrm{}`$; $`\beta _{ij}`$ are an isotropic parameters of the universe, which gives a matrix form of traceless and $`\chi ^i`$ are 1-forms. According to Misner , anisotropic parameters $`\beta _{ij}`$ are expressed in terms of two parameters $`\beta _+`$ and $`\beta _{}`$

$$\beta _{ij}=diag(\beta _++\sqrt{3}\beta _{},\beta _+\sqrt{3}\beta _{},2\beta _+).$$
(3)
The 1-forms obey

$$d\chi ^i=\frac{1}{2}C_{jk}^i\chi ^j\chi ^k,$$
(4)
where $`C_{jk}^i`$ are the structure constants.

The Bianchi models are characterized by the structure constants. The Bianchi models are classified into two groups from the structure constants, so-called Class A and B . It should be noted that, as was first indicated by Hawking , it cannot be reduced to the correct field equations by the ADM canonical formulation in the Bianchi Class B models except for the Bianchi type-V model . The reason is that the spatial divergence term does not vanish. Hence, we must discuss only the Bianchi Class A models because we use the canonical formulation. The characteristic of the Bianchi Class A models is that the structure constants obey the condition

$$C_{}^{i}{}_{ji}{}^{}=0.$$
(5)
We must therefore select the 1-forms so as to satisfy Eq. (5).

After the procedure for the ADM decomposition, we use the coordinates $`\alpha `$, $`\beta _+`$, $`\beta _{}`$ and the corresponding canonical momenta as $`p_\alpha ,p_+,p_{}`$. From the variation of $`N`$, one obtains the Hamiltonian constraint. The 3-geometry curvature $`^{\left(3\right)}R`$ is obtained from the structure constants. According to Misner , the relation of the curvature $`^{\left(3\right)}R`$ to the anisotropic potential $`V(\beta _+,\beta _{})`$ is introduced by

$$^{\left(3\right)}R=\frac{3}{2}e^{2\alpha }\left(1V(\beta _+,\beta _{})\right).$$
(6)
As a result, the Hamiltonian constraint becomes

$$H=p_\alpha ^2+p_+^2+p_{}^2+e^{4\alpha }(V1).$$
(7)
On the canonical quantization procedure, we treat the operators $`p_\alpha `$, $`p_+`$, $`p_{}`$ in the Schrödinger representation:

$$\widehat{p}_\alpha =i\frac{}{\alpha },\widehat{p}_+=i\frac{}{\beta _+},\widehat{p}_{}=i\frac{}{\beta _{}}.$$
(8)
Because the Hamiltonian constraint is the first-class constraint, the quantized Hamiltonian constraint must be satisfied by the supplementary condition

$$\widehat{H}\mathrm{\Psi }=0,$$
(9)
$$\widehat{H}\frac{^2}{\alpha ^2}\frac{^2}{\beta _+^2}\frac{^2}{\beta _{}^2}+e^{4\alpha }(V1),$$
(10)
where $`\mathrm{\Psi }=\mathrm{\Psi }(\alpha ,\beta _+,\beta _{})`$ is a state vector of the system, called by the wave function of the universe. This is the Wheeler-DeWitt equation for the Bianchi models.

## 3 Dirac factorization for the Wheeler-DeWitt equation

According to the Furusawa’s discussion , the probability density of Klein-Gordon type

$`\rho (\alpha ,\beta _+,\beta _{})={\displaystyle \frac{i}{2}}\left(\mathrm{\Psi }^{}_\alpha \mathrm{\Psi }\mathrm{\Psi }_\alpha \mathrm{\Psi }^{}\right)`$
becomes negative for some values $`\alpha `$. The negative values of the probability density are inconsistent with the physical theories and it is the failure of the structure in the Wheeler-DeWitt equation. To get rid of the problem, we adopt the similar procedure that Dirac factorized the Klein-Gordon equation in order to obtain the positive-definite probability density. That is, we try to factorize the Wheeler-DeWitt equation into the first-order differential equation. We denote the first-order equation as Dirac type equation. As a necessary condition of Dirac factorization, it is postulated that the resultant Hamiltonian is “self-adjoint”, because one tries to think various ways of the factorization into the first-order equation and there is not any rules for the factorization. We suppose that the Wheeler-DeWitt equation for the Bianchi models is replaced by the equation of Dirac type

$$\mathrm{\Psi }\left(\widehat{p}_\alpha +\sigma _1\widehat{p}_++\sigma _2\widehat{p}_{}+\sigma _3e^{2\alpha }v\right)\mathrm{\Psi }=0,$$
(11)
where $`v\sqrt{V1}`$ and the Pauli matrices $`\sigma _i`$ are expressed as

$`\sigma _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),\sigma _2=\left(\begin{array}{cc}0& i\\ i& 0\end{array}\right),\sigma _3=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right).`$
We here rewrite the Dirac type equation (11) in the Schrödinger form

$$i\frac{}{\alpha }\mathrm{\Psi }=\left[\sigma _1\widehat{p}_++\sigma _2\widehat{p}_{}+\sigma _3e^{2\alpha }v\right]\mathrm{\Psi }.$$
(13)
The self-adjointness of the Hamiltonian is easily proved:

$$^{}=\widehat{p}_\alpha +\sigma _1\widehat{p}_++\sigma _2\widehat{p}_{}+\sigma _3e^{2\alpha }v=,$$
(14)
where the momentum operators are self-adjoint and commute with the Pauli matrices.

With the help of relations (8), the Dirac type equation (11) is rewritten as

$$i\frac{\mathrm{\Psi }}{\alpha }i\sigma _1\frac{\mathrm{\Psi }}{\beta _+}i\sigma _2\frac{\mathrm{\Psi }}{\beta _{}}+\sigma _3e^{2\alpha }v\mathrm{\Psi }=0.$$
(15)
Taking the complex conjugate and multiplying $`\mathrm{\Psi }^{}`$ to Eq. (15) from the left-hand side, we have

$$i\mathrm{\Psi }^{}\frac{\mathrm{\Psi }}{\alpha }+i\mathrm{\Psi }^{}\sigma _1\frac{\mathrm{\Psi }}{\beta _+}+i\mathrm{\Psi }^{}\sigma _2\frac{\mathrm{\Psi }}{\beta _{}}+\mathrm{\Psi }^{}\sigma _3e^{2\alpha }v\mathrm{\Psi }=0.$$
(16)
Subtracting Eq. (16) from Eq. (15) multiplied $`\mathrm{\Psi }^{}`$ from the left-hand side, we obtain

$$\frac{}{\alpha }\mathrm{\Psi }^{}\mathrm{\Psi }+\frac{}{\beta _+}\left(\mathrm{\Psi }^{}\sigma _1\mathrm{\Psi }\right)+\frac{}{\beta _{}}\left(\mathrm{\Psi }^{}\sigma _2\mathrm{\Psi }\right)=0.$$
(17)
Setting the probability density $`\rho `$ and the current vectors $`j_+`$ and $`j_{}`$ as

$$\rho =\mathrm{\Psi }^{}\mathrm{\Psi },j_+=\mathrm{\Psi }^{}\sigma _1\mathrm{\Psi },j_{}=\mathrm{\Psi }^{}\sigma _2\mathrm{\Psi },$$
(18)
we get the equation of continuity

$$\frac{\rho }{\alpha }+\frac{j_+}{\beta _+}+\frac{j_{}}{\beta _{}}=0.$$
(19)
Then, one sees that the probability density is positive-definite. We apply our formalism to solve the Bianchi type-I, II, VI0, VII0 and VIII models below.

### 3.1 The Bianchi type-I model

The metric of the Bianchi type-I model is given in Eq. (2), in which the 1-forms are given by

$$\chi ^1=dx^1,\chi ^2=dx^2,\chi ^3=dx^3.$$
(20)
The structure constants in the model obviously supply

$$C_{}^{i}{}_{jk}{}^{}=0.$$
(21)
It leads to the zero curvature of 3-dimensional hypersurfaces $`^{\left(3\right)}R_\mathrm{I}=0`$, which implies the spatial flat structure. From this, the anisotropic potential becomes $`V_\mathrm{I}=1`$. Then, the Wheeler-DeWitt equation is reduced to

$$\left(\widehat{p}_\alpha ^2+\widehat{p}_+^2+\widehat{p}_{}^2\right)\mathrm{\Psi }=0.$$
(22)
We postulate here that the Dirac type equation for the Bianchi type-I model holds:

$$_\mathrm{I}\mathrm{\Psi }_\mathrm{I}(\alpha ,\beta _+,\beta _{})=\left(\widehat{p}_\alpha +\sigma _1\widehat{p}_++\sigma _2\widehat{p}_{}\right)\mathrm{\Psi }=0.$$
(23)
The solutions to this equation have been already found . They are plane wave solutions given by two-component spinor form

$$\mathrm{\Psi }_\mathrm{I}(\alpha ,\beta _+,\beta _{})=e^{i(k_+\beta _++k_{}\beta _{}\omega \alpha )}\left(\begin{array}{c}\psi _1\\ \psi _2\end{array}\right),$$
(24)
where $`\omega `$ and $`k_\pm `$ are arbitrary real constants. Substituting the solution (24) into the Dirac type equation (23), we get the components

$$\psi _1=\frac{k_+ik_{}}{\omega }\psi _2\mathrm{and}\psi _2=\frac{k_++ik_{}}{\omega }\psi _1.$$
(25)
It implies obviously $`|\psi _1|^2=|\psi _2|^2`$ and the relation among $`\omega `$ and $`k_\pm `$ is

$$\omega =\pm \sqrt{k_+^2+k_{}^2}.$$
(26)
Thus, the parameter $`\omega `$ corresponding to the frequency of the plane wave depends on the parameters $`k_\pm `$ corresponding to the wave number, which vary continuously.

The probability density is

$$\rho _\mathrm{I}(\alpha ,\beta _+,\beta _{})=\mathrm{\Psi }_\mathrm{I}^{}\mathrm{\Psi }_\mathrm{I}=|\psi _1|^2+|\psi _2|^2=2|\psi _1|^2,$$
(27)
which is positive-definite. The current vectors $`j_\pm `$ are respectively

$$j_+=\psi _2^{}\psi _1+\psi _1^{}\psi _2$$
(28)
and

$$j_{}=i\psi _2^{}\psi _1i\psi _1^{}\psi _2,$$
(29)
which take some constants. Therefore, the equation of continuity (19) holds and there is a conserved current. From the parameters (3) and the result of Eq. (28) and (29), we interpret that the evolution of the quantized universe is that whether one of the three spatial axes is expanding and two axes are contracting or one axis is contracting and two axes are expanding.

It should be noted that the quantity

$$\rho _\mathrm{I}(\alpha ,\beta _+,\beta _{})𝑑\beta _+𝑑\beta _{}=2|\psi _1|^2𝑑\beta _+𝑑\beta _{}$$
(30)
is not convergent and the wave function (24) is not the square integrable function. Therefore, it is impossible to normalize the wave function (24) without an appropriate boundary condition and we do not discuss the normalization constant in this paper. Here, the quantity (30) does not directly determine the probability, so that the probability must be regarded as a quantity proportional to this quantity (30). Moreover, since $`\rho /\alpha =0`$, the probability density is independent of the time-parameter and quantity (30) is conserved.

As mentioned above, it should be noted here that since the Wheeler-DeWitt equation and the Dirac type equation for the Bianchi type-II, VI0, VII0 and VIII models are reduced to that of the Bianchi type-I model at the singularity ($`\alpha \mathrm{}`$), all models have the same property that the behavior near the singularity becomes that of the plane wave solution. We shall discuss the quantized Hamiltonian system of the Bianchi type-I model in §4.

### 3.2 The Bianchi type-II model

The metric of the Bianchi type-II model is expressed as Eq. (2) and the 1-forms are given by

$$\chi ^1=dx^1x^3dx^2,\chi ^2=dx^2,\chi ^3=dx^3.$$
(31)
Substituting the 1-forms (31) into Eq. (4), we get the structure constants

$$C_{23}^1=C_{32}^1=1.$$
(32)
Other components are zero. The curvature is

$$^{\left(3\right)}R_{\mathrm{II}}=\frac{1}{2}e^{2\alpha +4\beta _++4\sqrt{3}\beta _{}},$$
(33)
and the anisotropic potential $`V_{\mathrm{II}}(\beta _+,\beta _{})`$ is obtained from the relation (6) as

$$V_{\mathrm{II}}1=\frac{1}{3}e^{4\beta _++4\sqrt{3}\beta _{}}.$$
(34)
Then, the Wheeler-DeWitt equation is expressed as

$$\left(\widehat{p}_\alpha ^2+\widehat{p}_+^2+\widehat{p}_{}^2+\frac{1}{3}e^{4\alpha +4\beta _++4\sqrt{3}\beta _{}}\right)\mathrm{\Psi }=0.$$
(35)
We postulate here that the factorized Hamiltonian is expressed as

$$_{\mathrm{II}}=\widehat{p}_\alpha +\sigma _1\widehat{p}_++\sigma _2\widehat{p}_{}+\frac{\sigma _3}{\sqrt{3}}e^{2\alpha +2\beta _++2\sqrt{3}\beta _{}}$$
(36)
and that the Dirac type equation holds:

$$_{\mathrm{II}}\mathrm{\Psi }_{\mathrm{II}}=0.$$
(37)
With the help of the relation (8), it becomes

$$i\frac{\mathrm{\Psi }_{\mathrm{II}}}{\alpha }=\left(\begin{array}{cc}\frac{1}{\sqrt{3}}e^{2\alpha +2\beta _++2\sqrt{3}\beta _{}}& i\frac{}{\beta _+}\frac{}{\beta _{}}\\ i\frac{}{\beta _+}+\frac{}{\beta _{}}& \frac{1}{\sqrt{3}}e^{2\alpha +2\beta _++2\sqrt{3}\beta _{}}\end{array}\right)\mathrm{\Psi }_{\mathrm{II}}.$$
(38)
Now we look for a solution to this equation. We take a trial function which gives the two-component spinor form and notices the form of the Dirac type equation as

$$\mathrm{\Psi }_{\mathrm{II}}=\mathrm{exp}\left[\frac{1}{6}e^{2\alpha +2\beta _++2\sqrt{3}\beta _{}}\right]\left(\begin{array}{c}\psi _1\\ \psi _2\end{array}\right),$$
(39)
where the spinor components $`\psi _1`$ and $`\psi _2`$ are arbitrary constants. Substituting the trail function (39) into Eq. (38), we get the relation of the spinor components

$$\psi _2=\psi _1.$$
(40)
Here we use the condition that the Dirac type equation (37) for the Bianchi type-II model is reduced to that for the Bianchi type-I model as the limit of the singularity ($`\alpha \mathrm{}`$). The trail function (39) would include such an asymptotic behavior. The trail function involving the plane wave part is expressed as

$$\mathrm{\Psi }_{\mathrm{II}}=e^{i(k_+\beta _++k_{}\beta _{}\omega \alpha )}\mathrm{exp}\left[\frac{1}{6}e^{2\alpha +2\beta _++2\sqrt{3}\beta _{}}\right]\left(\begin{array}{c}\psi _1\\ \psi _1\end{array}\right).$$
(41)
Although $`\omega `$, $`k_\pm `$ are arbitrary real constants, the coefficients $`\omega `$ and $`k_\pm `$ are not free at all. Substituting the solution (41) into Eq. (38), we have

$$\left(\begin{array}{cc}\omega & k_+ik_{}\\ k_++ik_{}& \omega \end{array}\right)\left(\begin{array}{c}\psi _1\\ \psi _1\end{array}\right)=0.$$
(42)
Therefore, the coefficients must satisfy

$$k_+=\omega ,k_{}=0.$$
(43)
Consequently, the solution is finally expressed in terms of the product of the amplitude and oscillation parts restricted by Eq. (43):

$$\mathrm{\Psi }_{\mathrm{II}}=e^{i\omega (\beta _++\alpha )}\mathrm{exp}\left[\frac{1}{6}e^{2\alpha +2\beta _++2\sqrt{3}\beta _{}}\right]\left(\begin{array}{c}\psi _1\\ \psi _1\end{array}\right).$$
(44)
The probability density is expressed, by using the wave function (44), as

$$\rho _{\mathrm{II}}(\alpha ,\beta _+,\beta _{})=\mathrm{\Psi }_{\mathrm{II}}^{}\mathrm{\Psi }_{\mathrm{II}}=2|\psi _1|^2\mathrm{exp}\left[\frac{1}{3}e^{2\alpha +2\beta _++2\sqrt{3}\beta _{}}\right].$$
(45)
This is always positive-definite. One sees that this probability density goes to zero as $`\beta _++\sqrt{3}\beta _{}\mathrm{}`$, while it approaches a certain constant as $`\beta _++\sqrt{3}\beta _{}\mathrm{}`$. The current vectors $`j_\pm `$ given by Eq. (18) are

$$j_+=2|\psi _1|^2\mathrm{exp}\left[\frac{1}{3}e^{2\alpha +2\beta _++2\sqrt{3}\beta _{}}\right]$$
(46)
and

$$j_{}=0.$$
(47)
Hence, the equation of continuity (19) holds and there is the conserved current. Moreover, it is shown that the wave part propagates to the direction of $`\beta _+`$ axis only and it agrees with the result of Eq. (43).

### 3.3 The Bianchi type-VI0 model

The metric of the Bianchi type-VI0 model is expressed as Eq. (2) and the 1-forms are given by

$$\chi ^1=\mathrm{cosh}x^3dx^1\mathrm{sinh}x^3dx^2,\chi ^2=\mathrm{sinh}x^3dx^1+\mathrm{cosh}x^3dx^2,\chi ^3=dx^3.$$
(48)
The structure constants are

$$C_{23}^1=C_{32}^1=1,C_{31}^2=C_{13}^2=1$$
(49)
and the rest are zero. The anisotropic potential becomes

$$V_{\mathrm{VI}_0}1=\frac{4}{3}e^{4\beta _+}\left(\mathrm{cosh}(4\sqrt{3}\beta _{})+1\right).$$
(50)
Then, the Wheeler-DeWitt equation for the Bianchi type-VI0 model is expressed as

$$\left[\widehat{p}_\alpha ^2+\widehat{p}_+^2+\widehat{p}_{}^2+\frac{4}{3}e^{4\beta _+}\left(\mathrm{cosh}(4\sqrt{3}\beta _{})+1\right)\right]\mathrm{\Psi }=0.$$
(51)
We postulate that the Dirac type equation for the Bianchi type-VI0 model is expressed as

$$_{\mathrm{VI}_0}\mathrm{\Psi }_{\mathrm{VI}_0}\left[\widehat{p}_\alpha +\sigma _1\widehat{p}_++\sigma _2\widehat{p}_{}+\frac{2}{\sqrt{3}}\sigma _3e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}+e^{2\sqrt{3}\beta _{}}\right)\right]\mathrm{\Psi }_{\mathrm{VI}_0}=0.$$
(52)
As in the case of the Bianchi type-II model, we assume a trial function of the form

$$\mathrm{\Psi }_{\mathrm{VI}_0}=\mathrm{exp}\left[\frac{1}{3}e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}e^{2\sqrt{3}\beta _{}}\right)\right]\left(\begin{array}{c}\psi _1\\ \psi _2\end{array}\right),$$
(53)
where the components $`\psi _1`$ and $`\psi _2`$ are arbitrary constants. Substituting the trial function (53) into the Dirac type equation (52), we get

$$\psi _2=\psi _1.$$
(54)
From the condition that the wave function (53) should be reduced to the plane wave solution at the limit of the singularity, we obtain

$$k_+=\omega ,k_{}=0.$$
(55)
The wave function involving the plane wave part is finally expressed as

$$\mathrm{\Psi }_{\mathrm{VI}_0}=e^{i\omega (\beta _++\alpha )}\mathrm{exp}\left[\frac{1}{3}e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}e^{2\sqrt{3}\beta _{}}\right)\right]\left(\begin{array}{c}\psi _1\\ \psi _1\end{array}\right).$$
(56)
The probability density becomes

$$\rho _{\mathrm{VI}_0}(\alpha ,\beta _+,\beta _{})=\mathrm{\Psi }_{\mathrm{VI}_0}^{}\mathrm{\Psi }_{\mathrm{VI}_0}=2|\psi _1|^2\mathrm{exp}\left[\frac{2}{3}e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}e^{2\sqrt{3}\beta _{}}\right)\right],$$
(57)
which is positive-definite. The current vectors are

$$j_+=2|\psi _1|^2\mathrm{exp}\left[\frac{2}{3}e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}e^{2\sqrt{3}\beta _{}}\right)\right],j_{}=0.$$
(58)
One sees that the equation of continuity holds and there is a conserved current. The plane wave part of the wave function (56) propagates to the direction of $`\beta _+`$ axis only. It is to be noted that since the wave function (56) and the probability density (57) are both divergent as $`\beta _{}\mathrm{}`$, this solution is unphysical in such a region.

### 3.4 The Bianchi type-VII0 model

The metric of the Bianchi type-VII0 model is expressed as Eq. (2) and the 1-forms are given by

$$\chi ^1=\mathrm{cos}x^1dx^2+\mathrm{sin}x^1dx^3,\chi ^2=\mathrm{sin}x^1dx^2+\mathrm{cos}x^1dx^3,\chi ^3=dx^1.$$
(59)
The structure constants are

$$C_{23}^1=C_{32}^1=1,C_{31}^2=C_{13}^2=1$$
(60)
and the rest are zero. The anisotropic potential is

$$V_{\mathrm{VII}_0}1=\frac{4}{3}e^{4\beta _+}\left(\mathrm{cosh}(4\sqrt{3}\beta _{})1\right).$$
(61)
Then, the Wheeler-DeWitt equation for the Bianchi type-VII0 model is expressed as

$$\left[\widehat{p}_\alpha ^2+\widehat{p}_+^2+\widehat{p}_{}^2+\frac{4}{3}e^{4\beta _+}\left(\mathrm{cosh}(4\sqrt{3}\beta _{})1\right)\right]\mathrm{\Psi }=0.$$
(62)
Here, we postulate that the Dirac type equation with a factorized Hamiltonian is expressed as

$$_{\mathrm{VI}_0}\mathrm{\Psi }_{\mathrm{VII}_0}\left[\widehat{p}_\alpha +\sigma _1\widehat{p}_++\sigma _2\widehat{p}_{}+\frac{2}{\sqrt{3}}\sigma _3e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}e^{2\sqrt{3}\beta _{}}\right)\right]\mathrm{\Psi }_{\mathrm{VII}_0}=0.$$
(63)
As in the case of previous models, we can find out the wave function

$$\mathrm{\Psi }_{\mathrm{VII}_0}=e^{i\omega (\beta _++\alpha )}\mathrm{exp}\left[\frac{1}{3}e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}+e^{2\sqrt{3}\beta _{}}\right)\right]\left(\begin{array}{c}\psi _1\\ \psi _1\end{array}\right).$$
(64)
The probability density becomes

$$\rho _{\mathrm{VII}_0}(\alpha ,\beta _+,\beta _{})=\mathrm{\Psi }_{\mathrm{VII}_0}^{}\mathrm{\Psi }_{\mathrm{VII}_0}=2|\psi _1|^2\mathrm{exp}\left[\frac{2}{3}e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}+e^{2\sqrt{3}\beta _{}}\right)\right],$$
(65)
which is positive-definite. The current vectors have

$$j_+=2|\psi _1|^2\mathrm{exp}\left[\frac{2}{3}e^{2\alpha +2\beta _+}\left(e^{2\sqrt{3}\beta _{}}+e^{2\sqrt{3}\beta _{}}\right)\right],j_{}=0.$$
(66)
One sees that the equation of continuity holds and there is a conserved current. The plane wave part of the wave function (64) propagates to the direction of $`\beta _+`$ axis only. At the limit of $`\beta _+\mathrm{}`$, the value of the probability density approaches a certain constant, while converges to zero as $`\beta _+\mathrm{}`$ or $`\beta _{}\pm \mathrm{}`$.

### 3.5 The Bianchi type-VIII model

The metric of the Bianchi type-VIII model is expressed as Eq. (2) and the 1-forms are given by

$$\{\begin{array}{c}\chi ^1=\mathrm{cosh}x^2\mathrm{cos}x^3dx^1\mathrm{sin}x^3dx^2,\hfill \\ \chi ^2=\mathrm{cosh}x^2\mathrm{sin}x^3dx^1+\mathrm{cos}x^3dx^2,\hfill \\ \chi ^3=\mathrm{sinh}x^2dx^1+dx^3.\hfill \end{array}$$
(67)
The structure constants are

$$C_{23}^1=C_{32}^1=1,C_{31}^2=C_{13}^2=1,C_{12}^3=C_{21}^3=1$$
(68)
and the rest are zero. The anisotropic potential becomes

$$V_{\mathrm{VIII}}1=\frac{1}{3}\left(e^{8\beta _+}2e^{4\beta _+}+4e^{2\beta _+}\mathrm{cosh}(2\sqrt{3}\beta _{})+2e^{4\beta _+}\mathrm{cosh}(4\sqrt{3}\beta _{})\right).$$
(69)
However, since the anisotropic potential (69) is a very complicated expression, we cannot present the Dirac type equation of the form as Eq. (11).

## 4 The behavior near the singularity

It is interesting for us to discuss the behavior of the wave function near the singularity. Misner discussed the Bianchi type- IX (Mixmaster) model which is represented by the classical dynamical system where an imaginary particle moves in the anisotropic potential of the approximated regular triangle and regarded the dynamical system as the billiard in the $`\beta _+`$-$`\beta _{}`$ plane. In this way Misner’s picture is very helpful to discuss the universe through the anisotropic parameters.

The vacuum Bianchi type-I universe is so-called Kasner Universe. The evolution of the model in classical region presents an anisotropically expanding universe in which two axises are expanding and the other is contracting. While, in quantum region the quadratic Hamiltonian is derived by the ADM approach and the evolution of the quantized universe is presented by the Wheeler-DeWitt equation like Klein-Gordon type. There is the problem that the probability density can become negative by the Klein-Gordon approach.

Let us consider to change the point of view to our formalism. Near the singularity ($`\alpha \mathrm{}`$), the Dirac type equations for the Bianchi Class A models are reduced to that of Bianchi type-I

$$i\frac{\mathrm{\Psi }}{\alpha }=\left(\sigma _1\widehat{p}_++\sigma _2\widehat{p}_{}\right)\mathrm{\Psi }\widehat{h}^{}\mathrm{\Psi }.$$
(70)
This dynamical system describes a massless imaginary particle moving in zero potential. The equations of motion in the Heisenberg representation for $`\widehat{\beta }_+`$ and $`\widehat{\beta }_{}`$ are

$$\frac{\widehat{\beta }_+}{\alpha }=i[\widehat{h}^{},\widehat{\beta }_+]=i\sigma _1[\widehat{p}_+,\widehat{\beta }_+]=\sigma _1$$
(71)
and

$$\frac{\widehat{\beta }_{}}{\alpha }=i[\widehat{h}^{},\widehat{\beta }_{}]=i\sigma _2[\widehat{p}_{},\widehat{\beta }_{}]=\sigma _2,$$
(72)
where the commutative relation $`[\widehat{\beta }_b,\widehat{p}_a]=i\delta _{ab}(a,b=+,)`$ is used. It is similar to the Zitterbewegung  (a rapidly oscillating motion) of electron, so that the imaginary particle moves a trembling motions near the singularity of the universe. This behavior in quantum mechanics with factorized Hamiltonian is very different from that for the Wheeler-DeWitt equation. We emphasize here that this kind of behavior cannot be obtained from the Wheeler-DeWitt equation.

Here, we set an orbital angular momentum $`\widehat{m}^{}`$ as

$$\widehat{m}^{}=\widehat{\beta }_+\widehat{p}_{}\widehat{\beta }_{}\widehat{p}_+.$$
(73)
The Heisenberg equation of motion for the orbital angular momentum is

$$\frac{d\widehat{m}^{}}{d\alpha }=i[\widehat{h}^{},\widehat{m}^{}]=\sigma _1\widehat{p}_{}\sigma _2\widehat{p}_+.$$
(74)
We have, furthermore,

$$\frac{d\sigma _3}{d\alpha }=i[\widehat{h}^{},\sigma _3]=2\sigma _2\widehat{p}_+2\sigma _1\widehat{p}_{}.$$
(75)
The quantity $`\widehat{m}^{}+\frac{1}{2}\sigma _3`$ is a constant of motion. This result one can interpret by saying a massless imaginary particle has a ‘spin-like’ angular momentum $`\frac{1}{2}\sigma _3`$. The Dirac type equation (70) is a similar type to the Weyl equation, which has two-dimension of the spatial part here. We know that the Weyl equation describes the behavior of the massless neutrino with a half spin. On the analogy of this characteristic with the Weyl equation, the massless imaginary particle would have a half spin. Accordingly, we think that the Dirac type equation for the Bianchi type-I model describes the universe as a massless imaginary particle with ‘spin-like’ degree of freedom. Thus, the ‘spin-like’ degrees of freedom cause the universe to go through an early quantum stage of evolution with agitated anisotropy-oscillation. However, the physical interpretation of the ‘spin-like’ degree of freedom is not comprehensible yet. In conclusion, we obtain that the Dirac type equation for the Bianchi type-I model shows the interesting results on quantum behavior, unlike Misner’s negative conclusion that the status of the wave function remains classical near the initial singularity .

## 5 Conclusions and Discussion

We have investigated the first-order equation using the Pauli matrices by factorizing the Wheeler-DeWitt equation for the vacuum Bianchi Class A models except for the Bianchi type-IX. We have derived the solution expressed in terms of two-component spinor on the Bianchi type-I, II, VI0 and VII0 models. It has been shown that the probability density becomes positive-definite and the equation of continuity holds with a conserved current. However, we could not find the solution to the Bianchi type-VIII model because the form of anisotropic potential is very complicated expression. It is to be noted that we could not normalize the wave function and the probability density for each model. In the Bianchi type-II model, it is able to factorize the quadratic Hamiltonian of the Wheeler-DeWitt equation into self-adjoint Hamiltonians other than the Hamiltonian (36). However, we cannot extract physically interesting solutions to the Dirac type equations by such Hamiltonians. The similar results are obtained for the Bianchi type-VI0 and type-VII0 models.

The Dirac type equations for the Bianchi type-II, VI0 and VII0 model are reduced to that for the Bianchi type-I model near the singularity ($`\alpha \mathrm{}`$). The quantized universe near the singularity is expressed by the plane wave solution with two-component spinor and the anisotropic parameters $`\widehat{\beta }_+`$ and $`\widehat{\beta }_{}`$ make the motion of Zitterbewegung. The Zitterbewegung can be related to the ‘spin-like’ degree of freedom and thus the ‘spin-like’ degree of freedom causes the universe to go necessarily through an early quantum stage of evolution with agitated anisotropy-oscillation. The ‘spin-like’ degree of freedom which does not appear when we deal with the Wheeler-DeWitt equation for the Bianchi models may play an important role when we discuss the anisotropy of the universe in the initial stage. In such a sense, the Zitterbewegung is a significant subject how the big-bang or the inflation is related. To solve the problem of the initial singularity we must find out the physical interpretation of the ‘spin-like’ degree of freedom newly found by the Dirac factorization.

As mentioned in §1, our method cannot be applied to the Bianchi type-IX model. Because, the anisotropic potential of the Bianchi type-IX $`V_{\mathrm{IX}}`$ is expressed in terms of

$$V_{\mathrm{IX}}=\frac{1}{3}e^{8\beta _+}\frac{4}{3}e^{2\beta _+}\mathrm{cosh}2\sqrt{3}\beta _{}+\frac{2}{3}e^{4\beta _+}(\mathrm{cosh}4\sqrt{3}\beta _{}1)+1>0.$$
(76)
The quantity $`v^2=V_{\mathrm{IX}}1`$ takes negative value as well as positive value here. It means that the quantity $`v=\sqrt{V_{\mathrm{IX}}1}`$ which appears in the factorized theory turns out to be imaginary as well as real. Therefore, we cannot factorize the quadratic Hamiltonian into a fixed form. Hence, we must not use the Pauli matrices but the Dirac matrices. To remedy this defect, we notice the fact that the anisotropic potential $`V_{\mathrm{IX}}`$ itself is non-negative, and carry out the Dirac procedure of factorization with anticommuting $`\gamma ^\mu `$ ($`\mu =0,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2},\mathrm{\hspace{0.17em}3}`$) and $`\gamma _5`$ matrices. One possibility is to make the factorization in the form

$$\left(\gamma ^0p_\alpha +\gamma ^1p_++\gamma ^2p_{}+\gamma ^3v^{}e^{2\alpha }+\gamma _5e^{2\alpha }\right)\mathrm{\Psi }=0,$$
(77)
where we set $`v^{}=\sqrt{V_{\mathrm{IX}}}`$. However, as we have seen, because the anisotropic potential $`V_{\mathrm{IX}}`$ is a very complicated expression, we could not have found out the solution to the Dirac type equation (77) with the factorized Hamiltonian.

We compare our formalism with the Dirac-Square-Root (DSR) formalism . On the DSR, it is postulated that the factorized Dirac type equation for the Bianchi type-I model is expressed as

$$\left(i\frac{}{\alpha }i\sigma _1\frac{}{\beta _+}i\sigma _2\frac{}{\beta _{}}\right)\mathrm{\Psi }(\alpha ,\beta _+,\beta _{})=0.$$
(78)
The Dirac type equation (78) then implies that the Wheeler-DeWitt equation for the Bianchi type-I model yields

$$\left(\frac{^2}{\alpha ^2}+\frac{^2}{\beta _+^2}+\frac{^2}{\beta _{}^2}\right)\mathrm{\Psi }(\alpha ,\beta _+,\beta _{})=0,$$
(79)
for each component of the wave function. According to Mallett , it is postulated that the Dirac type equation is expressed as

$$\left(i\frac{}{\alpha }i\sigma _1\frac{}{\beta _+}i\sigma _2\frac{}{\beta _{}}+iW(\alpha )\right)\mathrm{\Psi }(\alpha ,\beta _+,\beta _{})=0,$$
(80)
where $`W(\alpha )`$ is an unknown function depended on $`\alpha `$ only. In Kim and Oh  paper, it has been applied the first-order equation (80) to the Bianchi type-IX model.

However, it is not valid to apply the DSR formalism to the Wheeler-DeWitt equation including the anisotropic potential term like (76) and to factorize in the form of the Dirac type equation (80). Because, it is inconsistent with the original Wheeler-DeWitt equation in the point that the differential operators with the anisotropic parameters $`\beta _+`$ and $`\beta _{}`$ do not act on the unknown function $`W(\alpha )`$. The incorrect equation yields if one assume the equation (80) when the Wheeler-DeWitt equation has the anisotropic potential term included the anisotropic parameters. While, our formalism is valid to apply to the Bianchi Class A models because it is able to factorize the Wheeler-DeWitt equation included the anisotropic potential term into the Dirac type equation (11). This is the extended DSR formalism.

When we yield the original expression from the Dirac type equation discussed by our formalism and the iterating equation, we also obtain the extra terms which include the time parameter and the anisotropic parameters. Regarding the extra terms as a new constraint to the wave function and solving the constraint, we cannot get solutions which satisfy the Dirac type equation. As a consequence, we only discuss the Dirac type equation without imposing the constraint condition on the wave function.

## Acknowledgments

The author is grateful to Profs. Tetsuya Hara and Ikuo S. Sogami for the motivation of this work and helpful discussion.


format: text

text:
# Multicomponent mixture of charged hard-sphere chain molecules in the polymer mean-spherical approximation

## I Introduction

In recent years substantial amout of efforts has been focused on the development of theoretical approaches predicting the equilibrium structure and thermodynamical properties of charged polymer solutions. This information enables one to esteblish the stability limits and to build the phase diagrams of such fundamentally and technologically importent systems as polyelectrolyte solutions (for example DNA, polyacrilic acid) and polyampholyte solutions (for example proteins). One of the simplest continuum models of charged polymer solutions is represented by the fluid of charged hard-sphere linear chain flexible molecules. In spite of its simplicity this model incorporates several essential features of charged polymer fluids, such as excluded volume effects, chain connectivity and flexibility, distribution of the charge along the chain backbone and Coulomb interaction. Recently several off-lattice theories for the fluid of flexible charged chain molecules, based on the integral-equation techniques developed for the fluids of small molecules, have been proposed. These include extension  of the polymer reference interaction site model (RISM) integral equation theory of Schweizer and Curro , extension  of the simple interpolation scheme (SIS) of Stell and Zhou  and extensions  of the multidensity integral equation theory for associating fluids of Wertheim . The latter studies are based on the analytical solution of the mean spherical approximation (MSA) version of the multidensity theory, the so-called polymer MSA (PMSA) . In the case of uncharged system direct application of Wertheim’s original formulation by Chang and Sandler  gave rise to a polydisperse fluid of chains with a prescribe mean number of beads. A subsequent extension by Kalyuzhnyi and Cummings  yields a model that polymerizes to a fluid of chains of fixed length in the complete association limit. In that limit, two-density Wertheim’s Ornstein-Zernike (OZ) equation becomes identical to the Rossky-Chiles  version of the ”proper” site-site integral equations  (”proper-RISM”) as first noted by Kalyuzhnyi and Stell  and later discussed in more details by Stell  and by Kalyuzhnyi and Cummings . More recently complete association version of this extension was combined with the ideal chain approximation  and generalized for the systems with arbitrary long-range potential outside the hard-core region . As a result the corresponding version of the PMSA was formulated and H$`ø`$ye-Stell scheme  of calculating the simple fluids thermodynamical properties was extended in the frames of the PMSA . Subsequently von Solms and Chiew  combine PMSA and ideal chain approximation in a way similar to that proposed by Chang and Sandler . The resulting polydisperse chain mixture was used to describe thermodynamical properties of the two-component mixture of fixed-length chain molecules with equally charged beads and oppositelly charged monomer counterions. Results of the theory appeare to be in reasonable agreement with computer simulation results . We note in passing, that ideal chain approximation utilized by von Solms and Chiew  is somewhat different from the original one . In the corresponding Wertheim’s OZ equation the authors take into account all partial correlation functions, while due the original ideal chain approximation  correlations involving at least one doubly bonded particle are neglected. More recently general solution of the PMSA  for polymerizing charged hard spheres, supplemented by the ideal chain approximation, was used to derive closed form analytical expressions for thermodynamical properties of charged chain fluid with a certain restrictions imposed on the sizes of the chain beads . In the latter study the authors propose an improved version of the ideal chain approximation, which satisfies the Debye-Hückel limiting law.

In this paper we present an analytical solution of the ideal chain PMSA  for the multicomponent mixture of charged hard-sphere flexible chain molecules. Unlike the previous studies proposed version of the PMSA is quite general and apply for any mixture of the chain molecules with arbitrary distribution of the charge and size of the beads along the molecular backbone. Our solution represents a complete association limit of the general solution obtained earlier  for the fluid of polymerizing charged hard spheres. Similar, as in the case of primitive electrolyte model  this solution reduces to the solution of only one algebraic equation for the MSA-like screening parameter $`\mathrm{\Gamma }`$. We also generalize recently developed method of calculating PMSA thermodynamics via the energy route  for the multicomponent case and derived closed form analytical expressions for internal energy, Helmholtz free energy, chemical potentials and pressure in terms of the present solution. By way of illustration we consider several different versions of the fluid of charged hard-sphere chain molecules, including the molecules with uniform, diblock and alternating distribution of the charge along the chain backbone.

The paper is orgenized as follows. In the next section we discuss detailes of the model to be studied and formulate ideal chain PMSA theory. Section III gives solution of the ideal chain PMSA and Section IV contains expressions for thermodynamic and structure properties of the system. Numerical results and discussion can be found in Section V and in Section VI we collect our conclusions.

## II The model and ideal chain PMSA theory

The model fluid we consider consists of the $`M`$-component mixture of freely jointed tangent hard-sphere chain molecules with each molecule of species $`a`$ represented by $`m_a`$ charged hard-sphere sites of sizes $`\sigma _\alpha ^a`$ and charges $`ez_\alpha ^a`$. The total number density of the system is $`\rho =_a\rho _a`$, where $`\rho _a`$ is the number density of species $`a`$, and we assume charged neutrality conditions

$$\underset{a\alpha }{}\rho _az_\alpha ^a=0$$
(1)
We denote the molecular species by the small letters $`a,b,c,\mathrm{}`$ taking the values $`1,2,\mathrm{},M`$ and the site type in a given molecule $`a`$ by a small greek letters $`\alpha ,\beta ,\gamma ,\mathrm{}`$, which take the values $`1,2,\mathrm{},m_a`$. Thus the site has two indices, one denoting the species of the molecule and the other the site type in the molecule. For example the hard-sphere diameter of the site of type $`\alpha `$ belonging to the molecule of the species $`a`$ is denoted by $`\sigma _\alpha ^a`$.

The site-site pair potential $`U_{\alpha \beta }^{ab}(r)`$ between the sites of the type $`\alpha `$ and $`\beta `$ belonging to the molecules of species $`a`$ and $`b`$ can be written in the form

$$U_{\alpha \beta }^{ab}(r)=U_{\alpha \beta }^{(hs)ab}(r)+U_{\alpha \beta }^{(C)ab}(r)$$
(2)
where $`U_{\alpha \beta }^{(hs)ab}(r)`$ is the hard-sphere potential and $`U_{\alpha \beta }^{(C)ab}(r)`$ is the Coulomb potential

$$U_{\alpha \beta }^{(C)ab}(r)=\frac{e^2z_\alpha ^az_\beta ^b}{ϵ_0r}$$
(3)
with $`ϵ_0`$ being the dielectric constant of the continuum.

In this study we are using the multidensity version of the mean spherical approximation (MSA)  (or polymer MSA (PMSA)) supplemented by the so-called ideal chain approximation . Both PMSA and the ideal chain approximation have been discussed at length earlier and therefore we will omit any details and present here only the final expressions. The theory consists of the Ornstein-Zernike-like integral equation

$$\widehat{𝐡}_{\alpha \beta }^{ab}(k)=\widehat{𝐜}_{\alpha \beta }^{ab}(k)+\underset{c}{}\rho _c\underset{\gamma }{}\widehat{𝐜}_{\alpha \gamma }^{ac}(k)𝜶\widehat{𝐡}_{\gamma \beta }^{cb}(k)$$
(4)
and PMSA boundary conditions

$$\{\begin{array}{ccccc}𝐜_{\alpha \beta }^{ab}(r)\hfill & =\hfill & 𝐄\beta U_{\alpha \beta }^{(C)ab}(r)+\frac{𝐭_{\alpha \beta }^{ab}}{2\pi \sigma _{\alpha \beta }}\delta (r\sigma _{\alpha \beta }^{ab}),\hfill & r\hfill & >\sigma _{\alpha \beta }^{ab}=\frac{1}{2}\left(\sigma _\alpha ^a+\sigma _\beta ^b\right)\hfill \\ 𝐡_{\alpha \beta }^{ab}(r)\hfill & =\hfill & 𝐄,\hfill & r\hfill & <\sigma _{\alpha \beta }^{ab}\hfill \end{array}$$
(5)
where $`\widehat{𝐡}_{\alpha \beta }^{ab}(k)`$ and $`\widehat{𝐜}_{\alpha \beta }^{ab}(k)`$ are the matrices with the elements being the Fourier transforms of the elements of the matrices $`𝐡_{\alpha \beta }^{ab}(r)`$ and $`𝐜_{\alpha \beta }^{ab}(r)`$

$$𝐡_{\alpha \beta }^{ab}(r),𝐜_{\alpha \beta }^{ab}(r)=\left(\begin{array}{ccc}c_{\alpha _0\beta _0}^{ab}(r)& c_{\alpha _0\beta _A}^{ab}(r)& c_{\alpha _0\beta _B}^{ab}(r)\\ c_{\alpha _A\beta _0}^{ab}(r)& c_{\alpha _A\beta _A}^{ab}(r)& c_{\alpha _A\beta _B}^{ab}(r)\\ c_{\alpha _B\beta _0}^{ab}(r)& c_{\alpha _B\beta _A}^{ab}(r)& c_{\alpha _B\beta _B}^{ab}(r)\end{array}\right),$$

$`𝐭_{ij}`$, $`𝜶`$ and $`𝐄`$ are the following matrices

$$t_{\alpha _i\beta _j}^{ab}=\frac{\delta _{ab}}{2\rho _a}\left[\delta _{iA}\delta _{jB}\frac{\delta _{\alpha \beta +1}}{\sigma _{\alpha \alpha 1}}+\delta _{iB}\delta _{jA}\frac{\delta _{\alpha \beta 1}}{\sigma _{\alpha \alpha +1}}\right],\alpha _{ij}=1\delta _{ij}+\delta _{0i}\delta _{0j},E_{ij}=\delta _{0i}\delta _{0j}.$$

Here the lower indices $`i,j`$ each take the values $`0,A,`$ and $`B`$, and denote the bonding states of the corresponding particle, i.e. $`0`$ denotes the unbonded state and $`A`$ and $`B`$ denote the $`A`$-bonded and $`B`$-bonded states, respectively. The total partial correlation functions $`h_{\alpha _i\beta _j}^{ab}(r)`$ are related to the site-site total correlation functions $`h_{\alpha \beta }^{ab}(r)`$ by the following relation

$$h_{\alpha \beta }^{ab}(r)=\underset{ij}{}h_{\alpha _i\beta _j}(r)$$
(6)
The set of the OZ equations (4) together with the PMSA closure conditions (5) represent our ideal chain PMSA theory for the fluid of charged linear chain hard-sphere molecules.

## III Solution of the ideal chain PMSA for the fluid of charged chain molecules

Recently a general solution of the PMSA for the multicomponent mixture of polymerizing charged hard spheres was published . Here we elaborate on this general solution by utilizing additionally the ideal chain approximation  and specializing it for the case of the model at hand.

The general solution was obtained using Baxter’s tecnique, which factorizes the initial OZ equation (4) into a set of two equations

$$𝐒_{\alpha \beta }^{ab}(|r|)=𝐐_{\alpha \beta }^{ab}(r)\underset{c}{}\rho _c\underset{\gamma }{}𝑑r^{}𝐐_{\alpha \gamma }^{ac}(r^{})𝜶\left[𝐐_{\beta \gamma }^{bc}(r^{}r)\right]^T$$
(7)
$$𝐉_{\alpha \beta }^{ab}(|r|)=𝐐_{\alpha \beta }^{ab}(r)+\underset{c}{}\rho _c\underset{\gamma }{}𝑑r^{}𝐉_{\alpha \gamma }^{ac}(|r^{}r|)𝜶𝐐_{\gamma \beta }^{cb}(r^{})$$
(8)
where $`T`$ denotes the transpose matrix. The projections $`𝐒_{\alpha \beta }^{ab}(r)`$ and $`𝐉_{\alpha \beta }^{ab}(r)`$

$$𝐒_{\alpha \beta }^{ab}(r)=2\pi _r^{\mathrm{}}𝑑r^{}r^{}𝐜_{\alpha \beta }^{ab}(r^{}),𝐉_{\alpha \beta }^{ab}(r)=2\pi _r^{\mathrm{}}𝑑r^{}r^{}𝐡_{\alpha \beta }^{ab}(r^{})$$
(9)
satisfy the following boundary conditions

$$\{\begin{array}{ccccc}𝐉_{\alpha \beta }^{ab}(r)\hfill & =\hfill & \pi r^2𝐄+𝐉_{\alpha \beta }^{ab},\hfill & r\hfill & \sigma _{\alpha \beta }^{ab}\hfill \\ 𝐒_{\alpha \beta }^{ab}(r)\hfill & =\hfill & \frac{\beta e^2}{ϵ_0}\frac{z_\alpha ^az_\beta ^be^{\mu |r|}}{\mu }𝐄,\hfill & r\hfill & >\sigma _{\alpha \beta }^{ab}\hfill \end{array}$$
(10)
obtained from (5). Here $`𝐉_{\alpha \beta }^{ab}=𝐉_{\alpha \beta }^{ab}(0)`$ and the limit $`\mu 0`$ is to be taken at the end of the calculations.

From the analysis of the equations (7) (8) we get the following expression for the factor function $`𝐐_{\alpha \beta }^{ab}(r)`$

$$𝐐_{\alpha \beta }^{ab}(r)=\theta \left(\sigma _{\alpha \beta }^{ab}r\right)\left[𝐪_{\alpha \beta }^{ab}(r)+𝐭_{\alpha \beta }^{ab}\right]\left(\stackrel{~}{𝐳}_\alpha ^a\right)^T\stackrel{~}{𝐚}_\beta ^b,r>\lambda _{\alpha \beta }^{ab}$$
(11)
where $`\lambda _{\alpha \beta }^{ab}=\left(\sigma _\alpha ^a\sigma _\beta ^b\right)/2`$ and the function $`𝐪_{\alpha \beta }^{ab}(r)`$ is defined in the range $`\lambda _{\alpha \beta }^{ab}<r<\sigma _{\alpha \beta }^{ab}`$ by

$$𝐪_{\alpha \beta }^{ab}(r)=\frac{1}{2}\stackrel{~}{𝐄}^T\stackrel{~}{𝐀}_\beta ^b\left(r\sigma _{\alpha \beta }^{ab}\right)\left(r\lambda _{\beta \alpha }^{ba}\right)+\beta _{\alpha \beta }^{ab}\left(r\sigma _{\alpha \beta }^{ab}\right).$$
(12)
Here $`\stackrel{~}{𝐄},\stackrel{~}{𝐚}_\alpha ^a`$ and $`\stackrel{~}{𝐳}_\alpha ^a`$ are the row vectors

$$\stackrel{~}{𝐄}=(1,0,0),\stackrel{~}{𝐚}_\alpha ^a=(a_{\alpha _0}^a,a_{\alpha _A}^a,a_{\alpha _B}^a),\stackrel{~}{𝐳}_\alpha ^a=(z_\alpha ^a,0,0).$$

Coefficients of the Baxter $`q`$-function can be expressed in terms of only one parameter, i.e. the screening PMSA parameter $`\mathrm{\Gamma }`$ . We have

$$\beta _{\alpha _i\beta _j}^{ab}=\delta _{i0}\delta _{j0}\frac{\pi \sigma _\beta ^b}{\mathrm{\Delta }}+\frac{2\pi \beta ^{}}{\sigma _\alpha ^a\mathrm{\Gamma }}\left(X_{\alpha _0}^a\delta _{i0}z_\alpha ^a\right)X_{\beta _j}^b,$$
(13)
$$A_{\beta _j}^b=\frac{4\pi \beta ^{}}{\mathrm{\Gamma }}\eta ^BX_{\beta _j}^b+\delta _{j0}\frac{2\pi }{\mathrm{\Delta }}\left[1+\zeta _2\sigma _\beta ^b\frac{\pi }{2\mathrm{\Delta }}\right]\frac{\pi }{\mathrm{\Delta }}\left[\left(1\delta _{\beta m_a}\right)\frac{\delta _{jB}\sigma _{\beta +1}^b}{\sigma _{\beta +1\beta }^{bb}}+\left(1\delta _{\beta 1}\right)\frac{\delta _{jA}\sigma _{\beta 1}^b}{\sigma _{\beta 1\beta }^{bb}}\right]$$
(14)
$$\mathrm{\Delta }=1\frac{1}{6}\pi \zeta _3,\zeta _n=\underset{d}{}\rho _d\underset{k}{}\left(\sigma _k^d\right)^n,\chi _2=\underset{d}{}\rho _d\underset{k}{}z_k^d\left(\sigma _k^d\right)^2$$
(15)
$$a_{\alpha _i}^a=2\pi \beta ^{}X_{\alpha _i}^a/\mathrm{\Gamma },$$
(16)
$$X_{\alpha _0}^a=\left[z_\alpha ^a\eta ^B\left(\sigma _\alpha ^a\right)^2\right]\mathrm{\Gamma }_\alpha ^a,X_{\alpha _i}^a=\sigma _\alpha ^a\left[\tau _{\alpha _i}^a(z)\eta ^B\tau _{\alpha _i}^a(\sigma ^2)\right],(i0)$$
(17)
$$\eta ^B=\frac{\frac{\pi }{2\mathrm{\Delta }}\underset{d}{}\rho _d\underset{\gamma }{}\sigma _\gamma ^d\left\{z_\gamma ^d\mathrm{\Gamma }_\gamma ^d+\sigma _\gamma ^d\left[\tau _{\gamma _A}^d(z)+\tau _{\gamma _B}^d(z)\right]\right\}}{1+\frac{\pi }{2\mathrm{\Delta }}_d\rho _d_\gamma (\sigma _\gamma ^d)^2\left[\sigma _\gamma ^d\mathrm{\Gamma }_\gamma ^d+\tau _{\gamma _A}^d(\sigma ^2)+\tau _{\gamma _B}^d(\sigma ^2)\right]},\mathrm{\Gamma }_\alpha ^a=\left(1+\sigma _\alpha ^a\mathrm{\Gamma }\right)^1$$
(18)
$$\tau _{\alpha _A}^a(y)=\frac{1}{2}\underset{\gamma =2}{\overset{\alpha }{}}\mathrm{\Gamma }_{\gamma 1}^a\rho _a^{\alpha \gamma }\frac{y_{\gamma 1}^a}{\sigma _{\gamma \gamma 1}^{aa}}\underset{\delta =\gamma }{\overset{\alpha }{}}\mathrm{\Gamma }_\delta ^a\left[2^{\gamma \alpha }\left(1\delta _{\gamma \alpha }\right)\underset{\omega =\gamma +1}{\overset{\alpha }{}}\frac{\sigma _{\omega 1}^a}{\sigma _{\omega \omega 1}^a}+\delta _{\gamma \alpha }\right]\left(1\delta _{\alpha 1}\right)$$
(19)
$$\tau _{\alpha _B}^a(y)=\frac{1}{2}\underset{\gamma =\alpha }{\overset{m_a1}{}}\mathrm{\Gamma }_{\gamma +1}^a\rho _a^{\gamma \alpha }\frac{y_{\gamma +1}^a}{\sigma _{\gamma \gamma +1}^{aa}}\underset{\delta =\alpha }{\overset{\gamma }{}}\mathrm{\Gamma }_\delta ^a\left[2^{\alpha \gamma }\left(1\delta _{\gamma \alpha }\right)\underset{\omega =\alpha }{\overset{\gamma 1}{}}\frac{\sigma _{\omega +1}^a}{\sigma _{\omega \omega +1}^a}+\delta _{\gamma \alpha }\right]\left(1\delta _{\alpha m_a}\right)$$
(20)
Here $`\mathrm{\Gamma }`$ satisfies the following equation

$$\mathrm{\Gamma }^2=\pi \beta ^{}\underset{a\alpha }{}\rho _a\underset{\alpha }{}\stackrel{~}{𝐗}_\alpha ^a𝜶\left(\stackrel{~}{𝐗}_\alpha ^a\right)^T$$
(21)
where $`\stackrel{~}{𝐗}_\alpha ^a`$ is the row vector $`\stackrel{~}{𝐗}_\alpha ^a=(X_{\alpha _0}^a,X_{\alpha _A},X_{\alpha _B})`$.

Thus solution of the PMSA for the present model reduces to the solution of the algebraic equation (21) for the screening parameter $`\mathrm{\Gamma }`$.

## IV Structure and thermodynamics

### A Structure properties

The contact values of the regular part of the partial pair distribution functions $`g_{\alpha _i\beta _j}^{ab}(\sigma _{\alpha \beta }^{ab}+)`$ follow from the relation (8) after differentiating it with respect to $`r`$ and taking the limit of $`r\sigma _{\alpha \beta }^{ab}+`$

$$2\pi \sigma _{\alpha \beta }^{ab}g_{\alpha _i\beta _j}^{ab}(\sigma _{\alpha \beta }^{ab}+)=2\pi \delta _{0i}\delta _{j0}\sigma _{\alpha \beta }^{ab}g_{\alpha \beta }^{(hs)ab)}(\sigma _{\alpha \beta }^{ab}+)2\pi \beta ^{}X_{\alpha _i}^aX_{\beta _j}^b+T_{\alpha _i\beta _j}^{ab}$$
(22)
where $`g_{\alpha \beta }^{(hs)ab}(\sigma _{\alpha \beta }^{ab}+)`$ is the hard-sphere contact values and

$$T_{\alpha _i\beta _j}^{ab}=\delta _{ab}\{[\frac{(1\delta _{\alpha m_a})(1\delta _{\alpha m_a1})\delta _{iB}\delta _{jA}\delta _{\alpha \beta 2}}{4\rho _a\sigma _{\alpha \alpha +1}^{aa}\sigma _{\alpha +1\beta }^{aa}}+\frac{(1\delta _{\alpha 1})(1\delta _{\alpha 2})\delta _{iA}\delta _{jB}\delta _{\alpha \beta +2}}{4\rho _a\sigma _{\alpha \alpha 1}^{aa}\sigma _{\alpha 1\beta }^{aa}}]+$$

$$\frac{\pi }{2\mathrm{\Delta }}\delta _{j0}\sigma _\beta ^a\left[(1\delta _{\alpha m_a})\delta _{iB}\frac{\sigma _{\alpha +1}^a}{\sigma _{\alpha \alpha +1}^{aa}}+(1\delta _{\alpha 1})\delta _{iA}\frac{\sigma _{\alpha 1}^a}{\sigma _{\alpha \alpha 1}^{aa}}\right]+$$

$$\frac{\pi }{2\mathrm{\Delta }}\delta _{i0}\sigma _\alpha ^a[(1\delta _{\beta m_a})\delta _{jB}\frac{\sigma _{\beta +1}^a}{\sigma _{\beta +1\beta }^{aa}}+(1\delta _{\beta 1})\delta _{jA}\frac{\sigma _{\beta 1}^a}{\sigma _{\beta 1\beta }^{aa}}]\}$$
(23)
The values of $`g_{\alpha _i\beta _j}^{ab}(r)`$ for $`r>\sigma _{\alpha \beta }^{ab}`$ can be calculated using the following relation

$$g_{\alpha _i\beta _j}^{ab}(r)=\delta _{0i}\delta _{0j}+\frac{1}{2\pi ^2r}_0^{\mathrm{}}\left[\widehat{\gamma }_{\alpha _i\beta _j}^{ab}(k)\delta _{i0}\delta _{j0}\beta \widehat{U}_{\alpha \beta }^{(C)ab}(k)\right]k\mathrm{sin}(kr)𝑑k$$
(24)
where $`\widehat{U}_{\alpha \beta }^{(C)ab}(k)`$ is the Fourier transform of the Coulomb potential $`U_{\alpha \beta }^{(C)ab}(r)`$ and expression for the function $`\widehat{\gamma }_{\alpha _i\beta _j}^{ab}(k)=\widehat{h}_{\alpha _i\beta _j}^{ab}(k)\widehat{c}_{\alpha _i\beta _j}^{ab}(k)`$ follows from the set of equations (7) and (8), written in the Fourier $`k`$-space

$$\{\begin{array}{ccc}𝝆^1\widehat{𝐜}(k)\hfill & =\hfill & \widehat{𝐐}(k)𝝆\widehat{𝐐}^T(k)\hfill \\ 𝝆^1+\widehat{𝐡}(k)\hfill & =\hfill & \left[𝝆\widehat{𝐐}(k)𝝆\widehat{𝐐}^T(k)𝝆\right]^1\hfill \end{array},$$
(25)
which gives

$$\widehat{𝜸}(k)=\widehat{𝐐}(k)𝝆\widehat{𝐐}^T(k)+\left[𝝆\widehat{𝐐}(k)𝝆\widehat{𝐐}^T(k)𝝆\right]^12𝝆^1.$$
(26)
Here $`𝝆`$, $`\widehat{𝐡}(k)`$, $`\widehat{𝐜}(k)`$ and $`bQ(k)`$ are the matrices with the elements $`\rho _{\alpha _i\beta _j}^{ab}=\delta _{ab}\delta _{\alpha \beta }\alpha _{ij}\rho _a`$, $`\widehat{h}_{\alpha _i\beta _j}^{ab}(k)`$, $`\widehat{c}_{\alpha _i\beta _j}^{ab}(k)`$ and

$$\widehat{Q}_{\alpha _i\beta _j}^{ab}(k)=\left[𝝆^1\right]_{\alpha _i\beta _j}^{ab}_{\lambda _{\beta \alpha }}^{\mathrm{}}𝑑re^{ikr}Q_{\alpha _i\beta _j}^{ab}(r),$$
(27)
respectively. Substituting expression for the fuction $`Q_{\alpha _i\beta _j}^{ab}(r)`$ (11) into the right hand side of the Fourier transformation (27) we find that

$$\widehat{𝐐}_{\alpha \beta }^{ab}(k)=\delta _{ab}\delta _{\alpha \beta }\left(\rho _a𝜶\right)^1\{\stackrel{~}{𝐄}\stackrel{~}{𝐀}_\beta ^b\phi _2(k,\sigma _\alpha ^a)+(𝜷_{\alpha \beta }^{ab}\frac{1}{2}\stackrel{~}{𝐄}\stackrel{~}{𝐀}_\beta ^b\sigma _\alpha ^a)\phi _1(k,\sigma _\alpha ^a)+$$

$$+(𝜷_{\alpha \beta }^{ab}\sigma _\alpha ^a𝐭_{\alpha \beta }^{ab})\phi _0(k,\sigma _\alpha ^a)(\stackrel{~}{𝐳}_\alpha ^a)^T\stackrel{~}{𝐚}_\beta ^b\frac{i}{k}\}$$
(28)
where

$$\phi _2(k,\sigma )=\frac{1}{k^3}\left[\left(\frac{1}{2}k^2\sigma ^2+ik\sigma 1\right)e^{ik\sigma }+1\right],\phi _1=\frac{1}{k^2}\left[\left(1ik\sigma \right)e^{ik\sigma }1\right],$$

$$\phi _0(k,\sigma )=\frac{i}{k}\left(e^{ik\sigma }1\right)$$
(29)
This result allow us to derive the expression for the structure factor

$$\widehat{S}_{\alpha \beta }^{ab}(k)=\delta _{ab}\delta _{\alpha \beta }+\sqrt{\rho _a\rho _b}\widehat{h}_{\alpha \beta }^{ab}(k)$$
(30)
Using the second of the relations (25) we have

$$\widehat{S}_{\alpha \beta }^{ab}(k)=\frac{1}{\sqrt{\rho _a\rho _b}}\left\{\left[\widehat{𝐐}(k)𝝆\widehat{𝐐}^T(k)\right]^1\right\}_{\alpha _0\beta _0}^{ab}$$
(31)
### B Thermodynamic properties

Expression for the excess internal energy $`\mathrm{\Delta }E`$ follows from the standard relation

$$\frac{\mathrm{\Delta }E}{V}=2\pi \underset{ab}{}\rho _a\rho _b\underset{\alpha \beta }{}_0^{\mathrm{}}𝑑rr^2U_{\alpha \beta }^{(C)ab}(r)g_{\alpha \beta }^{ab}(r),$$
(32)
which after some calculations gives

$$\beta \frac{\mathrm{\Delta }E}{V}=\beta ^{}\underset{a}{}\rho _a\underset{\alpha }{}z_\alpha ^aN_\alpha ^a$$
(33)
where $`\sigma _\alpha ^aN_\alpha ^a=_{i=0}^BX_{\alpha _i}^az_\alpha ^a`$

The other thermodynamical properties of the system in question can be calculated following the method developed earlier . This method generalizes H$`ø`$ye-Stell MSA energy route to thermodynamics  in the frames of the present PMSA approach. Originally this generalization was developed for the one-component case. Extension of the method in the multicomponent case is rather straightforward and yields the following expressions for the Helmholtz free energy $`A`$, pressure $`P`$ and chemical potential $`\mu _\alpha ^a`$ in excess to their reference system values

$$\beta \frac{AA^{(ref)}}{V}=J+\frac{1}{2}\underset{ab}{}\rho _a\rho _b\underset{\alpha \beta }{}\left\{𝜶\left[\widehat{𝐜}_{\alpha \beta }^{ab}\widehat{𝐜}_{\alpha \beta }^{(0)ab}\right]𝜶\right\}_{00}\beta \frac{\mathrm{\Delta }E}{V}+$$

$$+\frac{1}{3}\pi \underset{ab}{}\rho _a\rho _b\underset{\alpha \beta }{}\left(\sigma _{\alpha \beta }^{ab}\right)^3Tr\left[𝐠_{\alpha \beta }^{ab}𝜶𝐠_{\beta \alpha }^{ba}𝜶𝐠_{\alpha \beta }^{(0)ab}𝜶𝐠_{\beta \alpha }^{(0)ba}𝜶\right]$$

$$\frac{1}{3}\underset{ab}{}\rho _a\rho _b\underset{\alpha \beta }{}\left(\sigma _{\alpha \beta }^{ab}\right)^2Tr\left\{𝐭_{\alpha \beta }^{ab}𝜶\left[𝐠_{\alpha \beta }^{ab}𝐠_{\alpha \beta }^{(0)ab}\right]𝜶\right\},$$
(34)
$$\beta \left(PP^{(ref)}\right)=\beta \frac{AA^{(ref)}}{V}\frac{1}{2}\underset{ab}{}\rho _a\rho _b\underset{\alpha \beta }{}\left\{𝜶\left[\widehat{𝐜}\widehat{𝐜}^{(0)}\right]𝜶\right\}_{00}+\beta \frac{\mathrm{\Delta }}{V},$$
(35)
$$\beta \rho _a\left(\mu _\alpha ^a\mu _\alpha ^{(ref)a}\right)=\frac{1}{2}\rho _a\underset{b}{}\rho _b\underset{\beta }{}\left\{𝜶\left[\widehat{𝐜}_{\alpha \beta }^{ab}\widehat{𝐜}_{\alpha \beta }^{(0)ab}\right]𝜶\right\}_{00}\beta \frac{\mathrm{\Delta }E_\alpha ^a}{V}$$
(36)
where the quantities with the superscript $`(ref)`$ denote the reference system quantities and the quantities with superscript $`(0)`$ denote the corresponding zero charge PMSA quantities

$$\widehat{𝐜}=\widehat{𝐜}(k=0),\widehat{𝐜}_{\alpha \beta }^{ab}=\widehat{𝐜}_{\alpha \beta }^{ab}(k=0),𝐠_{\alpha \beta }^{ab}=𝐠_{\alpha \beta }^{ab}(\sigma _{\alpha \beta }^{ab}+),𝐠_{\alpha \beta }^{ab}=\frac{𝐠_{\alpha \beta }^{ab}(r=\sigma _{\alpha \beta }^{ab}+)}{r},$$

$$\frac{\mathrm{\Delta }E_\alpha ^a}{V}=2\pi \rho _a\underset{b}{}\rho _b\underset{\beta }{}_0^{\mathrm{}}r^2g_{\alpha \beta }^{ab}(r)U_{\alpha \beta }^{(C)ab}(r)𝑑r$$
(37)
and

$$J=\frac{1}{6}\beta \underset{ab}{}\rho _a\rho _b\underset{\alpha \beta }{}_{r>\sigma _{\alpha \beta }^{ab}+}g_{\alpha \beta }^{ab}(r)𝐫\mathbf{}U_{\alpha \beta }^{(C)ab}(r)𝑑𝐫$$
(38)
The reference system is represented by the multicomponent mixture of uncharged hard-sphere chain molecules. Expressions (34), (35) and (36) are quite general and apply for any type of the potential outside the hard-core. The quantities, which enter these expressions can be written in terms of the present solution of PMSA for the model at hand. We have

$$\beta \frac{\mathrm{\Delta }E_\alpha ^a}{V}=\beta ^{}\rho _az_\alpha ^aN_\alpha ^a,J=\frac{1}{3}\frac{\mathrm{\Delta }E}{V},$$
(39)
$$\sigma _{\alpha \beta }^{ab}𝐠_{\alpha \beta }^{ab}=𝐡_{\alpha \beta }^{ab}+\underset{d}{}\rho _d\underset{\gamma }{}\sigma _{\alpha \gamma }^{ad}𝐡_{\alpha \gamma }^{ad}𝜶\left[𝐭_{\gamma \beta }^{db}(\stackrel{~}{𝐳}_\gamma ^d)^T\stackrel{~}{𝐚}_\beta ^b\right]𝐄𝜶\underset{d}{}\rho _d\underset{\alpha }{}\lambda _{\gamma \alpha }^{da}𝐭_{\gamma \beta }^{db}$$

$$\frac{1}{2}\underset{d}{}\rho _d\underset{\gamma }{}\sigma _\gamma ^d\left\{\left[\frac{1}{\pi }𝐭_{\alpha \gamma }^{ad}𝜶\frac{1}{6}(\sigma _\gamma ^d)^2\mathrm{𝟏}\right]\stackrel{~}{𝐄}^T\stackrel{~}{𝐀}_\beta ^b+\left[2\sigma _{\alpha \gamma }^{ad}𝐠_{\alpha \gamma }^{ad}𝐄\sigma _\gamma ^d\frac{1}{\pi \sigma _\gamma ^d}𝐭_{\alpha \gamma }^{ad}\right]𝜶𝜷_{\gamma \beta }^{db}\right\},$$
(40)
$$\underset{ab}{}\rho _a\rho _b\underset{\alpha \beta }{}\left[𝜶\widehat{𝐜}_{\alpha \beta }^{ab}𝜶\right]_{00}=\underset{ab}{}\rho _a\rho _b\underset{\alpha \beta }{}\left\{𝜶\left[(𝐌_{\beta \alpha }^{ba})^T+𝐌_{\alpha \beta }^{ab}+\underset{c}{}\rho _c\underset{\gamma }{}𝐌_{\alpha \gamma }^{ac}𝜶(𝐌_{\beta \gamma }^{bc})^T\right]𝜶\right\}_{00},$$
(41)
where

$$𝐌_{\alpha \beta }^{ab}=\frac{1}{12}(\sigma _\alpha ^a)^3\stackrel{~}{𝐄}^T\stackrel{~}{𝐀}_\beta ^b\frac{1}{2}(\sigma _\alpha ^a)^2𝜷_{\alpha \beta }^{ab}+\sigma _\alpha ^a𝐭_{\alpha \beta }^{ab}+\lambda _{\beta \alpha }^{ba}(\stackrel{~}{𝐳}_\alpha ^a)^T\stackrel{~}{𝐚}_\beta ^b.$$
(42)
Here expression (41) is obtained by differentiating of the equation (8) twice and considering the limit of $`r\sigma _{\alpha \beta }^{ab}+`$, and expression (42) follows from the first of the equations (25) at $`k0`$.

## V Results and discussion

To illustrate the above solution of the PMSA in this section we present the numerical results for thermodynamical properties of the several versions of charged chain fluid model. We consider two-component mixture of chain polyions with the beads of equal charge $`z_\alpha ^p=z_p=1`$ and oppositely charged counterions $`z_\alpha ^c=z_c=1`$ (model M1) and two versions of three-component mixture of chain polyions with diblock and alternating distribution of oppositelly charged beads, $`z_+^p=z_{}^p=1`$, and two types of counterions with opposite charges $`z_+^c=z_{}^c=1`$ (models M2 and M3). The densities of the counterions of these two models were choosed to be $`\rho _+^c=\rho _{}^c=\frac{1}{2}m_p\rho _p`$. Here the indices $`p`$ and $`c`$ denote polyion and counterion, respectively. We also consider the one-component version of charged chain fluid models with diblock and alternating distribution of the charge along the chain backbone. In all cases the hard-sphere diameters of the chain beads $`\sigma _\alpha ^p`$ and counterions $`\sigma _\alpha ^c`$ were choosed to be equal, $`\sigma _\alpha ^p=\sigma _\alpha ^c=\sigma =1`$. Schematic representation of the models studied are shown in figure 1. Thermodynamical properties of the uncharged versions of these models were calculated using TPTD theory of Chang and Sandler .

In figure 2 we compare osmotic coefficient results obtained using the present version of the PMSA, von Solms and Chiew version of PMSA  and molecular dynamic (MD) computer simulation method  for the model M1 at a different chain length ($`m_p=16,32`$ and $`64`$) and Bjerrum length $`\lambda _B=\beta ^{}=0.833`$. In the region of the intermediate densities our theory underestimate and the theory of von Solms and Chiew overestimate the value of the osmotic coefficient $`\varphi =\beta P/\rho `$, where $`\rho =_a\rho _a`$. In the region of the low densities both theories predict for the osmotic coefficient its ideal gas value of $`\varphi =0`$, while MD simulation gives the value, which is slightly above 1. This disagreement could be due to the uncertainties of the MD simulation in the diluted region. In general predictions of the present approach are more accurate than those of von Solms and Chiew approach .

In figures 3-5 we compare the density dependence of the osmotic pressure and osmotic coefficient at different values of Bjerrum length ($`\lambda _B=0.833,2.499`$) for different polyelectrolyte models. In the case of low Bjerrum length ($`\lambda _B=0.833`$) one can see the linear dependence of the osmotic pressure with respect to the density on the $`loglog`$ scale (figures 3a-5a). In this region the osmotic pressure is almost independent of the chain length and distribution of the charge along the chain. With the increase of $`\lambda _B`$ ($`\lambda _B=2.499`$) there is a departure from the linearity, which is substantial for the model M1, slightly smaller for the model M2 and almost negligible for the model M3. At the same time one can see the substantial chain length dependence of the osmotic pressure in the case of the model M2. This dependence becomes smaller in the case of the model M1 and is absent for the model M3 (figures 3b-5b). Thus for the systems with more random distribution of the charge along the chain backbone the chain length dependence is smaller. With the increase of the chain length this dependence becomes weaker and it is reasonable to expect, that for $`m_p>64`$ it will be negligible for all the models studied.

In figur 6 the liquid-gas phase diagrams for the one-component version of the models with diblock and alternating distribution of the charge are presented. We consider the systems with molecular chain length $`m_p=2,4,6,8,10`$. Computer simulation results are available only in the case of the chains length $`m_p=2`$ . In general theoretical results are in qualitative agreement with computer simulation results with certain disagreement for the slope steepness of the liquid branches of the coexistence curves; the slope of theoretical curve is steeper than that of the computer simulation curve. In addition, the position of the theoretical critical point is shifted towards lower values of the density and higher values of the temperature. According to the previous discussion the systems with diblock distribution of the charge have larger degree of nonideality in comparison with the systems with alternating distribution of the charge. This is reflected in the figure 6, where the critical temperature of the systems with diblock chains is higher than that of the systems with alternating chains of the same length. The corresponding values of the critical densities is larger in the case of systems with alternating chains. In both cases increase of the chain length leads to the increase of the critical temperature and slight decrease (increase) of the critical density for diblock (alternating) distribution of the charge. We note in passing, that due to the lack of the solution of the equation for screening parameter $`\mathrm{\Gamma }`$ (21) at lower densities we were not able to determine the coexistance curve for the systems with alternating chain charge of the chain length $`m_p=8,10`$ below the temperatures shown in figure 6b.

Finally in figure 7 we present the liquid-gas phase diagram for the system of trimers with the middle bead charge $`z_2^p=2z_1^p=2z_3^p`$, where $`z_1^p`$ and $`z_3^p`$ are the charges of the terminal beads, $`z_1^p=z_3^p=1`$. Similar as in the case of alternating charged chain system with $`m_p=8,9`$ there is no solution of the equation (21) in the low density region, therefore we were not able to calculate the coexistence curve for the temperatures lower than those shown in the figure. This model can be seen as a complete association limit of the primitive 2:1 electrolyte model. Recent computer simulation studies  of the electrolyte primitive model suggest that in the vicinity of the coexistance region the system is highly associated and fraction of charged clusters, including free ions, is negligible. Therefore the equilibrium properties of the system are determined by the presence of the neutral ionic clusters, to which the lower order clusters give the main contribution. Assuming that in the case of the electrolyte restricted primitive model (RPM) all such clusters are represented by the neutral ionic pairs the RPM liquid-gas phase diagram has been successfully reproduced by the Monte Carlo (MC) computer simulation method  and theoretically . Making similar assumption in the case of the 2:1 primitive electrolyte model the corresponding liquid-gas phase diagram can be modeled using the phase diagram obtained for trimers. Unfortunately computer simulation predictions for the entire phase diagram of the 2:1 electrolyte model currently is not available. Recently Camp and Patey  present MC estimate for the critical temperature, which appeares to be $`T_c^{}=0.1\pm 0.01`$. The rough estimate of the ranges for the critical density can be obtained from the figure 5 of ; our estimate is $`\rho _c^{}=0.115\pm 0.085`$. This estimate for the critical point together with Debye-Hückel (DH) theory estimate  and MSA phase diagram are shown in figure 7. In addition, for the sake of completeness, we present the liquid-gas phase diagram for the RPM of electrolyte, obtained from the MC simulation method , MSA and our PMSA for neutral diatomics . One can see that predictions of the PMSA is much better in comparison with usual MSA. In both cases DH predictions for the critical temperature is close to those of the PMSA, however the corresponding predictions for the critical density are much worse.

## VI Concluding remarks

In this study we obtain an analytical solution of the ideal chain PMSA for the multicomponent mixture of charged hard-sphere flexible linear chain molecules. This solution apply to any mixture of chain molecules with arbitrary distribution of the charge and size of the beads along the molecular backbone. We present closed form analytical expressions for thermodynamical and structure properties of the system. These expressions are used to calculate thermodynamics of several different versions of charged hard-sphere chain model. We consider two-component mixture of chains with equally charged beads and oppositely charged counterions, two versions of the three-component mixture of chains with diblock and alternating distribution of the charge and two types of oppositely charged counterions. Comparison of the theoretical predictions for the osmotic pressure of the two-component model with corresponding computer simulation predictions shows that the present version of the ideal chain PMSA theory is reasonably accurate and gives somewhat better agreement than that of the theory of von Solms and Chiew . Effects due to the difference in Bjerrum length, molecular chain length and distribution of the charge was also studied. The largest deviation from ideal behaviour was demonstrated by the system with uniform distribution of the charge; behaviour of the system with alternating distribution of the charge shows much less nonideality. This conclusion is varyfied by the comparison of the liquid-gas phase diagrams calculated for the one-component version of the models with diblock and alternating distribution of the charge; the latter model has substatially lower critical temperature than that of the former. With the increase of the chain length the critical temperature for both models increases, while the corresponding critical density slightly increases for alternating charge model and decreases for diblock charge model. Finally it was demonstrated, that liquid-gas phase diagram of the neutral trimers with doubly charged middle bead and singly charged terminal beads of the opposite sighn can be used to correct MSA results for the phase diagram of 2:1 primitive electrolyte model.

ACKNOWLEDGMENTS

CAPTIONS TO THE FIGURES

Figure 1. Schematic diagrams of the models with uniform distribution of the charge (M1), diblock distribution of the charge (M2) and alternating distribution of the charge.

Figure 2. Osmotic coefficient $`\varphi =PV/NkT`$ as a function of the packing fraction $`\eta =\pi \rho \sigma ^3/6`$ for the model M1 at $`\lambda _B=0.833`$ for $`m_p=16`$ (a), $`m_p=32`$ (b) and $`m_p=64`$ (c). Predictions of the present theory (solid lines), von Solms and Chiew theory  (dashed lines) and computer simulation predictions  (diamons).

Figure 3. Osmotic pressure $`P^{}=\beta P\sigma ^3`$ (a) and osmotic coefficient $`\varphi =PV/NkT`$ (b) as a function of packing fraction $`\eta =\pi \rho _T\sigma ^3/6`$ ($`\rho _T=_am_a\rho _a`$ for the model M1 from the present theory at $`\lambda _B=2.499`$ (solid lines), $`\lambda _B=0.833`$ (dashed lines) and $`\lambda _B=0`$ (dashed-dotted lines). For each set of lines from the top to the bottom at $`\mathrm{log}\eta =2`$ $`m_p=8,16,32,64`$. Symbols are computer simulation predictions  for $`\lambda _B=0.833`$ and for $`m_p=16`$ (diamonds), $`m_p=32`$ (circles) and $`m_p=64`$ (squares).

Figure 4. The same as in figure 3 for the model M2.

Figure 5. The same as in figure 3 for the model M3.

Figure 6. Liquid-gas phase diagram in $`T^{}=1/\beta ^{}`$ versus reduced density of the beads $`\rho _b^{}=m_p\rho \sigma ^3`$ coordinates for one component chain fluids with diblock distribution of the charge (a) and with alternating distribution of the charge (b). From the top to the bottom at $`\rho _b^{}=0.05`$ $`m_p=10,8,6,4,2`$. Solid lines are predictions from the present theory and diamonds are computer simulation predictions  for $`m_p=2`$.

Figure 7 Liquid-gas phase diagram in $`T^{}=1/\beta ^{}`$ versus reduced ionic density $`\rho ^{}=(\rho _++\rho _{})\sigma ^3`$ coordinates for restricted primitive model (RPM) of electrolyte (lower portions of the figure for $`T^{}>0.08`$) and for the primitive model (PM) of electrolyte with 2:1 assymetry in charge (upper portions of the figure for $`T^{}<0.08`$). Present theory (solid lines), MSA (dashed lines), Debye-Hückel theory  (squares), computer simulation estimate for the critical point of 1:2 electrolyte PM  (open diamonds) and for the phase diagram of electrolyte RPM (open circles , solid diamonds ).