format: text

text:
# The effect of synchronized area on SOC behavior in a kind of Neural Network Model 1footnote 11footnote 1The project supported by P.R.China National Basic Research Project ”Nonlinear Science”

## I Introduction

A few years ago, Bak, Tang and Wiesenfeld introduced the concept of the ”Self-Organized Criticality” (SOC) in sand pile model 1 . From then on, this concept has been widely studied in some extended dissipative dynamical systems, such as earthquake 2 , biology evolution 3 , and so on. It is shown that all these systems can naturally evolve into a ”critical state” with no intrinsic spatial and temporal scales through a self-organized process without the need to fine-tune parameters of the system. This critical state is characterized by a power-law distribution of avalanche sizes, where the size is the total number of toppling events or unstable units.

Now,the research on SOC has come into a new level. One has studied several factors’ influence on SOC behavior, such as network size, periodic or nonperiodic boundary conditions,local dynamics variable is conservative or not, and so on. Many investigators believe that it is intrinsic-stability(order) and variability(disorder)’s common action make the system evolves into a ”frozen disorder” (SOC) state 4 ; 5 ; 6 . It is the combination feature of stability and variability, and it’s complex spatial-temporal dynamic behavior, make the system in SOC state have maximum complexity and latent computing potency.

The brain is a complex system and its information process has the properties of stability and variability – on one hand, there are relative stable information stored mechanism, and stored area in brain(such as feature area in cortex) ; on the other hand, the brain is influenced by the environment and one should continuous update knowledge and concepts. The similarity between the SOC systems and the brain has lead us to study Artifical Neural Network(ANN) and SOC together. There is some SOC behavior shown in the neuron network model introduced by our group 7 .

The brain is also a complex system with highly complexity, highly order and special structure. The structure must have the effect on the brain’s dynamics behavior. Our neuron network can also produce some special structure, so we believe that the SOC behavior shown in our model must have it’s own special behavior and rule.

In this paper, the feature area(produced by self-organized process)’s definite effect on SOC behavior has been studied.

## II Model

Here we propose a two-dimensional neural network model of square lattice. This model is a kind of serial self-organized neural network model, based on the LISSOM model 8 . When a $`h`$\- dimensional vector $`\zeta `$ is inputted, the state $`\eta _{ij}(t)`$ of the neuron $`(i,j)`$ at time $`t`$ is changed according to the formula:

$`\eta _{ij}`$ $`=`$ $`\sigma \{{\displaystyle \underset{h}{}}\mu _{ij,h}\zeta _h+\gamma _e{\displaystyle \underset{kl}{}}E_{ij,kl}\eta _{kl}(t1)`$ (1)
$`\gamma _i{\displaystyle \underset{k^{}l^{}}{}}I_{ij,k^{}l^{}}\eta _{k^{}l^{}}(t1)\}`$
$`=`$ $`\sigma \{f_{ij}(t)\}`$
where $`\sigma (x)`$ is an active function, and we design it as sign function, i.e., if $`x0`$, then $`\sigma (x)=1`$ , otherwise $`\sigma (x)=1`$ . $`\mu _{ij,h}`$ is an afferent input weight vector; $`E_{ij,kl}`$ is the excitatory lateral connection weight on the connection from the neuron $`(k,l)`$ to the $`(i,j)`$ neuron; $`I(ij,kl)`$ is the inhibitory connection weight. $`f_{ij}(t)`$ is the local field of the neuron $`(i,j)`$ at time $`t`$. $`\gamma _e`$ and $`\gamma _i`$ are constant factors. The adjustment of those three connection weights is as following according to the dynamic Hebb rule:

$`\mu _{ij,h}(t+1)={\displaystyle \frac{\mu _{ij,h}(t)+\alpha \eta _{ij}\xi _h}{\{_h\left[\mu _{ij,h}(t)+\alpha \eta _{ij}\xi _h\right]^2\}^{1/2}}}\text{,}`$
$`E_{ij,kl}(t+1)={\displaystyle \frac{E_{ij,h}(t)+\alpha _E\eta _{ij}\eta _{kl}}{_{ij}\left[E_{ij,kl}(t)+\alpha _E\eta _{ij}\eta _{kl}\right]}}\text{,}`$
$`I_{ij,k^{}l^{}}(t+1)={\displaystyle \frac{I_{ij,k^{}l^{}}(t)+\alpha _I\eta _{ij}\eta _{k^{}l^{}}}{_{k^{}l^{}}\left[I_{ij,k^{}l^{}}(t)+\alpha _I\eta _{ij}\eta _{k^{}l^{}}\right]}}\text{.}`$ (2)
where $`\alpha ,\alpha _E,\alpha _I`$ are the learning rates.

Note the change of the neuron states is quick and the adjustment of the connection weights is slow. Usually, after over 10 iterations of the neuron states when any pattern is inputted (at this time, the network state becomes an attractor in the state space, often a fixed point.), all connection weights are updated once. Thus, after learning a while, the lateral connection weight self-evolves into the”Mexican hat” profile,the afferent input weight self-organizes into a topological map of the input space, the neuron network can produce some special feature areas, and the state of the neural network is evolved from disordered case to stable and topological case in state space.

Then we introduce the following interactive process between the neurons, similar with the pulse coupled interaction:

1) When the neuron $`(i,j)`$ is stable, i.e., $`\eta _{ij}(t)=\sigma \{f_{ij}(t1)\}`$ , it doesn’t influence the others;

2) If the neuron $`(i,j)`$ is unstable, i.e., $`\eta _{ij}\sigma \{f_{ij}(t1)\}`$ or $`f_{ij}(t1)=0`$ , then the nearest neighbors $`(i^{},j^{})`$ around this unstable neuron will receive a pulse respectively and their local fields will be changed. At the same time, the neuron $`(i,j)`$ becomes stable again, depending on the formula (1);

3) When all neurons of the neural network are stable, we choose the minimum $`g`$ among the absolute values of all local fields $`f_{ij}`$ and drive the local field of every neuron, i.e.,

$$f_{ij}f_{ij}c\eta _{ij}g$$
(3)
where $`c`$ is a constant.

Now, we present the computer simulation procedure of this model in detail:

1) Variable initialization. In the 2-dimensional $`n\times n`$ neural network model, let the initiatory state and local field equal $`0`$; random initialize each connection weight among $`[1,1]`$ ; and produce $`M`$ random input patterns, $`\zeta _j^i[1,1]`$, $`(i=1,2,\mathrm{}M;j=1,2,\mathrm{}h)`$.

2) Learning process. According to formula (1), we input the pattern circularly and iterate the neuron state and local field. After period of time, the space state of the network reaches stability and we consider the $`M`$ input patterns have been stored.

3) Associative memory. Input a new pattern, and then search the unstable neuron $`(i,j)`$ as defined above in whole neural network. Due to being unstable, the neuron $`(i,j)`$ discharges a pulse to the each nearest neighbor $`(i^{},j^{})`$ and thus causes the local fields of them to change as following:

$$f_{i^{}j^{}}f_{i^{}j^{}}\frac{\gamma }{2}\eta _{i^{}j^{}}(1+|f_{ij}|)$$
(4)
where $`\gamma `$ represents the pulse intensity, symbol $`||`$ denotes absolute value.

Simultaneously, according to formula (1), the neuron $`(i,j)`$ becomes stable as: $`\eta _{ij}\sigma (f_{ij})`$ , $`f_{ij}\eta _{ij}`$ ; where $`\sigma (x)`$ is sign function.

Repeat this procedure until all neurons of the model are stable. Define one avalanche as all unstable neurons in this process. Then begin drive process by formula (3) and new avalanche.

## III Simulation results

Recently, Bak and Sneppen have investigated the power law distribution $`P(X)`$ of the distances $`X`$ between subsequent unstable sites in lattice of BS biology evolution model 3 , and J.De.Boer et al. have found some interested result with it 9 .So in this paper we studied not only the distribution $`P(S)`$ of the avalanche sizes $`S`$ but also the distribution $`P(X)`$ of the distances $`X`$ between the subsequent unstable sites . We find the distribution will deviate from the power law in some conditions.

### III.1 The effect of synchronized area on SOC behavior

The size of our lattice is $`40\times 40`$.We find that in associative memory process, the distribution of the avalanche sizes has power-law behavior, $`P(S)S^\tau `$ , $`\tau 0.90`$ . It is shown in Fig.1.a. The distribution $`P(X)`$ of the distances $`X`$ between the subsequent unstable sites has power-law behavior too, $`P(X)X^\beta `$ , $`\beta 2.23`$ , it is shown in Fig.1.b.

The SOC behavior changed with the scope of lateral connection has been studied. We increase the excitatory lateral connection radius $`d_e`$ and the inhibitory connection radius $`d_i=3d_e`$. It can be seen in Fig 1.a that avalanche size and occurring probability of large scale avalanche decrease with the slope $`d_e`$ increasing, and the distribution of $`P(X)`$ deviates from power law more and more in large $`X`$, the probability of large distance $`X`$ between the subsequent unstable sites also increases, we consider it is a deviation from SOC behavior, it is shown in Fig 1.b.

By investigating , we think our system is in a ”partly-synchronized” state, hence the dynamics behavior mentioned above can be seen.

A.Corral et al. propose SOC state and synchronization state might be considered as two uttermost state of system (just like two sides of the same coin) 4 . The inhomogeneity introduced by boundary or initialization conditions can propagate into interior of network, hence makes the system evolve into SOC state 10 ; 11 . When inhomogeneity is not large enough, the system finds a compromise between synchronization and SOC 6 . It can be considered as a partly-synchronized state.

If Our Model is only a pure OFC model without learning process, the system will present a macroscopic SOC behavior among almost all the lattices, but it is also a neuron network model. As a kind of self-organized feature map model, after learning a while, it’s neurons will develop a unique lateral interaction ”Mexican hat” profile that represents its long-term associations with each other. the afferent input weights will self-organize into a topological map of the input space 8 , it can form some special topological feature regions. We consider that in these regions, as the connection weights become more topographically ordered, neuron’s synchronization effect between each other will be reinforced. When order in these regions is applied to the model, the system has a tendency from SOC state to synchronized state. At last, the system finds a compromise between synchronization and SOC, it could be seen as ”partly-synchronized” state. With the process, the distribution of avalanche varies from the continuous distribution to a discrete one, the possibility of large scale avalanche propagating into the interior of the synchronization regions will reduce, and occurring probability of large scale avalanche will decrease too, the one-off isolated avalanche(only one unstable site in an avalanche)in these regions will increase greatly, it makes the distances $`X`$ between the subsequent unstable sites have a stochastic spatial even distribution in the area. This distribution has more effect on probability of large distance between unstable sites than probability of small distance. (Because probability of small distance is larger than one of large distance.)

So at this moment, there are some areas in synchronized state and another regions in SOC state. We can consider that with the increasing of $`d_e`$, the feature region(synchronized region) produced by self-organized process becomes wider, which results in the whole system dynamics deviating largely from SOC state, it can be seen in Fig.1.

To verify the idea mentioned above, we draw the avalanche’s distribution map of the whole system. We draw one avalanche’s distribution snapshot every 1000 avalanches, then overlap all the snapshots in one picture. The result can be seen in Fig2. It can be clearly seen that the blank region (seldom avalanche area) expands with the increasing of $`d_e`$. It means that the synchronized region introduced by self-organized expands and large scale avalanches reduce more and more. Even though, there are still some isolated unstable neurons in the area, it indirectly verifies our deduction of $`p(X)`$ distribution mentioned above.

We investigate the relation between average avalanche size $`S`$ and radius $`d_e`$. From Fig.3 , we can see that with the decreasing of $`d_e`$ , $`S`$ will increases, and with $`d_e`$ approaching 0, the slope of the curve becomes quite large. This result is approximate to the phenomena with the increasing of pulse discharging intensity $`\gamma `$ 7 , and is consistent with the result in Ref. 12 . It implies that the network approaches to SOC state when $`d_e`$ approaching 0 or $`\gamma `$ approaching 0.5.

When the ratio $`T`$ of radius $`d_i`$ with radius $`d_e`$ is changed , the probability distribution $`P(S)`$ and $`P(X)`$ and avalanche distribution are changed too, the tendency is similar to Fig.1, Fig.2. With the increment of $`T`$, the partly-synchronized behavior of system becomes distinctness. The phenomena can also be explained by the expanding of synchronized region , but it leads us to think that the inhibitory lateral connection and excitatory lateral connection have what different effect on synchronized process? It is look like that the inhibitory lateral connection has more important effect, but it is not very clear, there are still a lot of work to do.

### III.2 The approximate behavior in a kind of quasi-OFC earthquake model

OFC earthquake model is a kind of SOC model which has been widely studied in recent years 2 .If we use transformation in our model : $`1\eta _{ij}f_{ij}F_{ij}`$ , then formula (4) become $`F_{i^{}j^{}}F_{i^{}j^{}}+\frac{\gamma }{2}F_{ij}`$ . It is the avalanche mechanism of the OFC earthquake model , in fact these two models belong to the same class.Therefore we add some synchronized regions in OFC model and examine that if the system has the similarity partly-synchronized behavior as in our neuron network model. Hence we introduced a cellular automaton model based on OFC model, the steps are as follows:

1) We define a $`N\times N`$ square lattice ,and a $`L\times L`$ square area in the lattice ( $`L=KN<N`$) .

2) Initialize all sites to a random value $`F_{ij}`$ between 0 and 1.

3) If any $`F_{ij}F_{th}=1`$,then redistribute the force on $`F_{ij}`$ to its neighbors according to the rule:

$`F_{i^{}j^{}}F_{i^{}j^{}}+{\displaystyle \frac{\gamma }{2}}F_{ij}\{\begin{array}{cc}\gamma =\alpha <0.5\hfill & ,whensite(i,j)L\times Lregion,\hfill \\ \gamma =0.5\hfill & ,otherwise.\hfill \end{array}`$ (7)
$`F_{ij}0`$

4) Repeat step 3 until no $`F_{ij}F_{th}`$,we define the avalanche is fully evolved.

5) Locate the site with the largest strain ,$`F_{max}`$ ,then drive all sites

$$F_{ij}F_{ij}+(F_{th}F_{max})$$
(8)
and return to step 3. Here we use open boundary conditions ,and serial working mode.It is same as our neuron network model but different from traditional OFC model.

We still focus on the spatial distribution of $`S`$ and $`X`$. Let $`\alpha =0.1`$, then change $`K`$, the result is shown in Fig 4. It can be seen that the probability distribution $`P(X)`$ increases and deviates from power-law more with $`K`$ increasing at larger $`X`$ in Fig 4.a. The tendency is similar to that in our previous neuron network model. In Fig 4.b , we can see the occurring probability of large scale avalanche reduces with $`K`$ increasing too, but the phenomena is not distinctness as that in previous model. The reason may be that the previous model has more complex structure than this model. It produces many pieces of synchronized region, but here, we only add a piece of synchronized region in model.

A. Corral, Grassberger et al. have investigated the influence of pulse discharging intensity $`\gamma `$ in OFC model. They propose the model present a macroscopic synchronization among all the elements of the lattice when $`\gamma `$ is small 4 ; 13 . Therefore in this model , with expanding of $`L\times L`$ area($`\gamma =\alpha =0.1`$), the synchronized area in the network expands, the synchronization behavior becomes distinctness, and the whole system evolves into a partly-synchronized state. The dynamics behavior of the previous model is so consistent with this model, it suggested that maybe they have the similar dynamics mechanism . Our results tend to agree with this idea.

For studying the partly-synchronized phenomena further, the relation between distribution $`P(x)`$ and $`\alpha `$ has been investigated. Let $`K=0.8`$ , the partly-synchronized behavior become distinctness with $`\alpha `$ decreasing from 0.5 . But when $`\alpha `$ less than a number (about 0.4), there is no more distinctness partly-synchronized behavior introduced by $`\alpha `$ changing (seen in Fig 5) . It can show that the $`L\times L`$ area have evolved into synchronized state. We also draw the whole network’s avalanches distribution map(Fig .6). It can be clearly seen that the isolated avalanches’ number in $`L\times L`$ area increases with $`\alpha `$ decreasing, and the large scale avalanche can’t propagate into the interior of $`L\times L`$ area. These result are consistent with the work of C.Tang , Grassberger et al. 10 ; 13 .

## IV conclusion

In this paper, we analyze the dynamics of the proposed neural network model and find the distribution of the avalanche sizes and the distances $`X`$ between subsequent unstable sites show the power-law behavior. More important, we find the function area and weight distribution produced by self-organized process in our Neural Network model will let the $`P(S)`$ and $`P(X)`$ distribution deviate from power law behavior , and the system evolves into a partly-synchronized state at this time. To verify the explanation , we study a quasi-OFC earthquake model containing synchronized region , and find it will deviate from power-law , evolve into a partly-synchronized state in some conditions too.

Our self-organized feature map Neuron Network model is just a very simple simulation of brain. The real brain has very complex structure and more specific feature regions in Cortex. So brain may express a quasi-SOC(partly-synchronized) behavior more than a pure SOC behavior. Therefore the stored pattern in brain might be designed as a quasi-SOC attractor, and the associate memory process might be designed as the process of the input pattern evolving into the attractor.

Now the neuron synchronization in brain has been observed in many experiments 14 . Rodriguez et al. have investigated the long distance synchronization of human brain activity 15 . We think it would be interesting to further investigate the relationship between synchronization and cognitive , associate process.


text:
# TURKU-FL-P34-00 Q-ball collisions in the MSSM: gravity-mediated supersymmetry breaking

## 1 Introduction

Stable non-topological solitons , Q-balls , can be present in several field theory models. In particular the supersymmetric extensions of the Standard Model may contain them. A Q-ball is a coherent state of a complex scalar field that carries a conserved charge, typically a $`U(1)`$-charge. In the sector of fixed charge a Q-ball is a ground state so that the conservation of charge assures that the Q-ball is stable. In the Minimal Supersymmetric Standard Model (MSSM) Q-balls carrying lepton or baryon number are present due to the existence of flat directions in the scalar sector of the theory .

The cosmological significance of Q-balls can present itself in many forms. Stable (or long living) Q-balls are natural candidates for dark matter  and the decay of Q-balls can explain the baryon to dark matter ratio of the universe . Q-balls can also protect the baryon asymmetry from sphalerons at the electroweak phase transition and decaying Q-balls may be responsible for the baryon asymmetry of the universe . Furthermore, Q-balls can play an important role in considering the stability of neutron stars .

The mechanism by which supersymmetry (SUSY) is broken in the theory is significant for the charges and stability of Q-balls in the theory. If SUSY is broken by a gauge-mediated mechanism, the baryon number carrying B-balls can have very large charges due to to the flatness of the potential. Assuming that the charge is large enough, they can then be stable against decay into nucleons . If, however, SUSY is broken by a gravitationally coupled hidden sector, the potential is not completely flat but Q-balls may still exist due to radiative corrections . In this case the Q-balls can decay (evaporate ) into baryons or supersymmetric particles.

The formation of Q-balls from an Affleck-Dine (AD) condensate in the early universe has been studied recently with numerical simulations . In these simulations both the gauge- and gravity-mediated SUSY breaking scenarios have been considered. In both cases it was found that Q-balls do form from the AD condensate. In the gravity-mediated case it was especially noted that the formed Q-balls have non-zero velocities and can hence collide with each other . Q-ball collisions were also simulated on a one dimensional lattice and it was found that Q-balls typically merge, exchange charge or pass through each other . Since the charge can change due to collisions, they may play an important role in the determination of the Q-ball charge distribution after their formation. On the other hand the charge distribution is important in evaluating the significance of Q-balls for the evolution of the universe. It is hence worthwhile to study Q-ball collisions in more detail. Q-ball collisions have also been studied previously in various potentials in - but not to our knowledge in either the gauge- or gravity-mediated scenarios.

In this paper we have studied numerically collisions of Q-balls in the gravity-mediated scenario on a two dimensional lattice. The gauge mediated scenario will be analyzed in a forthcoming paper.

## 2 Q-ball solutions

Consider a field theory with a U(1) symmetric scalar potential, $`U(\varphi )`$, with a global minimum at $`\varphi =0`$. The complex scalar field $`\varphi `$ carries a unit quantum number with respect to the $`U(1)`$-symmetry. The charge and energy of a given field configuration $`\varphi `$ are given by

$$Q=\frac{1}{i}(\varphi ^{}_t\varphi \varphi _t\varphi ^{})d^Dx$$
(1)
and

$$E=[|\dot{\varphi }|^2+|\varphi |^2+U(\varphi ^{}\varphi )]d^Dx.$$
(2)
The single Q-ball solution corresponds to the minimum energy configuration at a fixed charge. If it is energetically more favourable to store charge in a Q-ball compared to free particles, the Q-ball will be stable against radiative decays into $`\varphi `$-scalars. Hence, for a stable Q-ball, condition

$$E<mQ,$$
(3)
where $`m`$ is the mass of the $`\varphi `$-scalar, must hold.

Minimizing the energy is straightforward and is easily done by using Lagrange multipliers. The Q-ball solution can be shown to be of the form

$$\varphi (x,t)=e^{i\omega t}\varphi (r),$$
(4)
where $`\varphi (x)`$ is now time independent and real, $`\omega `$ is the Q-ball frequency, $`|\omega |[0,m]`$ and $`\varphi `$ is spherically symmetric. The charge of a spherically symmetric Q-ball in D-dimensions reads

$$Q=2\omega \varphi (r)^2d^Dr.$$
(5)
The equation of motion at a fixed $`\omega `$ is

$$\frac{d^2\varphi }{dr^2}+\frac{D1}{r}\frac{d\varphi }{dr}=\varphi \frac{U(\varphi ^2)}{\varphi ^2}\omega ^2\varphi .$$
(6)
To obtain the Q-ball profiles we must solve (6) with boundary conditions $`\varphi ^{}(0)=0,\varphi (\mathrm{})=0`$.

In the present paper we consider a potential of the form

$$U(\varphi )=m^2\varphi ^2(1K\mathrm{log}(\frac{\varphi ^2}{M^2}))+\lambda \varphi ^{10},$$
(7)
where the parameter values are chosen to be $`m=100`$ GeV, $`K=0.1`$ and $`\lambda =M_{\mathrm{Pl}}^6`$, where $`M_{\mathrm{Pl}}`$ is the reduced Planck mass. This choice of potential corresponds to a D-flat direction in the full scalar potential in the MSSM where supersymmetry has been broken by a gravitationally coupled hidden sector . The large mass scale M is chosen such that the minimum is degenerate

$$M=(\frac{1}{4}Km^2\lambda ^1\mathrm{exp}(1\frac{4}{K}))^{\frac{1}{8}}10^{11}\mathrm{GeV}.$$
(8)
We have calculated the charge and energy of Q-balls for different values of $`\omega `$. Energy vs. charge curves are shown in Figure 1(a). The axis scales are chosen differently for two and three dimensions; for two dimensions, $`Q_0=8000(M/\mathrm{GeV})^2`$, $`E_0=8\times 10^5M^2\mathrm{GeV}^1`$ and for three dimensions, $`Q_0=600(M/\mathrm{GeV})^2`$, $`E_0=6\times 10^4M^2\mathrm{GeV}^1`$. The dashed line is the stability line, $`E=mQ`$, that indicates that the Q-balls considered here are stable with respect to scalar decays.

From the figure it can be seen that the energy vs. charge curves are of a similar shape in two and three dimensions. Q-ball profiles are plotted in Figure 1(b) for different values of $`\omega `$ in two and three dimensions. The profiles appear very similar in these two cases. These figures suggest that the collision processes calculated in two spatial dimensions are likely to be similar to collisions in three spatial dimensions. Here it is worth noting that, as from the equation of motion (6) can be seen, a non-zero dissipation term is present for dimensions larger than one. This suggests that collisions in one dimension may differ quite significantly from collision processes in higher dimensions. This also seemed to be the case in the one dimensional simulations we have done.

## 3 Collisions

We have studied collisions of Q-balls with equal charges in the potential (7) . The studied values of $`\omega `$:s were $`\omega /m=0.50,0.60,0.75,0.90,0.99`$. These values correspond to charges $`1450,645,110,12.6,2.24`$ in units of $`(M/\mathrm{GeV})^2`$ in the case with two spatial dimensions so that in terms of $`\varphi `$-scalars the considered range of charges is $`10^{24}10^{26}`$. For the same range of $`\omega `$:s, the charges in the three dimensional space are in the range $`10^{23}10^{26}`$.

The relative phase of the Q-balls is also accounted for. By the relative phase, we mean the difference in individual phases at the point where the distance between Q-balls is at a minimum assuming there is no interaction between them. Position is defined as the point of the maximum value of the amplitude of $`\varphi `$. The relative phase is allowed to have values in the range $`0\mathrm{\Delta }\varphi 2\pi `$. The impact parameter is also varied to study the cross-sections. The Q-balls studied here are of the thick-wall type so that there is no natural definition of the Q-ball size. Therefore we have defined the size of a Q-ball by a Gaussian fit; we fit a Gaussian $`\varphi =Ae^{Br^2}`$ to the profiles and define the radius of the ball as $`R=B^{\frac{1}{2}}`$. The cross-sections quoted here are three dimensional cross-sections with the interaction radius taken from the two-dimensional simulations.

Collisions were simulated on a 2+1 -dimensional lattice. The lattice size typically used was $`200^2`$ with continuous boundary conditions. A 9-point Laplacian operator and a step size of $`5\times 10^3`$ was used in all calculations. Collisions were studied for different initial velocities of Q-balls, $`v=10^3`$ and $`v=10^2`$.

### 3.1 Numerical Results

Collisions can roughly be divided into three types; fusion, charge exchange and elastic scattering. Fusion is defined as a process where most of the initial charge is in a single Q-ball after the collision and the rest of the charge is lost either as radiation or as small Q-balls. By charge exchange we mean a process where Q-balls exchange some of their charge while the total amount of charge carried by the two balls is essentially conserved. An elastic scattering is defined to be a process where less than $`1\%`$ of the total charge is exchanged. After the collision the ratio of the charge in the largest Q-ball to the total initial charge as a function of the relative phase has been plotted in Figures 2 and 3.

From the figures the two different types of processes can be distinguished. In a fusion process typically $`1020\%`$ of the initial charge is lost as radiation and small Q-balls and the rest of the charge is in a single Q-ball. On the other hand, if charge is exchanged the larger Q-ball carries usually less than $`70\%`$ of the total charge. As from Figs. 2 and 3 can be seen, fusion occurs generally only when the relative phase is small and is more likely to occur with smaller $`\omega `$. The amount of exchanged charge decreases substantially with increasing relative phase. Increased velocity does not seem to have a large effect for the range of $`\omega `$ where fusion occurs. However, more charge is exchanged between balls of equal size when the initial velocity is larger. The relative changes in size are weakly dependent on $`\omega `$ (in both cases the standard deviation is less than $`1\%`$) and hence on the size of the Q-balls. Averaging over the relative phase (assuming a random distribution for the $`\mathrm{\Delta }\varphi `$:s) and $`\omega `$:s, the relative change in the size of a Q-ball is $`10\%`$ for $`v=10^3`$ and $`14\%`$ for $`v=10^2`$.

We are now ready to calculate the fusion cross-section, $`\sigma _\mathrm{F}`$, and the cross-section for the charge exchange, $`\sigma _\mathrm{Q}`$, as a function of $`\omega `$. These are plotted in Figure 4 with the geometrical cross-section, $`\sigma _\mathrm{G}`$.

Clearly the fusion cross-section is strongly dependent on $`\omega `$; larger Q-balls fuse more easily than smaller ones. This effect is clearly not explained by the different geometrical sizes of the Q-balls as can be noted from the geometrical cross-section. The fusion cross-section decreases with increasing $`\omega `$ because the Q-balls with higher $`\omega `$ have more regions with differing relative phases. The field dynamics cannot then even out the relative phase differences quickly enough to keep the colliding balls together. From the simulations it can be seen that balls with larger $`\omega `$ are less likely to fuse than balls with the same phase difference but smaller $`\omega `$. The effect of increasing the initial velocity on the cross-sections can also be noted from Fig. 4. The fusion cross-section is slightly decreased as velocity increases while the charge exchange cross-section increases. The increase in $`\sigma _Q`$ is due to the fact that now the Q-balls have more kinetic energy to overcome the repulsion resulting from the relative phase difference.

The total cross-section including all the cases i.e. when the balls fuse, exchange charge or scatter elastically, is also dependent on $`\omega `$, but only weakly. Averaging over the $`\omega `$:s, $`\sigma _{\text{tot}}=0.27\pm 0.01\mathrm{GeV}^2`$ ($`v=10^3`$) and $`\sigma _{\text{tot}}=0.19\pm 0.01\mathrm{GeV}^2`$ ($`v=10^2`$).

We have also studied Q-ball collisions with larger initial velocities for a more limited set of parameter values. When the initial velocity is increased to $`v=10^1`$, the fusion cross-section is reduced significantly from its value when $`v=10^3`$. Furthermore, at such high velocities we also see processes where the Q-balls pass through each other essentially without exchanging any charge. This is a similar process that was reported to occur in one dimension in  and which we have also observed in our one dimensional simulations.

Charge exchange also affects the final velocities of the Q-balls. As charge is exchanged the speed of the Q-balls typically increase. The final velocity of the smaller ball can be quite large; we have often noted final velocities ten times larger than the initial velocity.

## 4 Conclusions

In this paper we have studied Q-ball collisions in the MSSM with supersymmetry broken by a gravitational hidden sector. For the studied range of charges the total cross-section was found to be approximately constant. The cross-section for fusion, $`\sigma _\mathrm{F}`$, appeared to be smaller than the geometrical cross-section, $`\sigma _\mathrm{G}`$, whereas the cross-section for charge exchange, $`\sigma _\mathrm{Q}`$, was larger than $`\sigma _\mathrm{G}`$. In a collision it is hence more probable that a charge exchanging process occurs rather than a fusion process. This probability increases with increasing $`\omega `$ (or with decreasing charge). Averaging over the fusion and charge exchanging processes the average charge increase of the largest Q-ball emerging from a collision was found to be approximately constant. For the considered range of charges and velocities it was $`10\%(v=10^3)`$ and $`14\%(v=10^2)`$.

In a cosmological context, Q-ball collisions may have a significant effect on the charge distribution of Q-balls. Clearly for collisions to be important the number density of Q-balls must be high enough and the balls must have large enough velocities for the rate of interaction to be significant. In the early universe this obviously means that the interaction rate must be larger than the Hubble rate. If collisions typically do occur the resulting charge distribution can then be altered by the fusion and charge exchange processes. Based on the results presented in this paper, the relative phase, size and the initial velocity of the balls then play important roles in studying the evolution of the Q-ball charge distribution.

If the balls that are formed from the AD condensate are in the same phase, fusion processes will dominate and the average size of a Q-ball grows substantially in a collision. Since most of the charge is left in the remaining ball and the rest is in the form of several small, quickly evaporating , Q-balls and radiation, the number density of Q-balls reduces rapidly. Collisions can therefore freeze the Q-ball distribution quickly in the early phases of the universe. If, on the other hand, the phases are randomly distributed the probability for fusion is greatly reduced and the distribution will not change as significantly as in the previous case. Collision processes can then also continue for a longer period of time.

The typical size of Q-balls is obviously an important factor. The total scattering cross-section depends quite weakly on the size of the Q-balls but the fusion and charge exchange cross-sections do have a strong $`\omega `$-dependence.

The initial velocity of the Q-balls is also significant in the evolution of the Q-ball distribution. A larger velocity means that the interaction rate is increased but on the other hand if the initial velocity is too large the cross-sections decrease due to a decreased interaction time. Collisions can also significantly change the velocities of the Q-balls so that an initially uniform velocity distribution can be spread out by the collision processes.

The effect of collisions can be important in deciding the exact role and significance of Q-balls in the evolution of the universe. In determining their importance on cosmology, more information is needed about the Q-ball distribution after their formation and also about the effects of collisions on the initial distribution. To quantify the effects of the different collision processes described in this paper, a more detailed analysis is needed which gives motivation for future work.

Acknowledgements. We thank K. Enqvist for discussions and the Center for Scientific Computing for computation resources. This work has been supported by the Academy of Finland.


text:
# Enhanced Stability of Layered Phases in Parallel Hard-Spherocylinders Due to Addition of Hard-Spheres

## I Introduction

In hard particle fluids all allowed configurations have the same energy and therefore it is the number of states, or equivalently the entropy of a system that determines the equilibrium phase. Examples of well known phase transitions where the formation of ordered structures are driven solely by an increase in entropy are the liquid to crystal transition in hard spheres , the isotropic to nematic  and the nematic to smectic transition in hard-rods . Because of their high degree of monodispersity, and because of the dominant role of steric repulsion in the pair-potential, colloidal suspensions of polystyrene latex and rod-like viruses have often been used as experimental model systems for the study of entropy induced ordering in hard-sphere  and hard-rod systems, respectively .

A natural extension of the above work is to the phase behavior of mixtures, with a number of recent experimental and theoretical studies focusing on the phase behavior of binary mixtures of hard-spheres . We have recently begun work on less studied systems that closely approximate hard-rod/hard-sphere and hard-rod/polymer mixtures . As a model for hard-rods we used either $`\mathrm{𝑓𝑑}`$ or TMV virus, as hard-spheres we used polystyrene latex, and as polymers we used poly(ethylene-oxide) with varying molecular weights . The part of the phase diagram explored consisted of pure rods in either the isotropic, nematic, or smectic phase to which a small volume fraction of spheres or polymers was added. Remarkably, besides the expected uniform mixtures and bulk demixing, we also observed a variety of microphases for a wide range of sphere sizes and concentrations . In microphase separation the system starts separating into liquid-like regions that are rich in either spheres or rods. However, unlike bulk demixing where rod and sphere rich regions grow until reaching macroscopic dimensions, in microphase demixing these liquid-like regions increase only to a critical size after which they order into well defined three dimensional equilibrium structures. One of the micro-separated phases observed, named the lamellar microphase, consists of alternating two-dimensional liquid-like layers of rods and spheres and is the subject of theoretical analysis in this paper.

In this paper we use the second virial approximation first studied by Koda et. al.  to examine the influence of molecular parameters such as shape and size, on the phase behavior of rod/sphere mixtures. As the second-virial theory is approximate in nature, we validate the theoretical predictions by comparing them with either computer simulations or experimental results. The remainder of this paper is organized as follows: In section II we formulate the second virial approximation for the rod/sphere mixture. The general features of the phase diagram are discussed and a physical picture of the factors responsible for the enhanced stability of the layered phase due to the presence of spheres is presented. In Section III the influence of varying the spherocylinder length on the phase behavior of spherocylinder/sphere mixtures is studied using computer simulations and the results are compared to theoretical predictions. Section IV examines how changes of the sphere diameter influence the phase behavior of spherocylinder/sphere mixtures. Finally in Section V we present our conclusions.

## II General features of a phase diagram of a spherocylinder-sphere mixture

Although the equilibrium phases of all hard particle fluids are determined by maximizing the entropy, ordering transitions are still possible because the expression for the total entropy, or equivalently free energy, splits into two parts. The ideal contribution to the entropy is of the form $`\rho \mathrm{ln}\rho `$, where $`\rho `$ is the density distribution function. This contribution to the entropy attains a maximum for a uniform density distribution and therefore always suppresses transitions from uniform to modulated phases. In contrast, excluded volume entropy sometimes increases with increasing order and therefore drives the system towards a modulated phase. In this paper we use a highly simplified second virial approximation to calculate the excluded volume entropy.

The equilibrium phase in a spherocylinder/sphere mixture is determined by four parameters: length over diameter of a spherocylinder $`(L/D_{sc})`$, diameter of spherocylinder over diameter of sphere ($`D_{sc}/D_{sp}`$), total volume fraction of spheres and spherocylinders ($`\eta `$) and partial volume fraction of spheres ($`\rho _{sp}`$). To help us in interpretation of our results we first define the slope

$$\tau =\underset{\rho _{sp}0}{lim}\frac{\eta (\rho _{sp})\eta (0)}{\rho _{sp}}$$
(1)
where $`\eta (\rho _{sp})`$ is the total volume fraction of the rod-sphere mixture at the layering transition after the introduction of spheres at partial volume fraction $`\rho _{sp}`$. A positive value of $`\tau `$ implies that adding a second component stabilizes the nematic phase by displacing the smectic transition to higher densities. For the case when both components are spherocylinders of different lengths but with the same diameter, slope $`\tau `$ is positive if the ratio of lengths is less then approximately 7 . In the same manner, negative values of $`\tau `$ imply that the second component stabilizes the smectic phase. There are predictions of a negative value of $`\tau `$ in a bidisperse rod mixture when the ratio of rod lengths is large enough , or when added rods have a larger diameter . In this section we focus on the phase behavior of the spherocyinder-sphere mixture for the specific microscopic parameters $`L/D_{sc}=20`$ and $`D_{sc}/D_{sp}=1`$. We present a physical picture of excluded volume effects that are responsible for the enhanced stability of the lamellar phase. In the next two sections we extend our study on how changes in the molecular parameters $`L/D_{sc}`$ and $`D_{sc}/D_{sp}`$ modify the phase behavior and in particular, their influence on the magnitude and sign of the slope $`\tau `$.

### A Second virial approximation

The second virial approximation for a mixture of perfectly aligned spherocylinders and spheres of equal diameter was proposed by Koda, Numajiri and Ikeda  and is generalized for arbitrary $`L/D_{sc}`$ and $`D_{sc}/D_{sp}`$ in the appendix. It was previously shown that the second virial approximation described qualitatively the formation and various features of the smectic phase of hard rods . Here we study how the addition of spheres perturbs the formation of the smectic phase. Since the sphere volume fraction is very low we expect that the second virial approximation is still qualitatively correct for these mixtures. We consider a sinusoidal perturbation from the uniform density for both spherocylinders and spheres. From equations (7) and (16) in the appendix we obtain the free energy difference between the uniformly mixed and layered state in a spherocylinder/sphere mixture:

$$\delta F=a_1^2\left(S_{11}2\frac{a_1}{a_2}S_{12}+\left(\frac{a_1}{a_2}\right)^2S_{22}\right)=0$$
(2)
.

The phase diagram obtained within this approximation for microscopic parameters $`L/D_{sc}=20`$ and $`D_{sc}/D_{sp}=1`$ is shown in Fig. 1. From the phase diagram we see that the first prediction of the model is that spheres, upon addition to a smectic phase, will preferentially occupy space between smectic layers and therefore create a stable micro-separated lamellar phase. The second prediction is that the total volume fraction at which the system undergoes a transition from a uniform miscible state to a layered lamellar state is lowered by increasing the partial volume fraction of spheres. This implies that the slope $`\tau `$ is negative for this particular spherocylinder/sphere mixture and we conclude that in this case spheres enhance the layering transition.

We can assign a simple physical origin to every term given in Eq. (2) above and Eq.(12) of the Appendix. The parts of the spherocylinder-spherocylinder interaction term $`S_{22}`$ and sphere-sphere term $`S_{11}`$ that scale as $`\eta `$ are due to the ideal (id) contribution to the free energy, also known as the entropy of mixing and are denoted as $`S_{22}^{id}`$ and $`S_{11}^{id}`$, respectively. The terms having a $`\eta ^2`$ dependence in $`S_{22}`$, $`S_{12}`$, $`S_{11}`$ are due to the spherocylinder-spherocylinder, spherocylinder-sphere and sphere-sphere excluded volume (ex) interaction, respectively and are denoted as $`S_{22}^{ex}`$, $`S_{12}^{ex}`$ and $`S_{11}^{ex}`$. Since the instability is defined as $`\delta F(\eta _c,k_c)=0`$, at a critical density $`\eta _c`$ and at a critical wavevector $`k_c`$ all individual contributions to the free energy difference in Eq. (2) must add up to zero. In Fig. 2 we show the value of all terms with distinct physical origins at the instability density $`\eta _c`$ and wavevector $`k_c`$ as a function of partial volume fraction of spheres. Since from our analysis we cannot determine the absolute amplitude of $`a_1`$ we only plot the ratios of all free energy components to the absolute value of the spherocylinder-spherocylinder excluded volume $`S_{22}^{ex}`$. If we set the partial volume fraction of spheres to zero $`(\rho _{sp}=0)`$ in Eq.(2) we obtain an equation whose solution indicates the nematic-smectic stability limit in a pure suspension of aligned spherocylinders . For these conditions the only two nonzero components of free energy are $`S_{22}^{ex}`$, which is negative and therefore drives the transition and $`S_{22}^{id}`$, which is positive and therefore suppresses the transition. As we start increasing the partial sphere volume fraction $`\rho _{sp}`$, the spherocylinder-sphere free volume term $`S_{12}^{ex}`$ rapidly assumes large negative values as evidenced by the rapidly decreasing ratio of $`S_{12}^{ex}/S_{22}^{id}`$. This implies that layering the mixture significantly decreases the excluded volume that is due to the spherocylinder-sphere interaction.

We can use the information gained from the second virial approximation to obtain a clear physical picture of excluded volume effects in spherocylinder/sphere mixtures and explain the enhanced stability of the lamellar phase. Taking any single spherocylinder in a uniform spherocylinder/sphere mixture and replacing it by two spheres will leave the value of excluded volume virtually unchanged. The reason for this lies in the fact that the volume excluded to the spherocylinder due to the presence of a sphere with equal diameter, under the constraint of uniform packing, is a spherocylinder with diameter $`2D_{sc}`$ and length $`(L+2D_{sc})`$ where $`L`$ and $`D_{sc}`$ are defined in Fig. 3. However, the excluded volume between any two spherocylinders with large $`L/D_{sc}`$ is only about twice this value as illustrated in Fig. 3. Although replacing spherocylinders by spheres in such a manner leaves the excluded volume almost unchanged, it does significantly decrease the total volume fraction of the mixture since the volume of two spheres is much smaller then the volume of a spherocylinder with large $`L/D_{sc}`$. Therefore in the spherocylinder/sphere mixture we encounter excluded volume problems similar to those found in a pure spherocylinder solution, but at a lower total volume fraction. As in pure spherocylinders, the system reduces the excluded volume by undergoing a transition to a layered phase. The excluded volume is reduced in the lamellar state because a periodic density distribution forces spheres and spherocylinder into alternate layers thus decreasing the probability of the very unfavorable sphere-spherocylinder contacts as illustrated in Fig. 4. This explains the large decrease in the value of the $`S_{12}^{ex}`$ term at the lamellar transition that we observe in the second virial theory. This term is responsible for the enhanced stability of the lamellar phase in a sphere/spherocylinder mixture. In conclusion, it is the inability to efficiently pack a uniform mixture of spherocylinders and spheres, as reflected in the large spherocylinder/sphere excluded volume term, that destabilizes the nematic phase and enhances the formation of a layered phase.

An alternate way to think about the formation of a layered phase is to focus on the effects of spherocylinder ends . The nematic phase in our simplified model is characterized by random distribution of spherocylinders along their axial and radial direction as illustrated in Fig. 5. This end effect is responsible for the formation of the smectic phase, which is characterized by a periodic density distribution. In similar fashion, introducing a sphere into the nematic phase will have the same effect on the surrounding spherocylinders as another spherocylinder end. Therefore adding spheres very effectively increases the density of “spherocylinder ends” and decreases the total volume fraction. To resolve the difficulties in efficient packing due to these extra “spherocylinder ends”, the mixture layers at a lower total volume fraction.

### B Monte Carlo Simulation

In the previous section we discussed two predictions of the second virial theory for a spherocylinder/sphere mixture with $`L/D_{sc}=20`$ and $`D_{sp}/D_{sc}=1`$; the existence of the lamellar phase and the enhanced stability of the lamellar phase when compared to a smectic phase of pure spherocylinders. Our results are in agreement with previous studies of Koda et. al. . However, the second virial approximation is highly approximate and there is reasonable concern about the influence of higher terms on the topology of the phase diagram. To support their conclusions Koda et. al. performed computer simulations, which indicated the existence of an lamellar phase . Still, the question of whether spheres simply fill the voids between layers in an already formed smectic phase, or actually induce layering at lower total volume fraction was not addressed. In this section, using Monte Carlo simulations we address the question of the influence of adding spheres on the phase behavior of spherocylinders by determining the slope $`\tau `$ in Eq. 1 in a mixture of spherocylinders and spheres with parameters $`L/D_{sc}=20`$ and $`D_{sp}/D_{sc}=1`$

A Monte Carlo simulation of a mixture of hard-spheres and perfectly aligned hard-spherocylinders was performed at constant pressure and number of particles . Most simulations contained 392 spherocylinders and a variable number of spheres. To check for finite size effects we also ran simulations with 784 spherocylinders, but saw no significant difference in the results obtained. In one sweep, pressure was increased from a dilute homogeneous mixture up to a well ordered, dense smectic or lamellar phase. At each value of the pressure, the density of spheres and spherocylinders and their corresponding smectic order parameter were measured after the system was allowed to equilibrate. Identical results were obtained when the pressure was slowly decreased from a initially dense phase composed of alternating layers of spherocylinders and spheres to a dilute homogeneous mixture.

Besides lamellar transitions there is a possible demixing transition where spherocylinders and spheres phase separate into macroscopically distinct phases. However, once a layered phase is formed the exchange of spheres between layers drops to a negligible amount, leaving open the possibility that system would undergo a demixing transition, but is stuck in a lamellar phase, which is only a metastable state. To find out the location of the demixing transition it is necessary to measure the chemical potential of both spherocylinders and spheres in a spherocylinder/sphere mixture . This possibility was not examined in this work, primarily because we are only interested in how low concentrations of spheres perturb the formation of the layered phase. Therefore it is reasonable to expect that at a very low volume fraction of spheres, the lamellar transition is going to be more stable than the demixing transitions as predicted by the second virial theory.

A plot of the smectic order parameter for spherocylinders with $`L/D_{sc}=20`$ as a function of increasing total density for different partial volume fractions of spheres is shown in Fig. 6. As the system approaches a certain critical density we observe a rapid non-linear increase in the smectic order parameter that we interpret as a signature of the nematic to smectic phase transition. This critical density shifts to lower values of the total volume fraction as the partial volume fraction of spheres is increased. To reconstruct a phase diagram from the above data we define a phase as layered when its smectic order parameter reaches a value of 0.3 . For a pure spherocylinder suspension this value yields good agreement with previous studies of the volume fraction of the nematic-smectic phase transition . Since we are mostly interested in the qualitative behavior of a spherocylinder/sphere mixture this method should suffice our purposes. Using this phenomenological rule, the phase diagram for a mixture of spherocylinders and spheres ($`L/D_{sc}=20`$, $`D_{sc}/D_{sp}=1`$) is reconstructed and compared to the second virial theory in Fig. 1. An immediate conclusion drawn from Fig. 1 is that adding spheres to aligned spherocylinders enhances the stability of the lamellar phase, which is indicated by the negative value of slope $`\tau `$, in agreement with the prediction of the second virial approximation.

## III The effects of spherocylinder length on the phase diagram

Next we proceed to investigate the influence of varying the spherocylinder length on the magnitude of slope $`\tau `$. The predictions of the second virial theory for the nematic-lamellar instability are shown in Fig. 7a. The second virial theory clearly predicts increasing stability of the lamellar phase with increasing length of spherocylinder. To verify this prediction we repeated Monte Carlo simulations for spherocylinders with different $`L/D_{sc}`$ and used the same rule as before to identify the volume fraction of the nematic-lamellar transition. The simulation results for the location of the nematic to layered transition are shown in Fig. 7b. We can conclude that our simulations confirm predictions of the second virial model and that the length of the spherocylinder is an important parameter in forming the lamellar phase, with longer spherocylinders showing an increasing tendency to form a layered phase at a lower volume fraction of added spheres.

Using the physical picture of the excluded volume effects developed in the previous section provides a natural explanation for our simulation results in Fig. 7. With increasing spherocylinder length the excluded volume due to the spherocylinder-sphere interaction grows proportionally to the spherocylinder length and consequently the value of the $`S_{12}^{ex}`$ term increases in magnitude. As we have seen before, the larger the $`S_{12}^{ex}`$ term, the more likely it is for the system to form a layered phase.

It is interesting to consider the limit of spherocylinders with infinite aspect ratio. In the density regime of the nematic-smectic transition, this model can be mapped onto a system with skewed cylinders with an aspect ratio close to one. The nematic-smectic transition in this model has been studied numerically . If we consider the addition of spheres to this system, then the same affine transformation that maps the infinite spherocylinders onto squat, skewed spherocylinders, will map the spheres onto infinitely thin, parallel disks. As the disks are infinitely thin, they do not interact with each other but only with the cylinders. Inside the nematic phase, most volume is excluded for these disks. However, in the smectic phase, there is ample space for the disks between the layers. In fact, the stronger the layering, the larger the accessible volume. Hence in this limit, the addition of spheres will strongly stabilize the smectic phase.

## IV The effects of sphere diameter on the phase diagram

In this section we investigate the influence of sphere diameter on the value of slope $`\tau `$. Fig. 8 shows the prediction of the second virial theory for the dependence of slope $`\tau `$ on the ratio of spherocylinder to sphere diameter $`(D_{sc}/D_{sp})`$ for spherocylinders with different $`L/D_{sc}`$. In section A we examine the phase behavior of sphere/spherocylinder mixtures when the sphere diameter is smaller then spherocylinder diameter and in section B we examine the other case when the sphere diameter is larger then the spherocylinder diameter. In our model the presence of the spheres cannot alter the orientational distribution function of spherocylinders, which are always perfectly parallel to each other. It is reasonable to expect that this assumption holds for spheres smaller then the spherocylinder length, but as a sphere becomes larger then the spherocylinder length, long wavelength elastic effects start to dominate the behavior of the system and hard spherocylinders will tend to align parallel to the surface of the sphere . Therefore in Fig. 8 we plot the values of slope $`\tau `$ only for those values of $`D_{sc}/D_{sp}`$ for which our assumptions are at least qualitatively correct. As we increase the sphere size beyond this limit our model describes a highly artificial system of large spheres and parallel spherocylinders. In this regime we observe oscillations in the value of slope $`\tau `$ similar to what is observed in binary mixtures of parallel spherocylinders .

### A Sphere diameter smaller than spherocylinder diameter

In the regime where $`D_{sc}/D_{sp}>1`$ (for spherocylinders of any $`L/D_{sc}`$), decreasing the sphere size increases the stability of the lamellar phase as indicated by the increasing negative value of slope $`\tau `$ seen in the right hand side of Fig. 8. This prediction of the theory has a simple explanation in our picture of excluded volume in a sphere/spherocylinder mixture. If we halve the sphere radius $`D_{sp}`$, while keeping constant the volume fraction of spheres, we increase the number of spheres eight times. At the same time, the result of reducing the sphere size is to decrease the excluded volume of the spherocylinder-sphere interaction. However, the eightfold increase in the number of spherocylinder-sphere interactions more then compensates for the decrease in excluded volume between the sphere and spherocylinder and consequently the magnitude of $`S_{12}^{ex}`$ increases with decreasing sphere diameter. This leads to the increased stability of the layered phase with decreasing sphere size.

It becomes difficult to verify this prediction using computer simulations. As the sphere size decreases at constant total volume fraction $`\eta `$, the number of particles in a simulation rapidly reaches the order of thousands requiring simulation times that are prohibitively long. As the ratio of spherocylinder to sphere diameter $`(D_{sc}/D_{sp})`$ was varied within the accessible range between 0.5 to 2 we did not observe any changes in the value of slope $`\tau `$ that were larger than our measurement error. Larger and longer simulations are needed for a careful analysis of spherocylinder/sphere mixtures with extreme values of the ratio $`D_{sc}/D_{sp}`$.

### B Sphere diameter larger than spherocylinder diameter

For spherocylinders with small $`L/D_{sc}`$, Fig. 8 shows that the magnitude of slope $`\tau `$ uniformly decreases with increasing sphere size. Eventually the slope $`\tau `$ changes sign and becomes positive, implying that large spheres stabilize the nematic and not the smectic phase. The phase diagram under conditions where slope $`\tau `$ is positive is shown in Fig. 9. The wavevector associated with the layering transition, indicated with a solid line in Fig. 9, remains at an almost constant value. Another important point is that the amplitude ratio in Eq. (16) is positive. This means that the periodic density modulations of the spherocylinders and spheres are in phase, which implies that spheres no longer go into the gap between two spherocylinder layers, but rather fit into the spherocylinder layer. However, as the partial volume fraction of spheres $`(\rho _{sp})`$ is increased further we observe a discontinuous jump in the wavevector to zero value. This implies that there is a discontinuous change from a layering to a demixing transition. As the demixing transition is reached there is also a change in sign of the amplitude ratio, which becomes negative and the spherocylinders and spheres bulk separate. In contrast, the phase diagram for mixtures of small spheres and spherocylinders shown in Fig. 1 looks quite different. The amplitude ratio for this case is always negative implying formation of the lamellar phase. Another contrast is that in a mixture of small spheres and spherocylinders the wavevector associated with the layering transition decreases in a continuous fashion until it reaches zero value.

We now examine the behavior of individual terms in Eq. (2) for a mixture of large spheres and short spherocylinders shown in Fig. 9. Most notably, we find that at low volume fractions of spheres where the system undergoes the layering transition, the ratio $`S_{12}^{ex}/S_{22}^{ex}<<1`$. This implies that upon layering there is almost no reduction of the unfavorable sphere/sperocylinder interaction and that the spherocylinder/spherocylinder interaction alone drives the formation of the layered phase. In contrast, for small spheres this ratio was large and was responsible for enhanced stability of the lamellar phase as was shown in Fig. 2. At a higher volume fraction of large spheres where the mixture directly bulk phase separates we find that the ratio $`S_{12}^{ex}/S_{22}^{ex}>>1`$. This implies, as expected, that demixing very effectively reduces the unfavorable sphere/spherocylinder interactions. These results suggest a physical picture of the excluded volume effect. Unlike small spheres, large spheres can not fit into the gap between smectic layers and consequently there is no way to gain free volume by undergoing the layering transition. As an alternative, to gain free volume the system bulk phase separates at the lowest volume fraction of spheres possible.

While for short spherocylinders the magnitude of slope $`\tau `$ uniformly decreases with increasing sphere size, longer spherocylinders exhibit a qualitatively different behavior. For a mixture of spherocylinders with $`L/D_{sc}=100`$ and spheres with $`D_{sc}/D_{sp}=0.1`$ there is a pronounced increase in the stability of the lamellar phase as shown in Fig. 8. By increasing the length of spherocylinders to even larger values, the region of increased stability of the lamellar phase shifts to higher values of the sphere radius. Two conditions emerge, which when satisfied lead to enhanced stability of the lamellar phase. First, it is necessary for a sphere to fit between two smectic layers without disturbing them. This condition is satisfied when $`D_{sp}/L0.1`$. The second condition is that $`D_{sp}/D_{sc}>>1`$. It was argued before that under these condition large spheres are able to induce smectic correlations amongst neighboring spherocylinders , which in turn can enhance the formation of the lamellar phase.

Because of the large size asymmetry it was not feasible to carry out simulations for mixture of spherocylinders and spheres with $`L/D_{sc}100`$ and $`D_{sc}/D_{sp}0.1`$. However, these conditions are closely approximated by recent experiments on rod-like $`\mathrm{𝑓𝑑}`$ ($`L=1\mu `$m, $`L/D_{sc}`$ 100) and polystyrene spheres . Therefore, we compare theoretical results of slope $`\tau `$ for spherocylinders with $`L/D_{sc}=100`$ shown in Fig. 8 to these experimental results . When large spheres $`D_{sp}1\mu `$m, $`(D_{sc}/D_{sp}0.01)`$ are mixed with $`\mathrm{𝑓𝑑}`$ at any concentration for which the nematic phase is stable, we observe no formation of the layered phase. Instead, large spheres phase separate into dense aggregates elongated along the nematic director indicating that the value of slope $`\tau `$ is larger then zero. When the size of the sphere was decreased to $`D_{sp}=0.1\mu `$m, $`(D_{sc}/D_{sp}10)`$ we observed a transition to a layered state at a $`\mathrm{𝑓𝑑}`$ concentration of 20 mg/ml. The formation of a smectic phase in a pure fd suspension at the same ionic strength occurs at 65 mg/ml. The fact that adding spheres diminishes the rod density by a factor of three indicates a large negative value of slope $`\tau `$. As the sphere size was further decreased $`D_{sp}=0.022\mu `$m, ($`D_{sc}/D_{sp}=0.46`$) there was again indication of a lamellar phase, but this time at a much higher concentration of rods of about 50 mg/ml. Thus, although small spheres still stabilize the layering transition, implying a negative value of slope $`\tau `$, the magnitude of slope $`\tau `$ is much less for $`D_{sc}/D_{sp}0.46`$ then for $`D_{sc}/D_{sp}0.1`$. These qualitative trends of the non-monotonic behavior of slope $`\tau `$ with sphere size observed in experiments of fd-polystyrene mixtures are very similar to the theoretical prediction shown in Fig. 8 for spherocylinders with $`L/D_{sc}=100`$.

## V Conclusions

In this paper we have presented the predictions of the second virial theory for a mixture of parallel hard-spherocylinders and hard-spheres undergoing one dimensional microphase separation. We have been able to verify a number of these predictions using Monte Carlo simulations. We found that spheres induce layering, which implies a negative value of the slope $`\tau `$, which is the change in total volume fraction of the mixture at the point of nematic-smectic instability with respect to the partial volume fraction of added spheres (Eq. 1) . At the same time the magnitude of the slope $`\tau `$ increases with increasing spherocylinder length. In other words, spheres at the same partial volume fraction stabilize layering of longer spherocylinders more then shorter spherocylinders. Besides this, the theory predicts an unusual non-monotonic behavior in slope $`\tau `$ as a function of sphere to spherocylinder diameter. Although the physical origin of this effect is not clear, it is intriguing that similar qualitative trends are observed in experiments of mixtures of the spherocylinder-like $`\mathrm{𝑓𝑑}`$ and polystyrene spheres. However, in real experiments spherocylinders are free to rotate, are flexible, and have charge associated with them. Before quantitative comparisons with experiments are possible it will be necessary to perform simulations and formulate theories that take into account these effects mostly ignored in this highly idealized treatment.

## VI Acknowledgments

We acknowledge useful discussions with Richard Sear and Bulbul Chakraborty. We thank Tomonori Koda for critical reading of the manuscript. This research was supported by NSF DMR-9705336 and NSF INT-9113312. The work of the FOM Institute is supported by FOM (“Stichting Fundamenteel Onderzoek der Materie”) with financial aid from NWO (“Nederlandse Organisatie voor Wetenschappelijk Onderzoek”).

## VII Appendix

A general expression for the free energy of bidisperse mixture at the second virial level is

$`\beta F(\rho _1,\rho _2)={\displaystyle \underset{i=1,2}{}}{\displaystyle _V}d(𝐫)\rho _i(𝐫)\mathrm{ln}(\rho _i(𝐫))`$ (3)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1,2}{}}{\displaystyle \underset{j=1,2}{}}{\displaystyle _V}𝑑𝐫_\mathrm{𝟏}{\displaystyle _V}𝑑𝐫_\mathrm{𝟐}\rho _i(𝐫_\mathrm{𝟏})\rho _j(𝐫_\mathrm{𝟐})f_{i,j}(𝐫_\mathrm{𝟏},𝐫_\mathrm{𝟐})`$ (4)
where the function $`f_{i,j}`$ is the overlap function between two spheres, sphere and spherocylinder or two spherocylinders . It attains the value of -1 if two particles overlap, otherwise it is equal to 0. The terms involving $`\rho \mathrm{ln}\rho `$ represent the entropy of mixing while the terms involving $`f_{i,j}`$ represent the free volume entropy. Since we are interested in one dimensional layering we look at the response of the system to following density pertubation

$`\delta \rho _1(z)=a_1\mathrm{cos}(k_zz)`$ (5)
$`\delta \rho _2(z)=a_2\mathrm{cos}(k_zz)`$ (6)
The free energy difference between the uniform and perturbed state is

$$\delta F=F(1+\delta \rho _1(z),1+\delta \rho _2(z))F(1,1)=\stackrel{~}{𝐚}\mathrm{𝐒𝐚}$$
(7)
where $`\stackrel{~}{𝐚}=(a_1,a_2)`$ and $`𝐒`$ is a two dimensional stability matrix. To find the limit of stability we have to solve the equation $`\text{det}(𝐒)=0`$. For latter convenience we define the following function

$`S({\displaystyle \frac{L}{D_{sc}}},\sigma ,k)`$ $`=`$ $`{\displaystyle \frac{3\mathrm{sin}(k\sigma (2+2\frac{L}{D_{sc}}))}{4k^3}}`$ (9)
$`{\displaystyle \frac{2k\sigma \mathrm{cos}(k\sigma (2+2\frac{L}{D_{sc}}))\mathrm{sin}(k2\sigma \frac{L}{D_{sc}})}{4k^3}}`$
.

The above expression depends only on geometrical factors and is related to the Fourier transform of the spherocylinder which is specified by the excluded volume between a sphere of diameter $`D_{sp}`$ and a spherocylinder of length $`L`$ and diameter $`D_{sc}`$. Wavevector $`k`$ is dimensionless because it is rescaled with the spherocylinder diameter ($`D_{sc}`$). The parameter $`\sigma `$ is defined as ratio of sphere diameter to spherocylinder diameter $`(\sigma =D_{sp}/D_{sc})`$. In the limit of $`L/D_{sc}0`$ the above expression reduces to a Fourier transform of a sphere with unit diameter. The stability matrix $`𝐒`$ for a mixture of spherocylinders and spheres has the following form

$`𝐒=(\begin{array}{cc}\frac{\eta \left(1\rho _{sp}\right)\left(1+4\left(1\rho _{sp}\right)\eta S(0,1,k)\right)}{4}& \\ \frac{2\rho _{sp}\left(1\rho _{sp}\right)\eta ^2S({\displaystyle \frac{L}{D_{sc}}},1+\sigma ,k)}{\sigma ^6\left({\displaystyle \frac{2}{3}}{\displaystyle \frac{L}{D_{sc}}}+1\right)^2}& \end{array}`$ (12)
$`\begin{array}{cc}\frac{2\rho _{sp}\left(1\rho _{sp}\right)\eta ^2S({\displaystyle \frac{L}{D_{sc}}},1+\sigma ,k)}{\sigma ^6\left({\displaystyle \frac{2}{3}}{\displaystyle \frac{L}{D_{sc}}}+1\right)^2}& \\ \frac{\eta \rho _{sp}\left({\displaystyle \frac{\sigma ^6({\displaystyle \frac{3}{2}}{\displaystyle \frac{L}{D_{sc}}}+1)}{4}}+\eta \rho _{sp}S(2{\displaystyle \frac{L}{D_{sc}}},2\sigma ,k)\right)}{\sigma ^6\left({\displaystyle \frac{2}{3}}{\displaystyle \frac{L}{D_{sc}}}+1\right)^2}& \end{array})`$ (15)
where $`\rho _{sp}`$ denotes partial volume fraction of spheres and varies between 0 and 1 while $`\eta `$ denotes total volume fraction. Note that the terms in matrix elements $`S_{11}`$ and $`S_{22}`$ proportional to $`\eta `$ are due to configurational entropy while terms proportional to $`\eta ^2`$ are due to free volume entropy. As $`k0`$ the condition $`\text{det}(𝐒)=0`$ reduces to the usual thermodynamic condition for the stability of the system against bulk phase separation.

To reconstruct the stability diagram from the determinant we slowly increase the total volume fraction $`\eta `$. At a certain value of total volume fraction $`(\eta _c)`$ the determinant of $`𝐒`$ will equal zero for a specific wavevector $`(k_c)`$. If the wavevector $`k_c`$ obtained has a finite value it implies that system is undergoing a layering transition. On the other hand, the condition $`\text{det}(𝐒)=0`$ when $`k_c=0`$ implies complete demixing. Once we obtain values of $`\eta _c`$ and $`k_c`$ we can find out the ratio of amplitudes from the following formula

$$\frac{a_1}{a_2}=\frac{S_{12}(\eta _c,k_c)}{S_{11}(\eta _c,k_c)}.$$
(16)
A positive value of the amplitude ratio implies that the spheres and spherocylinders are in the same layer (the periodic modulations are in phase), while a negative value implies that the spheres and spherocylinders intercalate (the periodic modulations are out of phase).


text:
# CORRELATED MIXTURES OF ADIABATIC AND ISOCURVATURE COSMOLOGICAL PERTURBATIONS

## 1 Introduction

The CMB anisotropies can indirectly measure the cosmological parameters by looking at the evolution of the cosmological perturbations between the end of inflation and the recombination epoch. In order to do so, one must implicitly assume a simple form for the initial perturbations. Usually, one considers only adiabatic fluctuations, with a power law spectrum.

Adiabatic fluctuations arise in the context of the simplest inflationary scenario. However, as soon as one has a multiple inflation scenario, isocurvature fluctuations can be generated, the amplitude of which, as well as their correlation with the adiabatic part, depend on the parameters of the model. In many of the models already studied $`^\mathrm{?}`$, the isocurvature and the adiabatic parts of the fluctuations are uncorrelated, but it is possible to have a correlated mixture of such perturbations, as was already stressed by one of us $`^\mathrm{?}`$ in the study of a specific inflation model with two massive non interacting scalar fields.

In this communication, we study this issue in a phenomenological way, and look at the observable consequences of a correlated mixture of adiabatic and isocurvature cosmological perturbations. These perturbations happen to have a richer structure than the more usual uncorrelated adiabatic and isocurvature ones.

## 2 Definitions and notations

When considering a mixture of several fluids, one can define an entropy perturbation $`S_{A,B}`$ for any pair of components $`A`$ and $`B`$. It is non zero as soon as the different particle number density contrasts $`\delta n_X/n_X`$ ($`X=A,B`$), are not equal. It can also be written in terms of energy density contrasts $`\delta _X`$ :

$$S_{A,B}\frac{\delta n_A}{n_A}\frac{\delta n_B}{n_B}=\frac{\delta _A}{1+\omega _A}\frac{\delta _B}{1+\omega _B},$$
(1)
where $`\omega _Xp_X/\rho _X`$ is the equation of state parameter for the species $`X`$. (Note that this definition is gauge-invariant.) Adiabatic initial conditions are defined such that all the entropy perturbations are zero. In cosmology, the perturbations are actually considered as random fields, usually assumed to be Gaussian, and described by their power spectra. When one has to deal with several random fields, one must also impose the form of the cross-correlation between them. In what follows, we will consider only *totally correlated* adiabatic and isocurvature perturbations, where all the random fields are described in terms of a single random variable. We further assume that only one species deviates from adiabaticity. In this case, one simply has to define which species deviates from adiabaticity, say $`X`$ (i.e. $`S_{A,B}=0`$ when $`A,BX`$), and the relative initial amplitude between the entropy perturbation $`S_{X,Y}`$ and the Bardeen potential $`\mathrm{\Phi }`$ (where $`Y`$ is another species which does not deviate from adiabaticity) :

$$S_{X,Y}\lambda \mathrm{\Phi }.$$
(2)
We will then talk about “$`X`$ hybrid perturbation”.

## 3 An analytical estimate

Adiabatic scale invariant initial conditions make two predictions concerning the CMB anisotropies. First, they predict a flat (Sachs-Wolfe) plateau at low multipoles, which illustrates the fact that the gravitational potential is “frozen” as long as the modes have not yet entered into the Hubble radius. Second, one expects to find a serie of Doppler peaks at smaller angular scales, produced by acoustic oscillations in the photon-baryon plasma. The height of the first peak depends on almost all the cosmological parameters, but as soon as one considers adiabatic initial conditions and (very) conservative cosmological parameters, the peak is between $`2`$ and $`8`$ times higher than the Sachs-Wolfe plateau. This is a strong prediction of adiabatic models, and it is in good agreement with the current data. (In opposition, a pure isocurvature CDM model leads to a Sachs-Wolfe plateau higher than the first Doppler peak.) In the case of adiabatic and isocurvature mixtures, the CMB anisotropies and the matter power spectrum correspond, on large scales, to different combinations of the initial perturbations. Such a complementarity is all the more useful in our model that, contrarily to the adiabatic case, it does not generically predict the ratio of amplitude between the CMB and the matter power spectrum. Indeed, in the standard adiabatic case, it is a classic calculation to derive the temperature anisotropy as a function of the gravitational potential. In the long wavelength limit, the temperature anisotropies are one third of the gravitational potential, which has varied of a factor $`9/10`$ during the radiation-to-matter transition. In the case of CDM hybrid perturbations, this relation can be rewritten as :

$$\frac{\delta T}{T}|_{\mathrm{MD}}=\frac{3}{10}\left(1+\frac{4}{15}\mathrm{\Omega }_\nu ^{\mathrm{RD}}\frac{2}{5}\lambda \mathrm{\Omega }_c^{\mathrm{MD}}\right)\mathrm{\Phi }_{\mathrm{RD}},$$
(3)
where $`\mathrm{\Omega }_c`$ and $`\mathrm{\Omega }_\nu `$ are respectively the CDM and neutrino density parameters, and the indexes $`\mathrm{RD}`$ and $`\mathrm{MD}`$ mean that one considers the quantities during the radiation dominated and the matter dominated eras respectively. The relative amplitude of CMB anisotropies and matter power spectrum can therefore in principle be used to extract the isocurvature part of the initial conditions.

## 4 Numerical results

We have extensively studied the four hybrid perturbations in a recent paper $`^\mathrm{?}`$. The main result of our analysis is that the CMB anisotropy and the matter power spectra and rather strongly sensitive to the set of initial conditions we have considered. As a consequence, such models are already fairly well constrained : one cannot deviate strongly from adiabaticity. As an example, we have plotted on Fig. 1 the relative height of the first Doppler peak in the four hybrid models one can consider. Photon and CDM hybrid perturbations are the most constrained. As for the relative amplitude between the CMB anisotropies and the matter power spectrum, we have compared them on Fig. 2. It is clear that the amplitude of the first Doppler peak, which is already well measured by several ground and balloon experiments, is strongly sensitive to the parameter $`\lambda `$. This might help to discriminate between these models and the standard adiabatic model.

## References