text:
# Effective interactions and phase behaviour for a model clay suspension in an electrolyte *footnote **footnote *This paper is dedicated to Jean-Pierre Hansen on the occasion of his 60th birthday

## I Introduction

With the possibility to form an orientational ordering, the phase behaviour of anisotropic colloids is richer than its counterpart for spherical particles. Surprisingly, the isotropic-nematic transition expected on purely entropic grounds (excluded volume effects) has been extensively reported for rod-like colloids in the last sixty years , but could only be observed recently for discotic particles in apolar media . This experimental work directly points to the subtle effect of electrostatic interactions, given that in the widely studied model system of colloidal platelets, namely aqueous clay dispersions, the isotropic/nematic phase separation is hindered by a ubiquitous “fluid-solid” transition.

In spite of an important experimental and theoretical effort in the last 10 years, triggered by the emergence of Laponite as a model system for disc-like colloidal suspensions, the above transition is far from being well understood . In particular, the precise nature of the phases observed experimentally is debated , so that the relevance of the traditional terminology “sol-gel” to describe the transition is somewhat controversial.

On the other hand, theoretical approaches to describe the phase behaviour of charged disc-like particles in an electrolyte are in their infancy , essentially because there is no satisfactory generalization for discs of the isotropic DLVO potential  describing the effective interactions between a pair of spherical colloids. In the simplest situation of two coaxial and parallel uniformly charged platelets, the effective Coulomb force has been computed within linear  and non-linear  Poisson-Boltzmann (PB) theory. A screened electrostatic pair potential has been worked out , allowing to compute analytically the interaction energy for discs of arbitrary orientation; this approach holds at the level of linearized PB theory (weak electrostatic coupling), and is expressed as a perturbative expansion in the parameter $`\kappa r_0`$, $`r_0`$ being the disc radius and $`\kappa `$ the inverse Debye length in the electrolyte. It consequently becomes less reliable when $`\kappa r_0`$ becomes of order 1 or larger (the typical situation for Laponite solutions), or when the charge on the platelets becomes too large (high electrostatic coupling, also typical of clay systems). The latter limitation may be circumvented by introducing the concept of charge renormalization , while the former requires the re-summation of all the powers of $`\kappa r_0`$ involved in the “multi-polar” expansion propounded in . In this work, we discuss how to take both aspects into account and investigate the consequences on the phase behaviour. Our goal is to understand the effect of colloid anisotropy on a possible phase transition. We show that considering electrostatic effects only and discarding van der Waals interactions yields qualitative differences between spherical and discotic colloids and that the equivalent of the fluid-solid transition for spheres is nevertheless reminiscent of the “sol-gel” transition reported in .

The paper is organized as follows. In section II, we obtain a pair potential for discs valid at large distances within linearized PB theory. We show that at asymptotically large separations, this potential remains strongly anisotropic, unlike its counterpart in vacuum or in a plain dielectric medium (no salt). The approach is generalized in section III to an arbitrary “one-” or“two-dimensional” colloid (i.e. with vanishing internal volume). The non-linear phenomenon of counter-ion “condensation” leading to a renormalization of the bare platelet charge is then considered in section IV and the resulting effective charge successfully tested against “exact” numerical simulations reported in the literature. In section V, we propose a first and simplified investigation into the phase behaviour by mapping the interaction energy onto an isotropic Yukawa potential. Finally, the full anisotropic potential including charge renormalization is used to investigate the mesostructure of an assembly of interacting discs by means of Monte Carlo computer simulations (section VI), and conclusions are drawn in section VII.

## II Large distance screened electrostatic potential

We consider first a unique platelet $`𝒫`$ of radius $`r_0`$ and surface charge $`\sigma `$ in a 1:1 electrolyte of bulk density $`n`$ (infinite dilution limit). The solvent is assumed to be a dielectric continuum of permittivity $`\epsilon `$. Within linearized Poisson-Boltzmann theory (LPB), the dimensionless electrostatic potential obeys the following Poisson’s equation

$$^2\varphi =\kappa ^2\varphi ,$$
(1)
where $`\kappa ^2=8\pi \mathrm{}_Bn`$ is the inverse Debye length squared and $`\mathrm{}_B=e^2/(kT\epsilon )`$ denotes the Bjerrum length corresponding to the distance where the bare Coulomb potential felt by two elementary point charges becomes equal to the thermal energy $`kT`$ ($`\mathrm{}_B7`$Å in water at room temperature). The potential is chosen to vanish at infinity and $`e`$ denotes the elementary charge.

We are interested in the behaviour of $`\varphi `$ at large distances, which is relevant to describe the interactions in dilute suspensions. The solution of Eq. (1) has been worked out in the form of an integral representation in , or as a multi-polar-like expansion . It is however more convenient to obtain an explicit expression, that is a good starting point to derive the pair interaction. This can be achieved by writing the solution as a convolution of the surface charge of the platelets with the screened Coulomb potential:

$$\varphi (𝐫)=\mathrm{}_B_𝒫\frac{\sigma }{e}\frac{\mathrm{exp}(\kappa |𝐫𝐬|)}{|𝐫𝐬|}d^2𝐬.$$
(2)
It is important to realize that such a convolution procedure would give incorrect results for polyions of non vanishing excluded volume (e.g. spheres or rods of non zero radii, see the appendix). Introducing the unit vector $`\widehat{𝐧}`$ normal to the disc surface, and the angle $`\theta [0,\pi /2]`$ between the corresponding direction and the position vector $`𝐫`$ with origin at platelet centre, we obtain the leading contribution

$$\varphi (𝐫)\stackrel{\kappa r1}{}\mathrm{}_BZ_{\text{bare}}\mathrm{\hspace{0.17em}2}\frac{\text{I}_1(\kappa r_0\mathrm{sin}\theta )}{\kappa r_0\mathrm{sin}\theta }\frac{e^{\kappa r}}{r}\text{with}Z_{\text{bare}}e=\pi r_0^2\sigma ,$$
(3)
where $`\text{I}_1`$ denotes the modified Bessel function of the first kind, such that $`\text{I}_1(x)x/2`$ for $`|x|1`$. A peculiarity of screened Coulomb potential appears at this point: whereas the large $`r`$ potential becomes isotropic and behaves as $`Z_{\text{bare}}/r`$ in vacuum (situation corresponding to the limit $`\kappa 0`$), the anisotropy is present at all distances in an electrolyte ($`\kappa 0`$): the $`\theta `$ and $`r`$ dependence factorize in equation (3). In practice, the anisotropy of the potential is generically significant when the size of the object under consideration is larger than the Debye length (i.e. $`\kappa r_0>1`$ here). For Laponite discs of radius 150 Å this crossover corresponds to an ionic strength $`I^{}=10^4`$M. As noted in , $`I^{}`$ corresponds experimentally to a threshold value delimiting qualitatively different phase behaviours. At fixed distance $`r`$, the potential (3) is minimum for $`\theta =0`$ which corresponds to the configuration where the “average” distance between the point where $`\varphi `$ is computed and the platelet is maximal; on the other hand, the potential is maximum for $`\theta =\pi /2`$. From (3) the ratio between these extremal values is $`r`$-independent, and reads $`2\text{I}_1(\kappa r_0)/(\kappa r_0)`$ which can be as large as 10 for $`\kappa r_0=5`$ (important anisotropy).

This result may be used to compute the potential energy of interactions between two platelets $`A`$ and $`B`$ with arbitrary relative orientations, as shown in Fig. 1. This energy is obtained by integrating the screened potential created by disc $`A`$ defined by Eq. (3) over the surface charge distribution of the second disc ($`B`$):

$$V_{AB}(r,\theta _A,\theta _B)=_B\frac{\sigma }{e}\varphi _A(𝐫+𝐬)d^2𝐬,$$
(4)
where $`V_{AB}`$ is the dimensionless potential expressed in $`kT`$ units. Again, such a procedure would not give the proper interactions in the case of polyions with internal volume. At large distances, we obtain

$$V_{AB}(r,\theta _A,\theta _B)Z_{\text{bare}}^2\mathrm{}_B\mathrm{\hspace{0.17em}4}\frac{\text{I}_1(\kappa r_0\mathrm{sin}\theta _A)}{\kappa r_0\mathrm{sin}\theta _A}\frac{\text{I}_1(\kappa r_0\mathrm{sin}\theta _B)}{\kappa r_0\mathrm{sin}\theta _B}\frac{e^{\kappa r}}{r}.$$
(5)
Clearly, the relative orientation of the platelets is not completely specified by the three parameters $`r`$, $`\theta _A`$ and $`\theta _B`$, but the omitted Euler angles only appear in higher order terms such as $`\mathrm{exp}(\kappa r)/r^2`$, $`\mathrm{exp}(\kappa r)/r^3`$…In particular the energy (5) appears insensitive to a precession of the discs around their centre to centre direction.

At fixed centre to centre distance, the above energy (always repulsive ) is maximized for co-planar discs ($`\theta _A=\theta _B=\pi /2`$, which corresponds to the maximum overlap of electric double layers) and minimized when the discs are co-axial and parallel ($`\theta _A=\theta _B=0`$, see the configuration represented in Fig. 2). A situation of intermediate electrostatic energy is that of T-shape perpendicular discs ($`\theta _A=0`$ and $`\theta _B=\pi /2`$). These results are consistent with the numerical linearized Poisson Boltzmann pair potential reported in , where it was also shown that the interactions between platelets at constant surface charge or constant potential were qualitatively very similar.

In reference , the potential energy $`V_{AB}`$ was obtained in the form of a multi-polar expansion. Unlike its unscreened counterpart where the multipole of order $`l`$ has a $`1/r^{l+1}`$ large distance contribution to the potential, this expansion is such that the multipole of order $`l`$ contributes to orders $`\mathrm{exp}(\kappa r)/r,\mathrm{exp}(\kappa r)/r^2\mathrm{}\mathrm{exp}(\kappa r)/r^{l+1}`$. The large $`r`$ potential can in principle be obtained by re-summation of all multi-polar contributions, which should lead to expression (5). Alternatively, truncating the multi-polar expansion at a given order amounts to expanding (5) in powers of $`\kappa r_0`$. Including monopole-monopole, monopole-quadrupole and quadrupole-quadrupole interactions, the following expression was obtained in :

$$V_{AB}(r,\theta _A,\theta _B)=Z_{\text{bare}}^2\mathrm{}_B\left[1+\frac{\kappa ^2r_0^2}{8}(\mathrm{sin}^2\theta _A+\mathrm{sin}^2\theta _B)+\frac{\kappa ^4r_0^4}{64}(\mathrm{sin}^2\theta _A\mathrm{sin}^2\theta _B)\right]\frac{e^{\kappa r}}{r}.$$
(6)
From (5), we obtain

$$V_{AB}(r,\theta _A,\theta _B)=Z_{\text{bare}}^2\mathrm{}_B\left[1+\frac{\kappa ^2r_0^2}{8}(\mathrm{sin}^2\theta _A+\mathrm{sin}^2\theta _B)+\frac{\kappa ^4r_0^4}{64}(\mathrm{sin}^2\theta _A\mathrm{sin}^2\theta _B+\frac{1}{3}\mathrm{sin}^4\theta _A\mathrm{sin}^4\theta _B)+𝒪(\kappa ^6r_0^6)\right]\frac{e^{\kappa r}}{r}$$
(7)
Both expressions agree at order $`(\kappa r_0)^2`$, and the difference at order $`(\kappa r_0)^4`$ is the monopole-hexadecapole contribution which has not been included in (6). This comparison illustrates the perturbative nature of the potential derived in , and the fact that the corresponding multipolar contributions are implicitly re-summed in expression (5).

## III Generalization to a polyion of arbitrary shape

The method used in the previous section may be generalized to find the far field electrostatic potential of an arbitrary polyion with vanishing internal volume and bare charge $`Z=_{\text{polyion}}\sigma (𝐬)d^2𝐬`$. In general, the convolution solution (2) admits the large distance behaviour:

$$\varphi (𝐫)Z\mathrm{}_Bf(\widehat{𝐫},\kappa )\frac{e^{\kappa r}}{r}+𝒪\left(\frac{e^{\kappa r}}{r^2}\right)\text{for}\kappa r1.$$
(8)
In this expression, the anisotropic part of the potential again factorizes from the $`r`$-dependence and is given by

$$f(\widehat{𝐫},\kappa )=_{\text{polyion}}\frac{\sigma (𝐬)}{Ze}\mathrm{exp}\left(\kappa \widehat{𝐫}𝐬\right)d^2𝐬,$$
(9)
where $`\widehat{𝐫}`$ denotes a unit vector in the direction of the position $`𝐫`$ where the potential is computed, and $`\sigma (𝐬)`$ is the surface charge density at point $`𝐬`$ on the colloid. As expected, the isotropic bare Coulomb potential is recovered in the limit of low electrolyte density:

$$\underset{\kappa 0}{lim}f(\widehat{𝐫},\kappa )=\mathrm{\hspace{0.17em}1},$$
(10)
with $`Z`$ the bare charge.

For a circle of radius $`r_0`$ and uniform line charge we find

$$f(\theta ,\kappa r_0)=\text{I}_0(\kappa r_0\mathrm{sin}\theta ),$$
(11)
where $`\theta `$ is the angle between $`𝐫`$ and the normal to the plane containing the circle. Finally, for a uniformly charged rod of vanishing radius, and length $`2l`$:

$$f(\theta ,\kappa l)=\frac{\mathrm{sinh}(\kappa l\mathrm{cos}\theta )}{\kappa l\mathrm{cos}\theta },$$
(12)
where $`\theta `$ is now the angle between $`𝐫`$ and the rod direction.

## IV Towards charge renormalization for platelets

The bare charge of Laponite platelets can be considered to be of the order of $`Z_{\text{bare}}=700`$ to 1000 negative elementary charges . At the level of PB approximation, this corresponds to a high electrostatic coupling ($`\varphi `$ larger than unity ), where the linearization procedure underlying Eq. (1) fails. However, a few Debye lengths away from the colloid, the potential has sufficiently decreased so that (1) is recovered and the one-body potential takes the form

$$\varphi (r,\theta )=Z_{\text{eff}}f(\theta ,\kappa r_0)\frac{e^{\kappa r}}{r}.$$

The effective charge $`Z_{\text{eff}}`$ is (in absolute value) smaller than the bare one , and within PB theory, saturates to a value independent of $`Z_{\text{bare}}`$ if the latter quantity is large enough. The contribution embodied in $`f(\theta ,\kappa r_0)`$ results from the anisotropy of the colloid, and is a priori an unknown function. In the limit of small electrostatic coupling (small $`Z_{\text{bare}}`$) where $`Z_{\text{eff}}/Z_{\text{bare}}1`$ by definition, the results of section II show that $`f(\theta ,X)=2\text{I}_1(X\mathrm{sin}\theta )/(X\mathrm{sin}\theta )`$. On the other hand, at higher $`Z_{\text{bare}}`$ corresponding to the saturation regime, the functional form of $`f(\theta ,X)`$ is the signature of the effective charge distribution on the platelet and related to possible differences in counter-ion condensation around the centre of the discs or in the vicinity of the edges.

It has been shown within PB theory that in the colloidal limit $`\kappa r_0>1`$, highly charged colloids behave as far as their far field is concerned as constant potential particles whose value is close to $`4kT/e`$ (i.e. $`\varphi =4`$), irrespective of shape (planar, cylindrical, spherical) . This prescription may be applied to the present case, with the restriction that to our knowledge, the linearized Poisson-Boltzmann problem cannot be solved analytically with the boundary condition of constant surface potential. However, the potential associated with the constant surface charge boundary condition provides a reasonable first approximation providing the correct qualitative variation of the effective charge with physico-chemical parameters. The corresponding anisotropic contribution to the large distance field is that computed in section II and the calculation of $`Z_{\text{eff}}`$ follows from the knowledge of $`\varphi (r=0)`$. From the analytical expression reported in , we obtain the saturation value of $`Z_{\text{eff}}`$

$$Z_{\text{eff}}^{\text{saturation}}=Z_{\text{sat}}=\frac{r_0}{\mathrm{}_B}\frac{2\kappa r_0}{1\mathrm{exp}(\kappa r_0)}.$$
(13)
It is noteworthy that for a typical salt concentration corresponding to $`\kappa r_0=1`$, $`Z_{\text{sat}}100`$, which is an order of magnitude smaller than the bare charge. This asymmetry guarantees that the effective charge is in the saturation regimes where it no longer depends on the value of $`Z_{\text{bare}}`$  (the precise experimental determination of $`Z_{\text{bare}}`$ is therefore not required). This saturation picture is expected to be reliable in a 1:1 electrolyte (PB mean-field approximation generally fails in electrolytes of higher valence), when there is a clear scale separation between the Bjerrum length and the size of the charged object (this constraint being fulfilled for Laponite) .

The above prediction may be checked by inserting $`Z_{\text{sat}}`$ into the analytical expression giving the LPB force between two parallel platelets  (configuration depicted in Fig. 2):

$$\frac{r_0F_z}{kT}=\frac{4\pi \mathrm{}_B}{r_0}Z_{\text{sat}}^2_0^{\mathrm{}}J_1^2(x)\frac{1}{x}\mathrm{exp}\left\{(h/r_0)\sqrt{x^2+\kappa ^2r_0^2}\right\}𝑑x$$
(14)
and comparing the results with the Monte Carlo simulations of Meyer et al , who considered the situation of vanishing salt. In this limit, where the quality of our prescription for $`Z_{\text{eff}}`$ is expected to deteriorate, we get from (13) $`Z_{\text{sat}}2r_0/\mathrm{}_B42`$. The comparison is displayed in Fig. 3 which shows a good agreement. Neglect of charge renormalization effects lead to an overestimation of the force by more than two orders of magnitude \[more precisely, by a factor $`(Z_{\text{bare}}/Z_{\text{sat}})^2`$, see the difference between dashed and dotted curves in Fig 3\], which points to the prime importance of such a phenomenon.

## V Tentative investigation into the phase behaviour

At this point, it is interesting to investigate at least qualitatively the difference in the phase behaviour between spherical and disc-like charge stabilized colloids. For spheres, the fluid-solid transition driven by repulsive electrostatic interactions favors the isotropic fluid upon addition of salt in the solution . The opposite is observed for Laponite dispersions, where an increase of the ionic strength lowers the density where the “solid” phase appears . Given that van der Waals interactions are believed to be irrelevant in the corresponding parameter range , this effect (illustrated in Fig. 4), appears at first to contradict standard DLVO phenomenology, where a screening of electrostatic repulsion (decrease of Debye length $`\kappa ^1`$) is expected to promote the simple fluid phase (sol), as happens for spheres. However, the effective charge of an arbitrary charged object is generically a growing function of salt concentration \[see for instance expression (13), which increases with $`\kappa `$\], as a result of an enhanced screening of colloid/micro-ions attraction which diminishes the amount of counter-ion “condensation”. This increase of the effective charge favours the solid phase and is thus antagonistic to the abovementioned decrease of the range of electrostatic repulsion driven by the decrease of $`\kappa ^1`$. The competition between the increase of the amplitude of repulsion and the decrease of its range upon adding salt is a possible scenario to interpret the phase behaviour of discs, within the DLVO picture and without any attractive interactions. More quantitative results are given below.

In the following analysis, we treat spheres and discs on equal footings (both are considered to have the same radius $`r_0`$). From the above analysis, we consider the following pair potential for discs:

$$V_{AB}(r,\theta _A,\theta _B)=Z_{\text{sat}}^2\mathrm{}_B\mathrm{\hspace{0.17em}4}\frac{\text{I}_1(\kappa r_0\mathrm{sin}\theta _A)}{\kappa r_0\mathrm{sin}\theta _A}\frac{\text{I}_1(\kappa r_0\mathrm{sin}\theta _B)}{\kappa r_0\mathrm{sin}\theta _B}\frac{e^{\kappa r}}{r},$$
(15)
where $`Z_{\text{sat}}`$ is given by Eq. (13); its counterpart for charged spheres is the standard isotropic DLVO expression

$$V_{AB}(r)=Z_{\text{sat}}^2\mathrm{}_B\left(\frac{e^{\kappa r_0}}{1+\kappa r_0}\right)^2\frac{e^{\kappa r}}{r},$$
(16)
where the saturation value of the effective charge has been derived analytically in  as a function of salt concentration \[we again consider colloids with a high bare charge such that $`Z_{\text{eff}}`$ coincides with its saturation value\].

Given that the phase behaviour of particles interacting through Yukawa-like potentials has been extensively explored by computer simulations , we can readily obtain the melting density corresponding to a given salt concentration for the potential (16). The results (corresponding to spherical colloids) are shown in Fig. 5. A related procedure may be used to estimate qualitatively the position of the melting or freezing line for discs, by assuming that at low clay density (typically where the fluid/solid transition takes place for Laponite suspensions), the rotational motion of the platelets occurs on a shorter time scale than the translational one, so that these objects feel an effective potential resulting from the angular average of the expression given in Eq. (15), namely

$`V_{AB}^{\text{average}}(r)`$ $`=`$ $`V_{AB}(r,\theta _A,\theta _B)_{\theta _A,\theta _B}={\displaystyle \frac{1}{4}}{\displaystyle _0^\pi }\mathrm{sin}\theta _Ad\theta _A{\displaystyle _0^\pi }\mathrm{sin}\theta _Bd\theta _BV_{AB}(r,\theta _A,\theta _B)`$ (17)
$`=`$ $`Z_{\text{sat}}^2\mathrm{}_B\mathrm{\hspace{0.17em}4}\left({\displaystyle \frac{\mathrm{cosh}(\kappa r_0)1}{\kappa ^2r_0^2}}\right)^2{\displaystyle \frac{e^{\kappa r}}{r}}.`$ (18)
In doing so, we obtain an isotropic Yukawa potential where the energy scale –the term in parenthesis in (18)– reflects the original anisotropy of the pair potential. In the standard DLVO potential (16) for spheres, the energy scale \[the term $`\mathrm{exp}(\kappa r_0)/(1+\kappa r_0)`$\] also depends on Debye length, but has a different physical origin and results from the exclusion of micro-ions from the interior of the spheres. Making use of the numerical Yukawa phase diagram , we obtain the melting line represented in Fig. 5 for the averaged potential (18), where the threshold density decreases when the salt concentration (or equivalently $`\kappa `$) increases, at least for $`\kappa d<4`$. At the same level of description, spherical colloids show the opposite behaviour, for all values of $`\kappa d`$ (we also emphasize that for the spherical colloids used in the experiments of Monovoukas and Gast , the equivalent of the melting line reported in Fig. 5 is in excellent agreement with its experimental counterpart ).

For both spheres and discs, the prefactor of the Yukawa term $`\mathrm{exp}(\kappa r)/r`$ is an increasing function of $`\kappa `$ \[see Eqs. (16) and (18)\]. The subtle interplay between this increase and the decrease of Debye length is nevertheless able to produce a “re-entrant” melting line for discs only. In the limit of low $`\kappa r_0`$, both expressions (16) and (18) become $`\mathrm{}_BZ_{\text{sat}}^2\mathrm{exp}(\kappa r)/r`$; the melting lines for spheres and discs however do not merge in this limit in Fig. 5, due to the difference in the limiting values of $`Z_{\text{sat}}`$ for both geometries.

This approach consisting in averaging the two-body platelet potential over angular degrees of freedom predicts a qualitative change in the phase behaviour of platelet systems for a reduced density $`\rho ^{}=\rho d^3`$ of order 1, which corresponds to a clay mass fraction of 8% for Laponite with diameter $`d=300`$Å , see the upper $`x`$-label of Fig. 5. This density is approximately 4 times higher than the maximum density delimiting the fluid and solid regions in the phase diagram of Laponite suspensions (Fig. 4 and references ), so that the present approach does not allow a quantitative comparison with experiments. It is however noteworthy that a density $`\rho ^{}1`$ is much smaller than the isotropic/nematic coexistence density for uncharged plates ($`\rho ^{}4`$, ). At this stage, it is impossible to be more specific concerning the nature of the “solid” phase supplementing the fluid one at high densities. This question will be addressed in the following section by computer simulations.

## VI Mesostructure: Monte Carlo simulations

We have implemented standard Monte Carlo simulations with typically $`N=500`$ platelets interacting through the potential (15) where $`Z_{\text{sat}}`$ is given by Eq. (13). This allows to test the validity of the approach proposed in section V, where the initial anisotropic potential was mapped onto the isotropic Yukawa function (18). A typical run consisted in $`10^5`$ cycles (both random displacement and rotation of the $`N`$ particles).

For different values of salt concentration and platelet density, two diagnostics were used to characterize the simulated samples. First, the centre-to-centre pair distribution function $`g(r)`$ was computed. Second, the orientational ordering was quantified through the statistical average of (twice) the second Legendre polynomial $`P(\psi )=3\mathrm{cos}^2\psi 1`$ at a given centre to centre distance $`r`$, where $`\psi `$ is the angle formed by the normals to the two discs. The orientational pair correlation function denoted $`g_{or}(r)`$ follows from averaging over all pairs of platelets.

The corresponding information on the mesostructure is displayed in Figs. 6 and 7 for different densities and two values of $`\kappa `$. A striking observation is that a slight increase of the density induces an important increase of the structure (see the differences between the curves at $`\rho ^{}=0.82`$ and $`\rho ^{}=0.85`$ in Fig. 6; the same is seen at $`\kappa d=3`$ in Fig. 7 when the reduced density changes from 1.07 to 1.10). Such values of $`\rho ^{}`$ lie close to the threshold $`\rho ^{}=1`$ estimated in section V. The associated orientational distributions (not shown) are structure-less for the densities considered here. These results indicate a qualitative change in the structure of the fluid phase upon increasing the density, but it is still difficult to characterize the new “phase” emerging without a thermodynamical study.

We have tested the relevance of the scenario propounded at the beginning of section V (competition between an increase of the amplitude and a decrease of the range of the pair potential when $`\kappa `$ –or equivalently $`n_{\text{salt}}`$– increases). For a given density ($`\rho ^{}=1.1`$), the pair correlation functions are monitored for various values of $`\kappa `$. The corresponding charges change with $`\kappa `$, according to Eq. (13). Figure 8 shows an interesting feature. As expected, the structure is most pronounced for the lowest $`\kappa `$, and decreases when $`\kappa `$ increases (see the difference between the curves for $`\kappa d=1`$ and $`\kappa d=2`$). However, when salt concentration is further raised to $`\kappa d=3`$, the maximum of $`g(r)`$ and the whole structure are enhanced, before decreasing again when $`\kappa d=4`$. The effect evidenced in Fig. 8 is reminiscent of the “anomalous” slope of the phase diagram reproduced in Fig. 4, and in qualitative agreement with the simplified phase diagram drawn in Fig. 5.

Finally, we have performed exploratory runs at higher densities where a nematic phase would be observed in the uncharged system. The results of Fig. 9 for a high density $`\rho ^{}=5.0`$ show a rich local structure both for the $`g(r)`$ and the orientational $`g_{or}(r)`$, very different from that observed at lower densities. The reference uncharged system ($`q^{}=0`$) with its strong nematic plateau is displayed for comparison. Inclusion of electrostatic interactions preserves the long range nematic order, but decreases its strength (the height of the plateau). Such a nematic order was clearly absent at lower densities (see Fig. 7). The peaks of $`g(r)`$ appear correlated to those of $`g_{or}(r)`$, which points to the existence of oriented micro-domains with parallel platelets, and a higher nematic order than the mean one. The complex behaviour of the distribution functions of Fig. 9 calls for more thorough investigations at high densities. In particular the validity of the large distance expansion (15) may be questionable at high $`\rho ^{}`$.

## VII Conclusion

In this article we devise an electrostatic pair potential for both highly anisotropic and highly charged objects (of vanishing internal volume) based on a linearized Poisson-Boltzmann theory. On the example of infinitely thin platelets, we show that the anisotropic shape of the object results in an anisotropic Yukawa potential even at large distances. The anisotropy is an increasing function of $`\kappa r_0`$ and becomes very marked in the colloidal limit $`\kappa r_0>1`$. To account for the non-linearities in the Poisson-Boltzmann theory, when the bare charge of the platelets is large, we introduce an effective or renormalized charge. This approximation yields an effective pair potential for highly charged platelets.

The addition of salt generically results in two antagonists effects, irrespective of colloid geometry: the range of the effective potential decreases, but its amplitude increases. This interplay appears quite subtle, and able to discriminate between spherical and disc-like colloids by producing a re-entrant melting line for discs. This phenomenon appears at first incompatible with naive DLVO expectation. Moreover, the corresponding threshold density for discs is smaller than the isotropic/nematic coexistence for uncharged plates, but still higher than its experimental “Fluid/Solid” counterpart.

Preliminary results of Monte-Carlo simulations for particles interacting through an anisotropic potential provide informations on the mesoscopic structure of the dense phase. In particular, we show that the resulting structure is very sensitive upon increasing the density. We have also observed a “re-entrant” melting curve. These effects are in qualitative agreement with the experimental phase diagram of Laponite suspensions found in the literature, even if this topic is still under debate.

Our approach suffers from several weaknesses. a) the non-linearities in the Poisson-Boltzmann theory are accounted for at the level of charge renormalization, b) in this framework, we expect the effective charge to be better approximated with an anzatz where the particle effectively behaves as a constant electrostatic potential object , whereas we considered the platelets to be of constant surface charge, c) we work at the level of pair potential and correlation effects beyond the mean-field Poisson-Boltzmann approach might play a role, d) the precise role of van der Waals forces has not been assessed, and remains obscure for clay platelets. Further theoretical and numerical work is needed, but we hope that the arguments presented here provide a first hint into the full problem. At least, our results strongly suggest that the combined effect of anisotropy and charge condensation have a significant and non-trivial qualitative influence on the phase diagram of highly charged colloids. These features could well be of prime importance in our understanding of the thermodynamics of clay suspensions.

Acknowledgments. We thank J.-P. Hansen, H. Lekkerkerker, P. Levitz, A. Delville, D. Bonn, T. Nicolai, B. Jancovici and F. van Wijland for fruitful discussions.

## A

In this appendix, we show that the solution of Eq. (1) is not given by convolution product (2) for polyions having a non zero internal volume. For simplicity, we consider a charged hard sphere of radius $`r_0`$ (which consequently excludes micro-ions from its interior), with a uniform surface charge $`\sigma `$. The solution of Eq. (1) takes the well known DLVO form

$$\varphi (r)=Z\mathrm{}_B\frac{e^{\kappa r_0}}{1+\kappa r_0}\frac{e^{\kappa r}}{r},$$
(A1)
where $`Ze=4\pi r_0^2\sigma `$ denotes the total bare charge. Alternatively, the convolution route of Eq. (2) yields

$`\varphi _{\text{convol}}(r)`$ $`=`$ $`\mathrm{}_B{\displaystyle _𝒫}{\displaystyle \frac{\sigma }{e}}{\displaystyle \frac{\mathrm{exp}(\kappa |𝐫𝐬|)}{|𝐫𝐬|}}d^2𝐬`$ (A2)
$`=`$ $`Z\mathrm{}_B{\displaystyle \frac{\mathrm{sinh}(\kappa r_0)}{\kappa r_0}}{\displaystyle \frac{e^{\kappa r}}{r}}\text{for}rr_0.`$ (A3)
The origin of this discrepancy comes from the fact that upon summing up the elementary surface contributions in (A2), the micro-ions are allowed to enter the “interior” region $`rr_0`$. The field created is thus that of the shell with charge $`Z`$ (the colloid) plus that of the “plasma” inside. The corresponding interior charge $`Z_{\text{in}}`$ is easily computed from the expression of the electrostatic potential for $`rr_0`$, which follows from a simple permutation of $`r`$ and $`r_0`$ in expression (A3)

$$\varphi _{\text{convol}}^{\text{inside}}(r)=Z\mathrm{}_B\frac{\mathrm{sinh}(\kappa r)}{\kappa r}\frac{e^{\kappa r_0}}{r_0}.$$
(A4)
We obtain

$`Z_{\text{in}}`$ $`=`$ $`{\displaystyle \frac{\kappa ^2}{4\pi \mathrm{}_B}}{\displaystyle _0^{r_0}}\varphi _{\text{convol}}^{\text{inside}}(r)d^3𝐫`$ (A5)
$`=`$ $`Ze^{\kappa r_0}\left[\mathrm{cosh}(\kappa r_0){\displaystyle \frac{\mathrm{sinh}(\kappa r_0)}{\kappa r_0}}\right].`$ (A6)
Gathering results, a straightforward calculation allows to rewrite Eq. (A3) in the form

$$\varphi _{\text{convol}}(r)=Z\mathrm{}_B\frac{\mathrm{sinh}(\kappa r_0)}{\kappa r_0}\frac{e^{\kappa r}}{r}=\left(Z+Z_{\text{in}}\right)\mathrm{}_B\frac{e^{\kappa r_0}}{1+\kappa r_0}\frac{e^{\kappa r}}{r},$$
(A7)
so that the DLVO structure of expression (A1) is recovered.

The argument given here bears some similarities with the proof of the equivalence between two charged hard spheres models with a uniform background  (note that the equivalent of the background, the micro-ion density, is not uniform in the present situation). Finally, a related remark is that application of expression (4) for spheres does not give the DLVO pair potential (16) even if the correct one body potential (A1) is used \[use of expression (A2) is equally incorrect\]. A simple way to recover (16) is through the stress tensor .


format: text

text:
# Stacking-fault energies for Ag, Cu, and Ni from empirical tight-binding potentials

## I Introduction

The stacking-fault energy is an important property of a material since it determines to a high degree its deformation and failure behavior.Hirth and Lothe (1968) Therefore, for the computer modeling of processes involving plastic deformation and/or failure, a correct representation of the stacking-fault energy by the employed model is desirable. However, many models have been constructed without consideration of stacking-fault energies and the predicted values deviate often considerably from experimental results.Heino et al. (1999); Zimmerman et al. (2000); Mehl et al. (2000) A precise knowledge of the stacking-fault energies predicted by empirical models is therefore of high interest in order to correctly interpret results of computer simulations.

The empirical tight-binding potentials of Cleri and Rosato Cleri and Rosato (1993) are extensively used in atomistic simulations since they give a reasonable description of many structural and thermal properties of fcc metals. In this work we study the intrinsic stacking-fault energies and free energies predicted by these potentials for Ag, Cu, and Ni and compare them with the results of similar models for Cu and Ni. The stacking-fault energies were calculated at $`T=300`$ K from molecular-dynamics simulations and at $`T=0`$ K with the help of relaxation calculations (molecular-statics). The derivation of the vibrational density of states (VDOS) from the simulations made it also possible for us to calculate the stacking-fault free energies within the quasiharmonic approximation and to study the local effect of stacking-faults in fcc metals.

The aim of this paper is to provide the stacking-fault energies for Ag, Cu, and Ni modeled by the potentials of Cleri and Rosato. We expect that these values will contribute to proper discussions of results based on these potentials. Although in the case of Ag and Cu the values derived from the simulations are significantly lower than the experimental results, the trend between all three elements is well reproduced. It has to be noted that the occurrence of low stacking-fault energies is a general problem of central-force potentials. One reason for this is that without the consideration of angular forces significant contributions to the stacking-fault energy arise from the third and higher neighbor shells, only.Girshick et al. (1998) A related explanation is given by the nearest neighbor d-band tight-binding model where the energy difference between the fcc and hcp structure is dominated by the fifth and sixth moment of the density of states,Sutton (1993) whereas many empirical potentials, including those of Cleri and Rosato, are based on the second moment approximation. Better agreements of computationally simple empirical potentials with experimental values for the stacking-fault energy require therefore the explicit consideration of the stacking-fault energy during the construction of the potentials.

## II Computational Methods

### II.1 General

As stated in the introduction, the stacking-fault energies and free energies reported in this article have been derived from molecular-dynamics simulations employing the potentials of Cleri and Rosato.Cleri and Rosato (1993) In these simulations the equations of motions have been integrated with the help of the velocity form of the Verlet-algorithm Allen and Tildesley (1991) using a time-step of 2 fs. Isothermic-isochoric (NVT ensemble) and isobaric-isoenthalpic (NpH ensemble) conditions were achieved where necessary by the application of the Nosé-Hoover thermostat Hoover (1985) method and the Parrinello-Rahman scheme Parrinello and Rahman (1980) (restricted to orthogonal simulation boxes), respectively. Periodic boundary conditions were applied in all cases.

Calculations of the relaxed structures at $`T=\text{0\hspace{0.17em}K}`$ were done in two steps. First, the systems were cooled down to a temperature of $`T=\text{0.1\hspace{0.17em}K}`$. Afterwards, the final configurations were obtained with the help of a steepest-decent method.

For each of the three elements a stacking-fault configuration and a reference configuration were prepared. The stacking-fault configurations were composed as a stack of 110 triangular lattice layers building a nearly cubic system of $`N_{\mathrm{st}}=\mathrm{1\hspace{0.17em}018\hspace{0.17em}160}`$ atoms with one stacking-fault; the reference configurations on the other hand are regular fcc systems$`N_{\mathrm{fcc}}=\mathrm{1\hspace{0.17em}000\hspace{0.17em}188}`$ atoms ($`63^3`$ cubic cells).

At first glance, configurations with more than a million atoms may seem to be exaggerated for the simple calculation of stacking-fault energies. However, these system sizes guarantee the presence of a sufficiently large number of vibrational modes for the calculations of free energies in the quasiharmonic approximation. In addition, the large system sizes lead to large distances between the periodic copies of the stacking-fault, thereby minimizing possible self-interaction effects.

Before the calculation of the stacking-fault energies could begin, the equilibrium dimensions of all configurations had to be determined at $`T=0`$ and 300 K. In the case of the reference configurations these were derived from the $`T=\text{0\hspace{0.17em}K}`$ lattice constants given in Ref. Cleri and Rosato, 1993 and the average system volumes obtained from isobaric-isoenthalpic simulations at $`T=\text{300\hspace{0.17em}K}`$ and zero pressure.

For the determination of the equilibrium sizes of the stacking-fault configurations, it had to be taken into account that only volume relaxations perpendicular to the plane of the stacking fault are physically relevant. Relaxations in the directions parallel to the fault-plane are artifacts of the finite system size and diminish if the number of atomic layers parallel to the fault is increased. For this reason the in-plane dimensions of the simulation boxes of the stacking-fault configurations were derived from the lattice constants obtained for the perfect systems and only the out-of-plane dimensions had to be explicitly determined. At $`T=\text{300\hspace{0.17em}K}`$ these missing dimensions were again obtained with the help of zero-pressure isobaric-isoenthalpic simulations in which only the out-of plane direction of the simulation box was allowed to fluctuate. Similarly the volume relaxations at $`T=\text{0\hspace{0.17em}K}`$ were derived from relaxations of the stacking-fault configurations with fixed in-plane dimensions of the simulation box.

### II.2 Calculation of stacking-fault energies

After the dimensions of the configurations had been determined, the stacking-fault energies $`\gamma `$ of the three metals could be obtained. For this the simple formula

$$\gamma =\frac{N_{\text{sf}}(e_{\text{sf}}e_{\text{fcc}})}{A},$$
(1)
was used, where $`A`$ is the area of the stacking-fault plane, $`N_{\text{sf}}`$ the number of atoms in the stacking-fault configuration and $`e_{\text{sf}}`$, $`e_{\text{fcc}}`$ the potential energies per atom in the respective configuration. The kinetic energy does not contribute to the stacking-fault energy since its per particle value does not depend on the configuration.

Particular attention has to be paid to the system temperature. The usual method to attain the desired temperature in molecular-dynamics simulations is to rescale the particle velocities repeatedly until the system is equilibrated. In simulations run under conditions of constant energy or enthalpy, this procedure may lead to a small offset in the temperature which in turn affects the value of the potential energy. While these discrepancies usually are very small and can be neglected in most cases, they can have an dramatic impact on the values of $`\gamma `$ since the potential energy difference on the right hand side of Eq. (1) is much smaller than the potential energies themselves and can be comparable to — or even smaller than — the error in the potential energies induced by temperature deviations. In order to avoid this problem we derived the potential energies from simulations in the canonical (NVT) ensemble.

### II.3 Calculation of stacking-fault free energies

The exact calculation of free energies from molecular-dynamics simulations is a difficult task. In this work we employ the quasiharmonic approximationLewis and Klein (1990) to derive the free energies $`F`$ of the stacking-fault and reference configurations. Within this approach the free energy is given by

$$F=E_0(V)+3N_0^{\mathrm{}}g(\omega )\left[\frac{\mathrm{}\omega }{2}+\frac{1}{\beta }\mathrm{log}(1\mathrm{e}^{\beta \mathrm{}\omega })\right]𝑑\omega ,$$
(2)
where $`\beta `$ is the inverse temperature $`(k_\mathrm{B}T)^1`$, $`g(\omega )`$ the normalized VDOS and $`E_0(V)`$ denotes the ground-state energy of the system *for the equilibrium volume at temperature $`T`$*. In contrast to most molecular-dynamics simulations, Eq. (2) takes quantum effects fully into account. However, since we are mainly interested in the values of the stacking-fault energies in a classical computer-simulation context, it is more appropriate to calculate the free energy in the classical limit. Taking this limit and neglecting the purely quantum zero-point fluctuations, Eq. (2) becomes

$$F_{\mathrm{cl}}=E_0(V)+3N\beta ^1_0^{\mathrm{}}g(\omega )\mathrm{log}(\beta \mathrm{}\omega )𝑑\omega .$$
(3)
Given Eq. (3) the task of calculating the free energy of the configurations by molecular-dynamics simulations reduces to the determination of the vibrational density of states $`g(\omega )`$ and the ground-state energy $`E_0(V)`$. The latter quantity was obtained by a relaxation of the systems to $`T=\text{0\hspace{0.17em}K}`$ using the simulation box dimensions of $`T=\text{300\hspace{0.17em}K}`$. In order to derive the vibrational density of states $`g(\omega )`$, simulations in the NVE ensemble were performed at $`T=\text{300\hspace{0.17em}K}`$. From the trajectories of the particles during these simulations the velocity-autocorrelation function $`𝐯(t)𝐯(0)`$ was derived which is related to the normalized VDOS byLovesey (1986)

$$g(\omega )=_{\mathrm{}}^{\mathrm{}}\frac{𝐯(t)𝐯(0)}{𝐯(0)𝐯(0)}\mathrm{e}^{i\omega t}𝑑t.$$
(4)
In order to reduce statistical errors, the density of states $`g(\omega )`$ was averaged over ten or more simulation runs of 5 000 steps, each. From the density of states obtained this way the free energies of all six configurations could be calculated with the help of Eq. (3). These values were then used to calculate the intrinsic stacking-fault free energies $`\stackrel{~}{\gamma }`$ of the three elements by replacing the potential energies per particle in Eq. (1) by the corresponding free energies per particle. Since at $`T=\text{0\hspace{0.17em}K}`$ energy and free energy are identical these calculations had to be done only at $`T=\text{300\hspace{0.17em}K}`$.

## III Results

### III.1 Stacking-fault Energies

In Table 1 we present the stacking-fault energies and free energies obtained from molecular-dynamics simulations as described in the preceding section. From this table it can be seen that the general agreement with the experimental results is only qualitative, similar to comparable models. The calculated stacking-fault energies (and free energies) for Cu and Ag are much smaller than their experimental counterparts; the calculated stacking-fault free energy of Ag at room temperature is even negative. For Ni the situation is not clear due to the large variation of the experimental data. However, the trend among the three elements is well reproduced by the model. This means that while one has to be careful with the interpretation of quantitative results, the potentials of Cleri and Rosato are well suited for comparisons between materials with high, low and very low stacking-fault energies.

In accordance with the experimental data,Murr (1975) the results in Table 1 show a decrease of the stacking-fault energies and free energies with increasing temperature. In the case of the stacking-fault energies, the decrease is caused by anharmonic effects only, since the specific heat of a purely harmonic, classical system does not depend on structural details. In contrast to this, the decrease of the stacking-fault free energies in Table 1 contains harmonic and anharmonic contributions. However, due to the nature of the quasiharmonic approximation, the anharmonic contributions are treated only in an approximate fashion via the volume dependence of $`E_0`$ and the temperature dependence of the VDOS. This approximation gets in the way if one wants to compare the stacking-fault energies and free energies in order to determine the role of entropic contributions to the stacking-fault energy. Particularly striking is the large discrepancy between the stacking-fault energy and free energy of Ni at room temperature. Since our results indicate that the entropy difference favors the stacking-fault configuration, the stacking-fault free energy should be the lower quantity, as it is the case for Cu and Ag. This shows that, at least in the case of Ni, there are important anharmonic contributions to the stacking-fault free energy which are not accounted for by the quasiharmonic approximation.

### III.2 Vibrational Density of States

Figure 1 shows the VDOS calculated for the fcc structure as well as the difference of the VDOS between the stacking-fault and the perfect fcc configuration. It should be noted that the absolute value of this difference is physically irrelevant since it depends on the system size and goes to zero in the limit of an infinitely thick system. However, the shape of this function is system size independent and reveals the energy regions where the VDOS is affected by the stacking fault. From the lower panel of Fig. 1 it can be seen that the strongest effect of the stacking fault appears in the transverse peak where modes are transferred from the center of the peak to (mostly) lower energies. A similar but smaller transfer takes place in the longitudinal peak. It is these transfers which give rise to the entropy difference between the configurations, which in turn accounts for the temperature dependent change of the stacking-fault free energy within the quasiharmonic approximation. A similar behavior is found for Ag and Ni with a slightly more pronounced transfer of longitudinal modes in the case of Ni.

In order to see the spatial extent of the influence of stacking faults on the vibrational properties of the crystal lattice we have calculated the local VDOS for the first atomic layers near a stacking fault in Cu. The results of this calculation, which are presented in Fig. 2, show that atoms in the first and second layer near the fault account for most of the change of the total VDOS in Fig. 1, whereas the VDOS in the fourth layer differs only slightly from the bulk VDOS. Atoms in the third layer (whose VDOS is not given in order to avoid overcrowding of the figure) have an intermediate character. While their VDOS is still clearly different from the VDOS of the crystalline bulk, it is already much closer to this than to the VDOS of the first two layers.

## IV Summary and Conclusions

The stacking-fault energies predicted by the empirical tight-binding potentials of Cleri and Rosato reflect the experimental values of Ag, Cu, and Ni qualitatively. Similar to other empirical models they underestimate the values of Ag and Cu. In contrast to this and in rough agreement with experimental data, rather high stacking-fault energies are found for Ni. Since the experimental data for Ni show strong variations it is impossible to say whether the potential for this metal overestimates the stacking-fault energy or not. In any case, the potentials can be used for qualitative comparisons of materials with high, low, and very low stacking-fault energies.

For all three metals we find a decrease of the stacking-fault energy with increasing temperature of the same order of magnitude as it is experimentally observed.Murr (1975) This decrease has to be considered in the construction or evaluation of models which explicitly take the stacking-fault energies into account.

The discrepancy between the calculated room temperature stacking-fault energy and free energy for Ni shows that important anharmonic effects are missed by the quasiharmonic approximation. For this reason one should not put too much emphasis on the negative value found for the room temperature stacking-fault free energy of Ag. What remains is the question why the anharmonic effects are so much stronger for Ni. A possible explanation might be given by the relative importance of the many-body term vs. the pair-potential term in the model. An inspection of the parameters of the potentials in Ref. Cleri and Rosato, 1993 shows indeed that the pair-potential term has a lower weight and a shorter range for Ni than for Cu or Ag. Thus, the increased weight of the many-body term in the case of Ni might account for the increased occurrence of anharmonic effects in this metal.

Finally, the analysis of the local effect of stacking-faults on the vibrational properties in fcc metals reveals a strong transfer of modes to lower energies in the transverse peak and a similar but less pronounced transfer in the longitudinal peak of the VDOS. Inspection of the layer-resolved local VDOS shows that these effects are localized in a region of six to eight monolayers around the defect.

###### Acknowledgements.

This work has been supported by grants from the Canadian Natural Sciences and Engineering Research Council (NSERC), Québec’s Fonds pour la formation de chercheurs et l’aide à la recherche (FCAR) and Germany’s Deutsche Forschungsgmeinschaft (DFG). Most of the calculations have been performed on the facilities of the Réseau québécois de calcul de haute performance (RQCHP).


text:
# Hydrodynamics with spontaneous symmetry breaking: application to relativistic heavy ion collisions

## I Introduction

Relativistic hydrodynamics is a very useful tool in the field of heavy ion collisions when one wants to reproduce the space-time evolution of the hot and dense matter created in such events. Upon the hypothesis of local equilibrium, which can be questionable for such small and not long-lived systems, this approach directly gives in terms of usual quantities such as temperature, pressure…, an intuitive description of what happens during the whole evolution of the system until the freeze-out takes place. Moreover, it is usually assumed that the matter under consideration is almost perfect, so that standard hydrodynamics can be applied directly. However, the hot and dense matter produced is mainly constituted of pions and, with the exception of son , no attention has been drawn to the fact that the breaking of chiral symmetry can affect the theory of hydrodynamics itself. One of the aims of this paper is thus to give in a simple approach a quantitative estimate of this modification. More precisely, we will determine the temperature profile with respect to time in the framework of the (1+1) dimensional boost-invariant Bjorken model with and without the effect of the spontaneous breaking of the symmetry. In a second step we also introduce the dissipation process and compare the previous results with those obtained when Müller-Israel-Stewart second order theory of dissipation is included on top of standard hydrodynamics.

## II Hydrodynamics equations for chiral fluid

Following son , the hydrodynamics degrees of freedom in a chiral fluid with chiral $`SU(2)`$ symmetry spontaneously broken are the densities associated to conserved quantities, namely entropy density $`s`$, momentum density $`T^{0i}`$, baryonic number density $`n`$, left and right-handed charge densities which can be written as $`SU(2)`$-matrices $`\rho _L\rho _L^i\tau _i/2`$ and $`\rho _R\rho _R^i\tau _i/2`$ and finally the variables associated to the Goldstone modes introduced in the theory via a $`SU(2)`$-matrix $`\mathrm{\Sigma }e^{i\stackrel{}{\tau }.\stackrel{}{\pi }/f_\pi }`$. If we focus on the central rapidity region of a heavy ion collision, we can neglect in a first step the baryonic density. Thus, the hydrodynamic equations taken from son  become:

$$_\mu T^{\mu \nu }=0$$
(1)
$$i_\mu ((f_t^2f_s^2)u^\mu u^\nu \mathrm{\Sigma }_\nu \mathrm{\Sigma }^{}+f_s^2\mathrm{\Sigma }^\mu \mathrm{\Sigma }^{})=0$$
(2)
where the energy-momentum tensor $`T^{\mu \nu }`$ is given by $`T^{\mu \nu }=(ϵ+p)u^\mu u^\nu pg^{\mu \nu }+f_s^2/4\text{tr}(^\mu \mathrm{\Sigma }^\nu \mathrm{\Sigma }^{}+^\nu \mathrm{\Sigma }^\mu \mathrm{\Sigma }^{})`$.

In the above equations, $`f_t`$ ($`f_s`$), which can be interpreted as the temporal (spatial) pion decay constant, are actually not constants but have a non-trivial temperature dependence. Nevertheless, in order to have an estimation of the effects induced by the presence of the spontaneous breaking of the chiral symmetry, we will adopt two simple approaches. The first one, for the sake of simplicity, is to consider that $`f_t=f_s=f_\pi =93`$ MeV. Of course, at finite temperature, this is not true and, for instance, for small temperatures, one knows that $`f_t^2`$ and $`f_s^2`$ differ from $`f_\pi ^2`$ by a factor proportional to $`T^2`$ while the correction for $`f_t^2f_s^2`$ is of order $`T^4`$ so that, to first order, $`f_t^2=f_s^2=f_\pi ^2(1T^2/(6f_\pi ^2))`$. This second approach, more realistic, is based on this formula, extended to larger temperatures.

## III The Bjorken picture

Even if several important improvements concerning the geometric description of a heavy ion collision have been done since the pioneering paper of Bjorken Bjorken , we will use here this simplified version. The Bjorken picture of a collision is based on the hypothesis that initial conditions are invariant under Lorentz transformations throughout the central rapidity plateau. This implies that, for a one-dimensional flow in the $`z`$-direction, all the quantities depend actually not on the time $`t`$ and the position $`z`$ but only on the proper time $`\tau =\sqrt{t^2z^2}`$. The space-time evolution of the system can be then represented by hyperbola $`\tau =cst`$. If the hydrodynamic velocity $`u^\mu `$ is taken as $`u^\mu =x^\mu /\tau =(t/\tau ,0,0,z/\tau )`$, the hydrodynamic equations

$$_\mu T^{\mu \nu }=0$$
(3)
with $`T^{\mu \nu }=(ϵ+p)u^\mu u^\nu pg^{\mu \nu }`$ lead for an ultra relativistic gas to the following differential equations csernai ; CYW :

$$\frac{s}{\tau }+\frac{s}{\tau }=0$$
(4)
$$\frac{T}{\tau }+\frac{1}{3}\frac{T}{\tau }=0$$
(5)
$$\frac{ϵ}{\tau }+\frac{4}{3}\frac{ϵ}{\tau }=0$$
(6)
whose solutions give the following scaling behavior:

$$s(\tau )=s(\tau _0)\frac{\tau _0}{\tau }$$
(7)
$$T(\tau )=T(\tau _0)\left(\frac{\tau _0}{\tau }\right)^{1/3}$$
(8)
$$ϵ(\tau )=ϵ(\tau _0)\left(\frac{\tau _0}{\tau }\right)^{4/3}.$$
(9)
And, since $`p=ϵ/3`$:

$$p(\tau )=p(\tau _0)\left(\frac{\tau _0}{\tau }\right)^{4/3}.$$
(10)
## IV Results

We first take into account the effect of the spontaneous breaking of chiral symmetry where no renormalization of the pion decay constant occurs, i.e. $`f_t=f_s=f_\pi =93`$ MeV. We thus obtain from the hydrodynamic equations (1),(2) in the first non-vanishing order of the pion field :

$$\frac{s}{\tau }+\frac{s}{\tau }=0$$
(11)
$$\frac{T}{\tau }+\frac{1}{3}\frac{T}{\tau }=0$$
(12)
$$\frac{ϵ}{\tau }+\frac{Ts(\tau _0)\tau _0}{\tau ^2}\frac{1}{\tau }\underset{i=1}{\overset{3}{}}\left(\frac{\pi ^i}{\tau }\right)^2=0$$
(13)
$$\left(\frac{}{\tau }+\frac{1}{\tau }\right)\frac{\pi ^i}{\tau }=0.$$
(14)
and from the explicit expression of the pressure $`p`$ given in son :

$$p=aT^4+\frac{1}{2}\underset{i=1}{\overset{3}{}}\left(\frac{\pi ^i}{\tau }\right)^2$$
(15)
We see that there is no modification of the equation for the entropy density and for the temperature, that is the scaling behavior is exactly the same $`s(\tau )=s(\tau _0)\tau _0/\tau `$ and $`T(\tau )=T(\tau _0)\left(\tau _0/\tau \right)^{1/3}`$. On the contrary, for the energy density and for the pressure, one can observe a modification due to a term proportional to $`\pi ^i/\tau `$ which can be determined from the last equation (14). The solution for the time derivative is $`\pi ^i/\tau =A/\tau `$ where $`A`$ is a constant whose value is fixed by initial condition for $`\tau =\tau _0`$. If we now report in the above equations, we can solve analytically the resulting system to obtain:

$$ϵ(\tau )=ϵ(\tau _0)\left(\frac{\tau _0}{\tau }\right)^{4/3}\frac{3}{2}\frac{A^2}{\tau _0^2}\left(\frac{\tau _0}{\tau }\right)^2$$
(16)
$$p=p(\tau _0)\left(\frac{\tau _0}{\tau }\right)^{4/3}+\frac{3}{2}\frac{A^2}{\tau _0^2}\left(\frac{\tau _0}{\tau }\right)^2.$$
(17)
In order to fix the initial condition related to the value of $`A`$, let us now consider the density matrix $`\rho =\mathrm{exp}^{\beta H}`$. The part of the Hamiltonian relevant for chiral hydrodynamics degrees of freedom can be taken from son . An explicit calculation with Bjorken modelisation of the collision gives:

$$\rho =\mathrm{exp}\left(\beta \text{d}^3xT^{00}\right)\mathrm{exp}\left(\frac{\beta }{2}\text{d}^3x\underset{i=1}{\overset{3}{}}\left(\frac{\pi ^i}{\tau }\right)^2\right).$$
(18)
Assuming isospin symmetry, we then have: $`\rho \mathrm{exp}\left(\frac{3\beta }{2}\text{d}^3x\left(\frac{A}{\tau }\right)^2\right)`$. For $`\tau =\tau _0`$, we thus immediately see that $`\rho `$ is Gauss distributed. Actually, if we call $`V_0`$ the initial volume in the rest-frame of the two colliding nuclei and $`\gamma _{cm}`$ the Lorentz factor needed to go from that frame to the laboratory frame, we have:

$$\rho (\tau _0)\mathrm{exp}\left(\frac{3\beta }{2}\frac{V_0}{\gamma _{cm}}\left(\frac{A}{\tau _0}\right)^2\right).$$
(19)
We can easily estimate $`V_0`$ by geometric arguments: $`V_0=\pi R^22\tau _0`$ where $`R`$ is the radius of the colliding nuclei. If we now put typical numbers on all these constants: $`R5`$ fm, $`\tau _0=1`$ fm, $`\gamma _{cm}10`$ and $`T(\tau _0)=200`$ MeV, we can directly estimate the effect of the terms associated to the spontaneous symmetry breaking since $`<A^2>`$ is simply related to the Gaussian distribution: $`<A^2>=T(\tau _0)\tau _0\gamma _{cm}/(6\pi R^2)`$. Numerically, using these inputs, we obtain a correction of order $``$1% with respect to the perfect case.

We now introduce the explicit chiral perturbation $`T`$-dependence for $`f_t`$ and $`f_s`$, i.e. $`f_t^2=f_s^2=f_\pi ^2(1T^2/(6f_\pi ^2))`$. The evolution of the temperature is then governed by the following equation :

$$4a\tau T^3\frac{<A^2>}{\tau \left(1T^2/(6f_\pi ^2)\right)^2}\frac{T}{6f_\pi ^2}=s(\tau _0)\tau _0$$
(20)
where now $`<A^2>=<A^2>(1T_0^2/(6f_\pi ^2))`$. The numerical results are presented on figure 1. We can observe that the correction, which induce a faster cooling of the system is of order 5%.

Even if the approach can be considered too simplistic, it is not reasonably expected that an improved treatment can modify noticeably the order of magnitude. One can thus conclude that the effect due the spontaneous breaking of chiral symmetry is quite negligible and without any experimental consequences.

## V Dissipative effects

In the previous part we have seen that it is possible to neglect the effect of the breaking of the symmetry on the hydrodynamics theory itself, that is we can use ”normal” hydrodynamics. Nevertheless an important question which has already been addressed in murongaWinter  is the effect of dissipation on top of this model of perfect fluid.

There are two non equivalent possibilities to take into account dissipation. First, one can use the usual formalism and just say that when an inhomogeneity appears in a system, the latter reacts and a flow proportional to the gradient develops itself in order to restore equilibrium. This is first order theory. However it is well-known that this approach presents instabilities and that it is necessary to consider a relaxation time for the flow. This goal is achieved through Müller-Israel-Stewart second order theory. In the following we will explore both methods.

In murongaWinter , the author showed that only the dissipative effects related to the shear viscosity $`\eta `$ play a role for the temperature profile in the (1+1) dimensional Bjorken model. The main physical reason is that bulk viscosity vanishes for an ultra relativistic gas. He took the shear viscosity from Welke  where the temperature dependence for $`\eta `$ is $`\eta =\pi f_\pi ^4/(8T)`$. Then, the differential equations to be solved for the energy density $`ϵ=3aT^4`$ are:

$$\frac{ϵ}{\tau }=\frac{4}{3}\frac{ϵ}{\tau }+\frac{\mathrm{\Phi }}{\tau }$$
(21)
where $`\mathrm{\Phi }`$ is related to the stress tensor (see murongaWinter  for further details) and $`\mathrm{\Phi }=0`$ for the ideal (perfect) fluid, $`\mathrm{\Phi }=\frac{4}{3}\frac{\eta }{\tau }`$ for the first order and $`\tau _\pi \frac{\mathrm{\Phi }}{\tau }=\mathrm{\Phi }+\frac{4}{3}\frac{\eta }{\tau }`$ for the second order theory discussed before. In the last equation $`\tau _\pi `$ is the relaxation time associated to $`\mathrm{\Phi }`$.

We now use (21), using the same inputs as murongaWinter , namely massless pions and current algebra cross section. We can see on figure 2 that the correction due to viscosity compared to the ideal case is totally negligible not only for the (non-physical but usual) first order but also for the second order. Our results are different from those presented in the original paper murongaWinter . The difference lies in fact in a small computational mistake in murongaWinter . The main point is that, contrary to what is claimed in murongaWinter , there are, in this simple approach, no appreciable effects. If we would stop here, we could say, as in the first part dedicated to the breaking of chiral symmetry that a realistic description of a pion gas can be achieved with the perfect fluid model.

However, we stress that this conclusion only holds because of the hypothesis used in murongaWinter . The temperature dependences of the viscosity $`\eta `$ and of the relaxation time $`\tau _\pi `$ are only valid for massless pions and when the cross-section for $`\pi \pi `$ interaction is given by current algebra. If we take instead the realistic case of massive pions and the experimental cross-section, these dependences are completely different. For instance, $`\eta `$ increases with temperature Welke ; moi  instead of decreasing. With this realistic $`T`$-dependence for the viscosity and the relaxation time for the pions taken from moi , we obtain the curves plotted on figure 2. We can now observe that the first and the second order are still close but they are now around 20% above the ideal fluid case. This means that dissipative effects imply a slower cooling of the system and should be incorporated in the simulations of heavy ion collisions.

## VI Conclusion

In the first part of this paper we have examined in a simple model an application of hydrodynamic equations with spontaneous breaking symmetry. In order to get an estimate of the magnitude of the effects, we have used the one dimensional Bjorken description of a heavy ion collision. We have shown that the model of the perfect fluid is almost correct : the effect of the breaking of the symmetry is to give a correction of about 1% when no $`T`$-dependence of the pion decay constant is taken into account and a few percents when we consider a crude $`T`$-dependence. In a second part we have considered an other deviation from the perfect fluid model, namely the dissipation. The dissipative effects, when taken into account in a simple description murongaWinter  are not so important as claimed in this reference but represent actually a few percents. On the contrary, we have obtained sizable effects with a correct $`T`$-dependence for the shear viscosity and the shear flow relaxation time. For example, for an initial temperature of 200 MeV, the temperature goes from about 110 MeV in the ideal case to 140 MeV with our dissipative effects for $`\tau /\tau _06`$. The whole conclusion of this paper is that the temperature profile generated from the perfect fluid approach and used in many simulations is almost correct if one takes into account only the effects of symmetry breaking. On the contrary, if we consider dissipative effects the profile can be substantially modified, that is we have a slower cooling of the pion gas.

Acknowledgments : We thank M. Ericson, Emeritus Professor at University Claude Bernard and G. Chanfray for constant help and enlightening discussions. We also thank J. Serreau for many discussions.


text:
# Separability criterion for pure states in multipartite and high dimensional systems

## Abstract

We propose a sufficient and necessary separability criterion for pure states in multipartite and high dimensional systems. Its main advantage is operational and computable. The obvious expressions of this criterion can be given out by the coefficients of components of the pure state. In the end, we simply mention a principle method how to define and obtain the measures of entanglement in multipartite and high dimensional systems.

Entanglement is one of the most striking feature of quantum theory. It is very important how to quantify it. In our point of view, separability is a theoretical foot stone to define the measures of entanglement.

A $`n`$-partite quantum state $`\rho _{A_1A_2\mathrm{}A_n}`$ is called (fully) separable iff it can be written as a convex combination of product states, i.e.

$$\rho _{A_1A_2\mathrm{}A_n}=\underset{m}{}p_m\rho _{A_1}^{(m)}\rho _{A_2}^{(m)}\mathrm{}\rho _{A_n}^{(m)}$$
(1)
where any $`\rho _{A_i}^{(m)}(i=1,2,\mathrm{},n)`$ is a density matrix of $`r_i`$-dimensional. Otherwise, $`\rho _{A_1A_2\mathrm{}A_n}`$ is entangled. For low dimensional ($`2\times 2`$ and $`2\times 3`$) of bipartite systems, Peres and Horodecki gave out a necessary and sufficient condition for separability of arbitrary states, that is, $`\rho _{A_1A_2}`$ is separable if and only if it has a “positive partial transpose” (PPT) Peres ; Horodecki1 . Unfortunately, this criterion is only necessary for a sate to be separable in higher dimensions, but is not sufficient. Recently a lot of works about separability have been presented Lewenstein ; Rudolph ; Chen ; Horodecki2 ; Doherty ; Giedke ; Karnas ; Dur ; Acin ; My1 . These works explain a fact that the operational and computable separability criterion, in particular, in multipartite and high dimensional systems becomes very interesting now. Here, our letter will devote to this problem.

First of all, we need the following lemma:

Lemma If a pure state $`\rho _{A_1A_2\mathrm{}A_n}^\mathrm{P}`$ is separable, thus it must be “purely” separable, that is

$$\rho _{A_1A_2\mathrm{}A_n}^\mathrm{P}=|\psi _{A_1A_2\mathrm{}A_n}\psi _{A_1A_2\mathrm{}A_n}|=\underset{i=1,}{\overset{n}{}}\rho _{A_i}$$
(2)
Proof Sine $`\rho _{A_1A_2\mathrm{}A_n}^\mathrm{P}`$ is separable, it must be able to write as the form of eq.(1). Obviously, for a given $`m`$, $`\rho _{A_1}^{(m)}\rho _{A_2}^{(m)}\mathrm{}\rho _{A_n}^{(m)}`$, as a Hermit operator in Hilbert space $`_{A_1}_{A_2}\mathrm{}_{A_n}`$, has a complete set of orthogonal and normalized eigenvectors, which is denoted by $`\{|\psi _{A_1A_2\mathrm{}A_n}^m(\alpha ),(\alpha =1,2,\mathrm{},N)\}`$. When a value of $`m`$ is chosen, we can take them as a set of basis of Hilbert space $`_{A_1}_{A_2}\mathrm{}_{A_n}`$, it implies we can write

$$|\psi _{A_1A_2\mathrm{}A_n}=\underset{\alpha =1}{\overset{N}{}}C_\alpha ^m|\psi _{A_1A_2\mathrm{}A_n}^m(\alpha )$$
(3)
Note that a $`m`$-th component state $`\rho _{A_1}^{(m)}\rho _{A_2}^{(m)}\mathrm{}\rho _{A_n}^{(m)}`$ has been taken as a pure state, without loss of generality, we set $`_{i=1,}^n\rho _{A_i}^{(m)}=|\psi _{A_1A_2\mathrm{}A_n}^m(1)\psi _{A_1A_2\mathrm{}A_n}^m(1)|`$. Thus,

$`1`$ $`=`$ $`\psi _{A_1A_2\mathrm{}A_n}|\rho _{A_1A_2\mathrm{}A_n}^\mathrm{P}|\psi _{A_1A_2\mathrm{}A_n}`$ (4)
$`=`$ $`{\displaystyle \underset{m}{}}p_m\psi _{A_1A_2\mathrm{}A_n}|{\displaystyle \underset{i=1,}{\overset{n}{}}}\rho _{A_i}^{(m)}|\psi _{A_1A_2\mathrm{}A_n}`$
$`=`$ $`{\displaystyle \underset{m,\alpha ,\alpha ^{}}{}}p_mC_\alpha ^mC_\alpha ^{}^m\psi _{A_1A_2\mathrm{}A_n}^m(\alpha )|{\displaystyle \underset{i=1,}{\overset{n}{}}}\rho _{A_i}^{(m)}|\psi _{A_1A_2\mathrm{}A_n}^m(\alpha ^{})`$
$`=`$ $`{\displaystyle \underset{m}{}}p_m{\displaystyle \underset{\alpha ,\alpha ^{}}{}}C_\alpha ^mC_\alpha ^{}^m\delta _{\alpha \alpha ^{}}\delta _{\alpha ^{}1}={\displaystyle \underset{m}{}}p_mC_1^m{}_{}{}^{}C_{1}^{m}`$
It implies that $`C_1^m{}_{}{}^{}C_{1}^{m}=1`$ and then the other $`C_\alpha ^m(\alpha =2,3,\mathrm{},N)`$ are zero. Furthermore, we have $`\rho _{A_1}^{(m)}\rho _{A_2}^{(m)}\mathrm{}\rho _{A_n}^{(m)}=\rho _{A_1A_2\mathrm{}A_n}^\mathrm{P}`$, that is, this lemma is valid.

Without loss of generality, the density matrix $`\rho _{A_1A_2\mathrm{}A_n}`$ in a $`n`$-partite system can be expanded as

$`\rho _{A_1A_2\mathrm{}A_n}`$ $`=`$ $`{\displaystyle \frac{1}{2^n}}{\displaystyle \underset{\mu _1=0}{\overset{r_1^21}{}}}{\displaystyle \underset{\mu _2=0}{\overset{r_2^21}{}}}\mathrm{}{\displaystyle \underset{\mu _n=0}{\overset{r_n^21}{}}}`$ (5)
$`a_{\mu _1\mu _2\mathrm{}\mu _n}\lambda _{A_1}^{\mu _1}\lambda _{A_2}^{\mu _2}\mathrm{}\lambda _{A_n}^{\mu _n}`$
where $`a_{\mu _1\mu _2\mathrm{}\mu _n}`$ is a $`n`$-rank real tensor, $`\lambda _{A_i}^0`$ is proportional to $`r_i`$-dimensional identity matrix and its proportional factor can be taken as $`2/r_i`$ in usual. While $`\lambda _{A_i}^{m_i}(m_i=1,2,\mathrm{},r_i^21)`$ are generators of $`SU(r_i)`$.

For simplicity, we will limit our focus on such a multipartite system in which every partite only has a qubit or/and a qutrit, and then the corresponding generators are respectively Pauli matrices $`\sigma _i(i=1,2,3)`$ and Gell-Mann matrices $`\lambda _i(i=1,2,3,4,5,6,7,8)`$.

Theorem One. The sufficient and necessary conditions for separability of an arbitrary pure state in $`n`$-partite systems made up of $`n`$ qubits (every partite has a qubit) are

$$𝝃_{A_1}^2=𝝃_{A_2}^2=\mathrm{}=𝝃_{A_n}^2=1$$
(6)
where $`𝝃_{A_i}`$ is a polarized vector of the reduced density matrix $`\rho _{A_i}=\mathrm{Tr}_{_{j=1,ji}^nA_j}(\rho _{A_1A_2\mathrm{}A_n})=\frac{1}{2}(\sigma _0+𝝃_{A_i}𝝈)`$ (tracing off $`n1`$ partites).

Proof In order to prove it, let us start from the well-known bipartite systems. If denoting the pure state in a bipartite system as

$$|\psi _{A_1A_2}=a|00+b|01+c|10+d|11$$
(7)
we can obtain

$$𝝃_{A_1}^2=𝝃_{A_2}^2=\underset{i=1}{\overset{3}{}}a_{i0}a_{i0}=\underset{i=1}{\overset{3}{}}a_{0i}a_{0i}=14|adbc|^2$$
(8)
Thus, the sufficient and necessary condition for separability of an arbitrary pure state in bipartite systems of two qubits can be written in form

$$ad=bc$$
(9)
It is easy to prove that theorem one is a necessary creiterion. In fact, if the pure state $`|\psi _{A_1A_2}`$ is separable, in terms of our lemma, without loss of generality we have

$$|\psi _{A_1A_2}=(a_1|0+b_1|1)(a_2|0+b_2|1)$$
(10)
Comparing it with the equation (7), up to an undetermined overall phase factor, we can obtain

$$a=a_1a_2,b=a_1b_2,c=a_2b_1,d=b_1b_2$$
(11)
Substituting eq.(11) into eq.(8), it immediately following that our separability criterion is satisfied.

Now let us prove it to be sufficient. If there is only one non-zero element in the coefficient set $`\{a,b,c,d\}`$ of the concerning pure state (7) then this state is obviously separable. Actually, at least, there are two non-zero element in the coefficient set, the case is not trivial. It is easy to obtain the product states in all of possible cases as following

(1) $`a,b0`$

$`|\psi _{A_1}`$ $`=`$ $`{\displaystyle \frac{|a|}{\sqrt{|a|^2+|c|^2}}}\left(|0+{\displaystyle \frac{c}{a}}|1\right)`$ (12a)
$`|\psi _{A_2}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{|a|^2+|b|^2}}}(a|0+b|1)`$ (12b)
(2) $`a,c0`$

$`|\psi _{A_1}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{|a|^2+|c|^2}}}\left(a|0+c|1\right)`$ (13a)
$`|\psi _{A_2}`$ $`=`$ $`{\displaystyle \frac{|c|}{\sqrt{|c|^2+|d|^2}}}\left(|0+{\displaystyle \frac{d}{c}}|1\right)`$ (13b)
(3) $`a,d0`$ or $`b,c0`$ , thus $`b,c`$ or $`a,d`$ must be non-zero. Product states can be taken the same form as above.

(4) $`b,d0`$

$`|\psi _{A_1}`$ $`=`$ $`{\displaystyle \frac{|b|}{\sqrt{|b|^2+|d|^2}}}\left(|0+{\displaystyle \frac{d}{b}}|1\right)`$ (14a)
$`|\psi _{A_2}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{|a|^2+|b|^2}}}(a|0+b|1)`$ (14b)
(5) $`c,d0`$, it is the same as the case (2).

It is simple to verify that for all of the cases

$$|\psi _{A_1A_2}=|\psi _{A_1}|\psi _{A_2}$$
(15)
by means of two conditions. One is separable condition $`ad=bc`$, another is the normalized condition $`aa^{}+bb^{}+cc^{}+dd^{}=1`$. So, theorem one is a sufficient criterion.

Furthermore, we can use mathematical induction to prove the case of $`n`$-partite systems. Suppose this separability criterion is sufficient and necessary for $`n1`$ partite systems. In the case of $`n`$-partite systems, by tracing off any single partite, a separable pure state in a $`n1`$ partite system is obtained. Thus, the norm square of polarized vector of each partite for this state is always 1. Since the partite traced off is arbitrary, the criterion for $`n`$-partite systems must be valid. Theorem one is then necessary. Have assumed that for any $`n1`$ partite the forms of product states can be written out. If we can not obtain the forms of product states for $`n`$ parties system which obey our criterion, it must be conflict with our precondition. Therefore, our criterion is sufficient. As an example , let us verify the case of a three partite system. Denoting

$`|\psi _{A_1A_2A_3}`$ $`=`$ $`a|000+b|001+c|010+d|011`$ (16)
$`+e|100+f|101+g|110+h|111`$
we can evaluate out

$`𝝃_{A_1}^2`$ $`=`$ $`14|afbe|^24|agce|^24|ahde|^2`$ (17a)
$`4|bgcf|^24|bhdf|^24|chdg|^2`$
$`𝝃_{A_2}^2`$ $`=`$ $`14|adbc|^24|agce|^24|ahcf|^2`$ (17b)
$`4|bgde|^24|bhdf|^24|ehfg|^2`$
$`𝝃_{A_3}^2`$ $`=`$ $`14|adbc|^24|afbe|^24|ahbg|^2`$ (17c)
$`4|cfde|^24|chdg|^24|ehfg|^2`$
Thus, the obviously independent expressions of our separability criterion in terms of coefficients of component of pure state are

$`ad`$ $`=`$ $`bc,af=be,ah=bg,ag=ce,bg=cf,`$
$`bh`$ $`=`$ $`df,cf=de,ch=dg,eh=fg`$ (18)
If $`|\psi _{A_1A_2A_3}`$ is separable, without loss of generality, in terms of our lemma we have

$$|\psi _{A_1A_2A_3}=(a_1|0+b_1|1)(a_2|0+b_2|1)(a_3|0+b_3|1)$$
(19)
Comparing with eq.(16), it follows that

$`a=a_1a_2a_3,b=a_1a_2b_3,c=a_1b_2a_3,d=a_1b_2b_3`$
$`e=b_1a_2a_3,f=b_1a_2b_3,g=b_1b_2a_3,h=b_2b_2b_3`$ (20)
Therefore, we have $`𝝃_{A_1}^2=𝝃_{A_2}^2=𝝃_{A_3}^2=1`$. When $`a,b0`$, we obtain

$`|\psi _{A_1A_2}`$ $`=`$ $`{\displaystyle \frac{|a|}{\sqrt{|a|^2+|c|^2+|e|^2+|g|^2}}}`$ (21a)
$`\left(|00+{\displaystyle \frac{c}{a}}|01+{\displaystyle \frac{e}{a}}|10+{\displaystyle \frac{g}{a}}|11\right)`$
$`|\psi _{A_3}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{|a|^2+|b|^2}}}(a|0+b|1)`$ (21b)
Again, when the coefficients $`c,e,g`$ are all zero,

$$|\psi _{A_1}=|0|\psi _{A_2}=|0$$
(22)
and if at least one element among the set $`\{c,e,g\}`$ is not zero, thus we can write down the product states like two partite case. For example $`c0`$, we have

$`|\psi _{A_1}={\displaystyle \frac{|c|}{\sqrt{|c|^2+|g|^2}}}\left(|0+{\displaystyle \frac{g}{c}}|1\right)`$ (23a)
$`|\psi _{A_2}={\displaystyle \frac{|a|}{\sqrt{|a|^2+|c|^2}}}\left(|0+{\displaystyle \frac{c}{a}}|1\right)`$ (23b)
For the other cases, we can write down the product states in a similar way. In particular, if only one partite is separable with the other two partite, then the corresponding separable partite has a normalized polarized vector. For example, $`A_3`$-partite is separable with $`A_1A_2`$ partites, then $`𝝃_{A_3}^2=1`$. It must be emphasized that $`A_2`$ partite is separable with $`A_1A_3`$ partites means

$$(I_{A_1}S_{A_2A_3})\rho _{A_1A_2A_3}(I_{A_1}S_{A_2A_3}^{})=\rho _{A_1A_3}\rho _{A_2}$$
(24)
where $`S`$ is an exchanging operator for the near partites,

$$S=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right)$$
(25)
Therefore, we think that the partially separability, which means that some partites are separable with the other partites, has the following criterion:

Conjecture One The sufficient and necessary condition for partially separability of an arbitrary pure state in multipartite systems are that the polarized vectors of separable parties obey

$$𝝃_i^2=1,\text{ }A_i\text{ takes over all of separable partites}$$
(26)
We have checked it up to four partite systems, but the partially separable problem is still open to prove it for arbitrary partite systems in a strict method.

Theorem Two The sufficient and necessary conditions for separability of an arbitrary pure state in $`n`$-partite systems made up of $`n`$ qutrits (every partite has a qutrit) are

$$𝝃_{A_1}^2=𝝃_{A_2}^2=\mathrm{}=𝝃_{A_n}^2=\frac{4}{3}$$
(27)
where $`𝝃_{A_i}`$ is a coherence vector of the reduced density matrix $`\rho _{A_i}=\mathrm{Tr}_{_{j=1,ji}^nA_j}(\rho _{A_1A_2\mathrm{}A_n})`$ (tracing off $`n1`$ partite), that is, its component is

$$\xi _{A_i}^m=\left(\frac{\sqrt{6}}{2}\right)^{n1}a_{0_10_2\mathrm{}0_{i1}m0_{i+1}\mathrm{}0_n}(m=1,2,\mathrm{},8)$$
(28)
where subindex $`0_j`$ indicates 0 in the position $`j`$ and we have used the expansion of the density matrix of a qutrit

$$\rho _{A_i}=\frac{1}{3}I_{3\times 3}+\frac{1}{2}\underset{m=1}{\overset{8}{}}\xi _{A_i}^m\lambda _m$$
(29)
Proof Let us start from bipartite systems. If denoting the pure state in a bipartite system as

$`|\psi _{A_1A_2}`$ $`=`$ $`x_1|11+x_2|12+x_3|13+x_4|21+x_5|22`$ (30)
$`+x_6|23+x_7|31+x_8|32+x_9|33`$
we can obtain

$`𝝃_{A_1}^2`$ $`=`$ $`𝝃_{A_2}^2=4/34|x_1x_5x_2x_4|^2`$ (31)
$`4|x_1x_6x_3x_4|^24|x_1x_8x_2x_7|^2`$
$`4|x_1x_9x_3x_7|^24|x_2x_6x_3x_5|^2`$
$`4|x_2x_9x_3x_8|^24|x_4x_8x_5x_7|^2`$
$`4|x_4x_9x_6x_7|^24|x_5x_9x_6x_8|^2`$
Thus, the sufficient and necessary conditions for separability of an arbitrary pure state in bipartite systems of two qutrits can be written in form

$`0`$ $`=`$ $`x_1x_5x_2x_4=x_1x_6x_3x_4=x_1x_8x_2x_7`$ (32)
$`=`$ $`x_1x_9x_3x_7=x_2x_6x_3x_5=x_2x_9x_3x_8`$
$`=`$ $`x_4x_8x_5x_7=x_4x_9x_6x_7=x_5x_9x_6x_8`$
It is easy to prove theorem two to be necessary. In fact, if the pure state $`|\psi _{A_1A_2}`$ is separable, in terms of our lemma, without of generality we have

$$|\psi _{A_1A_2}=(a_1|1+b_1|2+c_1|3)(a_2|1+b_2|2+c_2|3)$$
(33)
Comparing it with the equation (42), up to an undetermined overall phase factor, we can obtain

$`x_1=a_1a_2,x_2=a_1b_2,x_3=a_1c_2,`$ (34)
$`x_4=b_1a_2,x_5=b_1b_2,x_6=b_1c_2,`$ (35)
$`x_7=c_1a_2,x_8=c_1b_2,x_9=c_1c_2,`$ (36)
Substituting eq.(34) into eq.(31), it immediately following that our separable criterion is satisfied.

It is a trivial case if only if there is one non-zero element in the coefficient set $`\{x_i,i=1,2,\mathrm{},8\}`$. Consider the case that at least two coefficients are not zero, for example $`x_1,x_20`$, we have: when $`x_3,x_6,x_9=0`$

$`|\psi _{A_1}`$ $`=`$ $`{\displaystyle \frac{|x_1|}{\sqrt{|x_1|^2+|x_4|^2+|x_7|^2}}}`$ (37)
$`\times \left(|1+{\displaystyle \frac{x_4}{x_1}}|2+{\displaystyle \frac{x_7}{x_1}}|3\right)`$
$`|\psi _{A_2}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{|x_1|^2+|x_2|^2}}}\left(x_1|1+x_2|2\right)`$ (38)
when at least one element in the set $`x_3,x_6,x_9`$ is not zero, for example $`x_30`$, $`|\psi _{A_1}`$ can be taken the same form as above, but

$$|\psi _{A_2}=\frac{1}{\sqrt{|x_1|^2+|x_2|^2+|x_3|^2}}\left(x_1|1+x_2|2+x_3|3\right)$$
(39)
It must be emphasized that since separability conditions give out the relation among the coefficients, the expressions of product states can be written as the different forms. Likewise, we can prove the other cases. In order to save space, we omit them. So, theorem two is a suffcient criterion. As to the proof for $`n`$-partite systems, it is similar to one of theorem one. We do not intend to repeat it. We also have checked the three partite systems made up of three qutrits. Except for the computation detail, there is no any new idea comparing with the three partite systems made up of three qubits.

In general, we have

Conjecture Two In $`n`$-partite systems and the $`i`$th-partite with $`r_i`$ dimensional, the sufficient and necessary conditions for separability of an arbitrary pure state are

$$𝝃_{A_i}^2=2\left(1\frac{1}{r_i}\right),(i=1,2,\mathrm{},n)$$
(40)
where $`𝝃_{A_i}`$ is a coherence vector of the reduced density matrix $`\rho _{A_i}=\mathrm{Tr}_{_{j=1,ji}^nA_j}(\rho _{A_1A_2\mathrm{}A_n})`$ (tracing off $`n1`$ partite) defined by

$$\rho _{A_i}=\frac{1}{r_i}I_{r_i\times r_i}+\frac{1}{2}\underset{m_i=1}{\overset{r_i^21}{}}𝝃_{A_i}^{m_i}\lambda _{m_i}$$
(41)
$`I_{r_i\times r_i}`$ is $`r_i`$-dimensional identity matrix and $`\lambda _{A_i}^{m_i}(m_i=1,2,\mathrm{},r_i^21)`$ are generators of $`SU(r_i)`$.

By using of our above idea and method, we can discuss the multipartite systems made up of qudits as well as the multipartite systems made up of mixture of qubits, qutrits and qudits. We can find our separability criterion to be valid. For example, $`A_1`$ partite is a qubit and $`A_2`$ partite is a qutrit. If denoting

$$|\psi _{A_1A_2}=x_1|01+x_2|02+x_3|03+x_4|11+x_5|12+x_6|13$$
(42)
It is easy to obtain

$`1𝝃_{A_1}^2`$ $`=`$ $`{\displaystyle \frac{4}{3}}𝝃_{A_2}^2=4|x_1x_5x_2x_4|^2`$ (43)
$`+4|x_1x_6x_3x_4|^2+4|x_2x_6x_3x_5|^2`$
Of course, in the same way we can prove the above conjecture. However, how to prove our conjecture two strictly and generally is still open.

In the end, we would like to point out that our result is helpful for defining measures of entanglement in multipartite and high dimensional systems. Actually, we can verify that for any cat states in multipartite systems, the norm of coherent vector of reduced density matrix of any $`A_i`$-partite takes the minimum value. Moreover, we have proved that for a (fully) separable pure state, the norm of coherent vector of reduced density matrix of any $`A_i`$-partite arrives at the maximum value. Therefore, the entropy functions or their combinations or the other monotonic continuous functions varying with the norms of coherent vectors of reduced density matrices of all of partites or their functions can be defined as the measures of entanglement for pure state in multipartite and high dimensional systems. Just like done and doing by us (limit within a scalar description and including measure of strength of entanglement degree of all of concerned partites) My2 ; My3 ; My4  based on Bennett, Wootters et.al’s as well as Vedral’s ideas Bennett ; Wootters ; Vedral  . Perhaps, not only the norm but also components of coherent vectors in a tensor (including vector) description of entanglement are required. While, distillation and preparation of entanglement for multipartite systems will act a role to judge what form is more reasonable and applicable.