text:
# Spatial Coherence Resonance near Pattern-Forming Instabilities

## Abstract

The analogue of temporal coherence resonance for spatial degrees of freedom is reported. Specifically, we show that spatiotemporal noise is able to optimally extract an intrinsic spatial scale in nonlinear media close to (but before) a pattern-forming instability. This effect is observed in a model of pattern-forming chemical reaction and in the Swift-Hohenberg model of fluid convection. In the latter case, the phenomenon is described analytically via an approximate approach.

The ability of noise to induce temporal coherence in nonlinear systems is a well documented fact. In stochastic resonance (SR), for instance, random fluctuations enhance the response to weak periodic driving, as observed in many different physical, chemical and biological scenarios wiesenfeld . From a more intriguing perspective, even systems showing no explicit time scale exhibit an enhancement of temporal coherence due to noise gang ; rappel . This autonomous stochastic resonance, where noise extracts a hidden, intrinsic time scale of the system’s dynamics, has been termed coherence resonance (CR) pikovsky , and has been predicted theoretically in a wide variety of models (see revexc  for a review) and observed experimentally in fields as diverse as plasma physics i , laser dynamics jorge , electronics postnov , atom trapping wilkowski  and neuroscience gu .

Parallel to the advances in the fields of stochastic and coherence resonance, recent years have witnessed an increasing interest in the dynamics of systems with spatial degrees of freedom. Much work has been devoted to analyze the enhancement of spatiotemporal order by noise nises ; alonso , and the phenomena of standard temporal SR and CR have been studied in detail in arrays of dynamical elements sr-ext ; cr-ext , where it has been shown that coupling enhances the noise-induced temporal coherence. But so far, up to our knowledge, no analogues of CR (or SR, for that matter) for spatial degrees of freedom exist. Namely, we ask ourselves if noise is able to extract and optimize an intrinsic spatial scale from a deterministically homogeneous extended medium. In this Letter we show that such spatial coherence resonance (SCR) exists close to pattern-forming instabilities, mimicking one of the mechanisms of standard (temporal) CR.

There are different mechanisms giving rise to temporal CR. In the most extended situation, CR arises as the result of a matching between the time scale associated to noise-induced escape from a stable steady state and the characteristic time scale of the phase-space trajectories outside from the immediate vicinity of the fixed point (such as what happens in excitable media pikovsky ; benji  or close to a saddle-node bifurcation in phase models rappel ). A second, less studied mechanism takes place near the onset of dynamical instabilities giving rise to periodic behavior, where noise excites precursors of the bifurcation in the form of stochastic limit cycles with a well-defined peak in the power spectrum, even before the bifurcation really happens kurt85 ; omberg . Such noisy precursors exhibit an optimal coherence at an intermediate noise strength, in a way typical of coherence resonance shura . This effect has been confirmed experimentally kiss .

Analogously to the temporal case, noisy precursors of pattern-forming bifurcations have been observed, in the form of noise-induced spatial structures, in electroconvection rehberg , Rayleigh-Bénard convection wu  and optics miguel , for instance. In this case, the spatial power spectrum (now structure function) of the system exhibits a well defined peak due to noise, even below the bifurcation point. In what follows, we show that noise controls the shape of this spectral peak in a manner characteristic of coherence resonance, but with the role of temporal frequency being played by spatial wave number. This allows us to define a pure spatial analogue of temporal coherence resonance. The generality of the phenomenon is illustrated through its occurrence in two different pattern-forming systems, namely an activator-inhibitor model of chemical pattern formation and the Swift-Hohenberg model of fluid convection. In the latter case the effect is interpreted analytically.

We first illustrate the phenomenon of SCR in chemical pattern formation, in particular in the case of the chlorine dioxide-iodine-malonic acid (CDIMA) reaction, which is well known to exhibit both temporal and spatial dynamics Eps98 . Experimental and theoretical studies of this reaction have concluded that the malonic acid serves only as a source of iodide, and thus the system can be described by just two variables, in the way of an activator-inhibitor model: Eps98 ; ancherlengues

$`\dot{u}`$ $`=`$ $`^2u+acu4{\displaystyle \frac{uv}{(1+u^2)}}\varphi `$
$`\dot{v}`$ $`=`$ $`b\left\{^2v+cu{\displaystyle \frac{uv}{(1+u^2)}}+\varphi \right\},`$ (1)
where $`u`$ and $`v`$ play the role of activator and inhibitor fields, corresponding to the concentrations of iodine and chlorine dioxide ions, respectively (all variables and parameters are in dimensionless units). This model describes the photosensitivity of the CDIMA reaction, with $`\varphi `$ accounting for the illumination level, which is an externally controlled parameter. Additionally, for large differences in the diffusion coefficients of the two species this system presents a Turing instability. In particular, the model has an homogenous stable state given by $`u_h=(a5\varphi )/(5c)`$ and $`v_h=a(1+u_h^2)/(5u_h)`$, which becomes unstable versus spatial perturbations of nonzero wavenumber for $`\varphi <\varphi _c`$. For $`a=16`$, $`b=301`$, $`c=0.6`$ and $`d=1.07`$, the critical illumination is $`\varphi _c=2.3`$ ancherlengues .

This chemical system allows for an experimental study of stochastic effects, by introducing noise in the illumination profile. In this way, for instance, it has been shown experimentally that a quenched spatial noise sustains Turing patterns in this system well beyond the deterministic bifurcation threshold ancherlengues . We now consider a spatiotemporally stochastic illumination profile, $`\varphi \varphi +\xi `$, where $`\xi (\stackrel{}{x},t)`$ is a Gaussian noise with zero mean, white in time and with spatial autocorrelation given by

$$\xi (\stackrel{}{x},t)\xi (\stackrel{}{x}^{},t^{})=\mathrm{\hspace{0.17em}2}\sigma ^2C(\stackrel{}{x}\stackrel{}{x}^{})\delta (tt^{}).$$
(2)
In order to examine the spatial response of the system in the presence of noise, we evaluate the structure function of the activator variable, defined as

$$S(\stackrel{}{k},t)=\frac{\widehat{u}(\stackrel{}{k},t)\widehat{u}(\stackrel{}{k},t)}{V},$$
(3)
where $`\widehat{u}(\stackrel{}{k},t)`$ is the spatial Fourier transform of the activator field, $`V`$ is the d-dimensional volume of the system, and $`\mathrm{}`$ represents an ensemble average over noise realizations. Taking into account the spherical symmetry of this system, we compute the spherical average of the structure function as $`s(k)=_{\mathrm{\Omega }_k}𝑑\mathrm{\Omega }_kS(\stackrel{}{k},t)`$, where $`k=\left|\stackrel{}{k}\right|`$ and $`\mathrm{\Omega }_k`$ is a hyperspherical shell of radius $`k`$. The result is shown in Fig. 1(a) for increasing noise intensities dclevel , evaluated for a two-dimensional system and in the case of spatially white noise, i.e. $`C(\stackrel{}{x}_i\stackrel{}{x}_j)=\delta _{ij}/\mathrm{\Delta }x`$. The average illumination $`\varphi `$ to which the system is subjected is above threshold, thus for vanishing or very weak noises the system is homogeneous and the structure function does not show a well-defined peak as $`k`$ increases. However, increasing the noise intensity induces a peak in the structure function at $`k=k_{\mathrm{max}}`$, which represents a spatial noisy precursor of the pattern-forming bifurcation that occurs at $`\varphi =\varphi _c`$: fluctuations are capable to anticipate the instability, making the structure manifest even before the threshold has been crossed.

If fluctuations become too large, disorder inevitably comes into play, and the peak in the structure function becomes less well defined \[see triangles in Fig. 1(a)\]. It is then clear that an optimal noise intensity exists at which the structure function peak is best resolved from the background fluctuations. In order to quantify the efficiency of the noise, an adequate measure must be used. There are different ways of quantifying such a response depending on the peculiarities of the problem under study. Here, we estimate the spatial coherence of the system through the quantity $`\mathrm{\Delta }Ss(k_{\mathrm{max}})s(0)`$, taking into account that $`s(0)`$ measures the level of fluctuations existing in the system.

Figure 1(b) shows how the SND varies with noise intensity $`\sigma ^2`$ for different illumination levels. It can be seen that there is an optimal noise strength $`\sigma _{\mathrm{opt}}^2`$ for which the structure function peak is resolved best. This is the signature of a spatial coherence resonance driven by the external stochastic forcing, which extracts the intrinsic spatial distribution of the system. As shown in Fig. 1(b), the closer the system is to the bifurcation point, the more pronounced the resonance is, and the position of the maximum $`\sigma _{\mathrm{opt}}^2`$ shifts to smaller values.

Figure 2 shows several patterns produced by model (1) for different values of the noise strength at $`\varphi =2.35`$. While the system is quasi-homogenous for small noise (left plot), it exhibits a rather self-organized configuration for $`\sigma ^2\sigma _{\mathrm{opt}}^2`$ (middle plot), and becomes gradually more disordered for large noise strengths (right plot). Although the effect may not appear very clear in the spatial profiles, the relevant fact is the appearance and enhancement of the peak in the structure function, since this is a measurable effect through, for instance, scattering experiments.

We may illustrate analytically this phenomenon in the case of a simpler pattern-forming model, namely the Swift-Hohenberg equation, which describes the onset of hydrodynamic convection in Rayleigh-Bénard cells cross :

$$\frac{\varphi (\stackrel{}{x},t)}{t}=r\varphi \varphi ^3\left(^2+k_0^2\right)^2\varphi +\xi (\stackrel{}{x},t),$$
(4)
where the control parameter $`r`$ is proportional to the temperature gradient driving the fluid, and the additive noise $`\xi (\stackrel{}{x},t)`$ is again considered Gaussian and white in both time and space. As in the CDIMA model, in the absence of noise and for $`r>0`$, this system has a homogeneous steady state $`\varphi =0`$, representing conduction, which upon decrease of $`r`$ becomes unstable at $`r=r_c=0`$ versus static perturbations of nonzero wavenumber $`k=k_0`$, which represents convection cross . In the presence of a small amount of noise, the structure function can be calculated by linearizing Eq. (4) around $`\varphi =0`$, and transforming the system to Fourier space nises . The result is:

$$S_{st}(k)=\frac{\sigma ^2}{(k^2k_0^2)^2+r}$$
(5)
Since we want to operate in regimes where, even though the system is deterministically homogeneous, noise can be large, we generalize the previous expression (5) by allowing the parameters $`\sigma `$ and $`r`$ to depend on noise; we will call them $`\sigma ^{}(\sigma )`$ and $`r^{}(\sigma )`$ in what follows. In order to estimate how these renormalized parameters depend on the noise intensity, we rederive expression (5) without neglecting completely nonlinear terms. To do that, we estimate the cubic term in (4) by applying the Gaussian approximation langer  $`\varphi (x,t)^33\varphi (x,t)\varphi ^2_{st}`$. This renders back the equation linear, so that the structure function obeys again expression (5) with $`r`$ being replaced by $`r^{}r+3\varphi ^2_{st}`$. We can now calculate the value of the average squared field $`\varphi ^2_{st}`$ from the self-consistency relation

$`\varphi ^2_{st}`$ $`=`$ $`{\displaystyle \frac{1}{(2\pi )^2}}{\displaystyle _0^{\mathrm{}}}S_{st}(k)2\pi k𝑑k`$ (6)
$`=`$ $`{\displaystyle \frac{\pi \sigma ^2}{(2\pi )^2\sqrt{r^{}}}}\left[{\displaystyle \frac{\pi }{2}}+tg^1\left({\displaystyle \frac{k_0^2}{\sqrt{r^{}}}}\right)\right]{\displaystyle \frac{\sigma ^2}{8\sqrt{r^{}}}}.`$
Since we will be operating close to threshold, $`r`$ will be small and $`r^{}3\varphi ^2_{st}`$, which leads from (6) to $`\varphi ^2_{st}(\sigma ^2)^{2/3}`$. Therefore, from the definition of $`r^{}`$ above we finally obtain that this effective parameter also scales with noise intensity as $`r^{}(\sigma ^2)^{2/3}`$. In order to verify this result, we perform numerical simulations of model (4) in the homogeneous regime for different noise intensities, and fit the resulting structure functions with expression (5), with effective parameters $`\sigma ^{}`$ and $`r^{}`$ which are allowed to depend on noise. The results of this analysis are shown in Fig. 3, and confirm that $`\sigma ^{}`$ scales linearly with noise intensity $`\sigma `$ (as expected, since nonlinear corrections did not affect this parameter in the analysis made above), whereas $`r^{}`$ scales with noise intensity with an exponent fairly close to 2/3.

The numerical results of Fig. 4(a) have been fitted with the analytic approximation given by Eq. (5)) with renormalized parameters, leading to a quite good agreement. This is shown in Fig. 4(a) for four different noise intensities. Similarly to what occurs in the chemical model, we can see how noise excites a peak in the structure function (in this case at $`k=k_0`$) even within the homogeneous regime. This peaks increases in size as noise intensity increases, both in an absolute way and relative to the value of the structure function at $`k=0`$. However, for large noise intensities disorder kicks in again, and the peak becomes less pronounced.

This effect can be seen quantitatively, in terms of the dependence of $`\mathrm{\Delta }S`$ on $`\sigma ^2`$, for two different distances to threshold in Fig. 4(b), where the numerical results are again compared to the theoretical prediction coming from the renormalized version of (5), which in this case gives

$$\mathrm{\Delta }S(\sigma ^2)=\frac{k_0^4\sigma ^2}{r^{}(k_0^4+r^{})}$$
(7)
The agreement between the numerical results and the theoretical prediction is rather satisfactory. In all cases, a clear enhancement of spatial coherence (in terms of spatial structures with wavenumber $`k_{\mathrm{max}}=k_0`$) is evident.

In conclusion, we have demonstrated that spatiotemporal noise is able to extract and enhance spatial coherence from deterministically homogeneous nonlinear media. This effect constitutes a pure spatial analogue of temporal coherence resonance, and relies on the generation of spatial noisy precursors of pattern-forming instabilities. This phenomenon is different from what happens in noise-induced pattern formation nipf , where special types of multiplicative noise induce (or displace) the pattern-forming transition, destabilizing the homogeneous state. Here that state is still stable, and the noise can be simply additive. Additionally, our results lead us to expect that spatial analogues of both stochastic resonance and other mechanisms of coherence resonance (such as the one that takes place in excitable systems pikovsky ) could be found. The latter perspective, given the ubiquity of excitable systems in all areas of science revexc , is in our opinion very attractive. Excitable neural tissue, for instance, combines the features of being highly noisy and intrinsically spatially extended kleinfeld . Elucidating whether noise is able to enhance the spatial coherence of the system in this context would be extremely interesting.

We thank L. Schimansky-Geier for useful comments. This research was supported by the Ministerio de Ciencia y Tecnología (Spain) and FEDER under projects BFM2000-0624, BFM2001-2159, and BFM2002-04369. J.G.O. is partially supported by the NSF IGERT Program of Nonlinear Sciences (Cornell).


format: text

text:
# 1 Introduction

## 1 Introduction

The Standard Model (SM) has been highly successful in explaining all experimental observations in the energy regime up to a few hundred GeV. However, it is believed to be an effective field theory valid only up to a cutoff scale $`\mathrm{\Lambda }`$. Non-renormalizable operators which are gauge invariant but suppressed by appropriate inverse powers of $`\mathrm{\Lambda }`$ should then be considered in the low energy effective theory . The dimension 5 operator $`\mathrm{}\mathrm{}HH/\mathrm{\Lambda }_L`$ ($`\mathrm{}`$ is the lepton doublet) which violates lepton number ($`L`$) by two units is the lowest dimensional of such operators. Experimental evidence for neutrino masses suggests the effective scale of $`L`$-violation is around $`\mathrm{\Lambda }_L10^{14}10^{15}\mathrm{GeV}`$. The $`d=6`$ operator $`QQQ\mathrm{}/\mathrm{\Lambda }_B^2`$ violates both baryon number ($`B`$) and lepton number and leads to the decay of the nucleon. The current limits on proton lifetime are $`\tau _p>5\times 10^{33}\mathrm{yrs}`$ for $`pe^+\pi ^0`$ and $`\tau _p>1.6\times 10^{33}\mathrm{yrs}`$ for $`p\overline{\nu }K^+`$ . These limits imply that $`\mathrm{\Lambda }_B>10^{15}\mathrm{GeV}`$. Grand Unified Theories with or without supersymmetry generate such $`B`$-violating operator with $`\mathrm{\Lambda }_B10^{14}10^{16}\mathrm{GeV}.`$ These theories are currently being tested through nucleon decay. Any new physics with a threshold $`\mathrm{\Lambda }`$ less than the GUT scale will thus be constrained by both proton lifetime and neutrino masses .

It is widely anticipated that new physics will show up around the TeV scale. This expectation is based on the stability of the Higgs mass. Supersymmetric theories with threshold of a few hundred GeV have been extensively studied as a solution to the Higgs mass problem (or the hierarchy problem). Indeed, these theories allow for both $`B`$ and $`L`$ violating operators. Care must be (and usually is) taken to satisfy the experimental constraints.

A second scenario that has been widely discussed recently to address the hierarchy problem is the possible existence of large extra dimensions . Here the fundamental scale of gravity $`M_{\mathrm{Pl}}`$ is lowered from its 4 dimensional value of $`10^{19}\mathrm{GeV}`$ to $`M_{}`$ a few TeV. Non-renormalizable operators associated with quantum gravity which are only suppressed by inverse power of $`M_{}`$ will then be present in the low energy theory. $`L`$ and $`B`$ violation would put major constraints on these theories. Similar remarks will apply to other solutions to the hierarchy problem, such as composite Higgs, little Higgs , etc.

In this paper we address the role discrete gauge symmetries  can play in suppressing baryon number and lepton number violating processes even when the scale of new physics is low. The SM effective lagrangian does not have a continuous anomaly-free symmetry that can suppress these processes. This is our reason for focusing on discrete symmetries. It is preferable that such symmetries have a gauge origin since all global symmetries are expected to be violated by the quantum gravitational effects . We find that the SM lagrangian, including small neutrino masses, has an anomaly-free $`Z_6`$ symmetry. This $`Z_6`$ acts as discrete baryon number and forbids potentially dangerous $`\mathrm{\Delta }B=1`$ and $`\mathrm{\Delta }B=2`$ nucleon decay processes. $`\mathrm{\Delta }B=3`$ operators are allowed, but they arise only as $`d=15`$ operators. “Triple nucleon decay” can proceed through these operators; we estimate its rate and find $`\mathrm{\Lambda }_B`$ can be as low as $`10^2\mathrm{GeV}`$. We show that the $`Z_6`$ symmetry has a natural gauge origin in the anomaly-free $`U(1)`$ corresponding to $`I_R^3+L_i+L_j2L_k`$, where $`L_i`$ is $`i`$th lepton number and $`I_R^3`$ is the third component of the righthanded isospin. No new particles are introduced to cancel anomalies except the righthanded neutrinos. We show how an additional $`Z_5`$ or $`Z_7`$ gauge symmetry generates the desired neutrino masses and mixings even with a low scale for $`\mathrm{\Lambda }`$.

Discrete gauge symmetries have been utilized in suppressing nucleon decay  as well as in addressing other aspects of physics such as solving the $`\mu `$ problem of supersymmetry, fermion mass hierarchy problem  and the stability of the axion . A $`Z_3`$ baryon parity was found in Ref.  that suppresses nucleon decay. In order for it to have a gauge origin, extra particles were introduced.444The charges of the SM particles under the $`U(1)`$ symmetry in Ref.  are $`\{Q,u,d,l,e\}=(0,1,1,1,1)`$. The model has extra leptons $`l^{}(3)`$, $`\overline{l}^{}(0)`$, $`\overline{e}^{}(3)`$, $`\overline{e}^{}(0)`$ (where the $`U(1)`$ charges are given in parentheses) as well as exotic leptons with electric charge of $`\pm 1/3`$ and $`0`$: $`2\times F^{1/3}(3)`$, $`1\times F^{1/3}(3)`$, $`2\times F^{1/3}(0)`$, $`1\times F^{1/3}(0)`$ and $`5\times N^0(3)`$, $`N^0(1)`$, $`N^0(3)`$. There have been also attempts to gauge baryon number by introducing a fourth family of fermions , or with the use of extra dimensions . The spirit of this paper is similar to that of Ref. , our models are however much more economical.

## 2 Gauged Baryon Parity and Suppression of Proton Decay

In this section we show that the SM lagrangian has a discrete $`Z_6`$ gauge symmetry which forbids all $`\mathrm{\Delta }B=1`$ and $`\mathrm{\Delta }B=2`$ baryon violating effective operators. This can be seen as follows. The SM Yukawa couplings incorporating the seesaw mechanism to generate small neutrino masses is

$$_\mathrm{Y}=Qu^cH+Qd^cH^{}+\mathrm{}e^cH^{}+\mathrm{}\nu ^cH+M_R\nu ^c\nu ^c.$$
(1)
Here we have used the standard (lefthanded) notation for the fermion fields and have not displayed the Yukawa couplings or the generation indices. This lagrangian respects a $`Z_6`$ discrete symmetry with the charge assignment as shown in Table 1. Also shown in Table 1 are the charge assignments under the $`Z_3`$ and $`Z_2`$ subgroups of $`Z_6`$. The $`Z_3`$ assignment is identical to that in Ref.

From Table 1 it is easy to calculate the $`Z_6`$ crossed anomaly coefficients with the SM gauge groups. We find the $`SU(3)_C`$ or $`SU(2)_L`$ anomalies to be

$`A_{[SU(3)_C]^2\times Z_6}=3N_g`$
$`A_{[SU(2)_L]^2\times Z_6}=N_g`$ (2)
where $`N_g`$ is the number of generations. The condition for a $`Z_N`$ discrete group to be anomaly-free is

$$A_i=\frac{N}{2}modN$$
(3)
where $`i`$ stands for $`SU(3)_C`$ and $`SU(2)_L`$. For $`Z_6`$, this condition reduces to $`A_i=3mod6`$, so when $`N_g=3`$, $`Z_6`$ is anomaly-free. Obviously, the $`Z_3`$ and $`Z_2`$ subgroups are also anomaly-free. The significance of this result is that unknown quantum gravitational effects will respect this $`Z_6`$. It is this feature that we utilize to stabilize the nucleon. Absence of anomalies also suggests that the $`Z_6`$ may have a simple gauge origin.

To see how the $`Z_6`$ forbids $`\mathrm{\Delta }B=1`$ and $`\mathrm{\Delta }B=2`$ processes, we note that it is a subgroup of $`U(1)_{2YB+3L}`$ where $`Y`$ is SM hypercharge . We list in Table 2 the charges under the three $`U(1)`$ symmetries. Comparing the last line of Table 2 with the charge assignment of Table 1 it is clear that the $`Z_6`$ is a subgroup of $`U(1)_{2YB+3L}`$. Any $`Z_6`$ invariant effective operator must then satisfy

$$2\mathrm{\Delta }Y\mathrm{\Delta }B+3\mathrm{\Delta }L=0mod6.$$
(4)
Invariance under $`U(1)_Y`$ implies $`\mathrm{\Delta }Y=0`$. Consider $`\mathrm{\Delta }B=1`$ effective operators which must then obey (from Eq. (4)) $`3\mathrm{\Delta }L=1mod6`$. This has no solution, since $`3\mathrm{\Delta }L=0mod3`$ from Table 2. Similarly, $`\mathrm{\Delta }B=2`$ operators must obey $`3\mathrm{\Delta }L=2mod6`$ which also has no solution. $`\mathrm{\Delta }B=3`$ operators, which corresponds to $`3\mathrm{\Delta }L=0mod6`$, are allowed by this $`Z_6`$. Such operators have dimension 15 or higher and have suppression factors of at least $`\mathrm{\Lambda }^{11}`$. These will lead to “triple nucleon decay” processes where three nucleons in a heavy nucleus undergo collective decays leading to processes such as $`pnne^+\pi ^0`$. We estimate the rates for such decay in Section 4 and find that $`\mathrm{\Lambda }`$ can be as low as $`10^2\mathrm{GeV}`$.

## 3 Embedding $`Z_6`$ in $`I_R^3+L_i+L_j2L_k`$

It is interesting to see if the $`Z_6`$ symmetry of Table 1 can be realized as an unbroken subgroup of a gauged $`U(1)`$ symmetry. Although the $`Z_6`$ is a subgroup of the $`U(1)_{(2YB+3L)}`$, this $`U(1)`$ would be anomalous without enlarging the particle content. We have found a simple and economic embedding of $`Z_6`$ into a $`U(1)`$ gauge symmetry associated with $`I_R^3+L_i+L_j2L_k`$. Here $`L_i`$ is the $`i`$th family lepton number and $`ijk`$. No new particles are needed to cancel gauge anomalies. With the inclusion of righthanded neutrinos $`I_R^3=Y(BL)/2`$ is an anomaly-free symmetry. $`L_i+L_j2L_k`$, which corresponds to the $`\lambda _8`$ generator acting in the leptonic $`SU(3)`$ family space, is also anomaly-free.

The charges of the SM particles under this $`U(1)`$ are

$$\begin{array}{cccc}Q_i=(0,0,0),& u_{i}^{}{}_{}{}^{c}=(1,1,1),& d_{i}^{}{}_{}{}^{c}=(1,1,1),& \\ \mathrm{}_i=(4,2,2),& e_{i}^{}{}_{}{}^{c}=(5,1,1),& \nu _{i}^{}{}_{}{}^{c}=(3,3,3),& H=1.\end{array}$$

This charge assignment allows all quark masses and mixings as well as charged lepton masses. When the $`U(1)`$ symmetry breaks spontaneously down to $`Z_6`$ by the vacuum expectation value of a SM singlet scalar field $`\varphi `$ with a charge of 6, realistic neutrino masses and mixings are also induced. The relevant lagrangian for the righthanded neutrino Majorana masses is

$$_{\nu ^\mathrm{c}}=M_{12}\nu _1^c\nu _2^c+M_{13}\nu _1^c\nu _3^c+\nu _3^c\nu _3^c\varphi +\nu _2^c\nu _2^c\varphi +\nu _1^c\nu _1^c\varphi ^{}+\nu _3^c\nu _2^c\varphi .$$
(5)
After integrating out the heavy righthanded neutrinos we obtain the following $`\mathrm{\Delta }L=2`$ effective operators:

$$_{\mathrm{\Delta }L=2}=\frac{1}{\mathrm{\Lambda }}(\mathrm{}_1\mathrm{}_2HH+\mathrm{}_1\mathrm{}_3HH+\mathrm{}_1\mathrm{}_1HHϵ+\mathrm{}_2\mathrm{}_2HHϵ^{}+\mathrm{}_2\mathrm{}_3HHϵ^{}+\mathrm{}_3\mathrm{}_3HHϵ^{}).$$
(6)
Here $`\mathrm{\Lambda }M_{12}M_{13}`$ is the scale of $`L`$-violation and we have defined $`ϵ\varphi /\mathrm{\Lambda }`$. For $`ϵ1`$, this lagrangian leads to the inverted mass hierarchy pattern for the neutrinos which is well consistent with the current neutrino oscillation data. This neutrino mass mixing pattern is analogous to the one obtained from $`L_eL_\mu L_\tau `$ symmetry . However, here the $`U(1)`$ is a true gauge symmetry.

We have also investigated other possible $`U(1)`$ origin of the $`Z_6`$ symmetry and found the $`I_R^3+L_i+L_j2L_k`$ combination to be essentially unique. To see this, let us assign a general $`U(1)_X`$ charge for $`i`$th generation of the SM fermions consistent with the $`Z_6`$ symmetry as

$$\{Q_i,u_i^c,d_i^c,\mathrm{}_i,e_i^c,\nu _i^c\}=\{6m_1^{(i)},5+6m_4^{(i)},1+6m_3^{(i)},2+6m_2^{(i)},5+6m_5^{(i)},3+6m_6^{(i)}\}$$

where $`m_j^{(i)}`$ are all integers. The Higgs field has a charge $`H=1+6m_0`$. If we impose the invariance of the Yukawa couplings of the charged fermions and Dirac neutrinos for each generation, the anomaly coefficients from $`i`$th generation become

$`A_{[SU(3)_C]^2\times U(1)_X}^{(i)}`$ $`=`$ $`0`$
$`A_{[SU(2)_L]^2\times U(1)_X}^{(i)}`$ $`=`$ $`1+9m_1^{(i)}+3m_2^{(i)}`$
$`A_{[U(1)_Y]^2\times U(1)_X}^{(i)}`$ $`=`$ $`(1+9m_1^{(i)}+3m_2^{(i)})`$
$`A_{[U(1)_X]^2\times U(1)_Y}^{(i)}`$ $`=`$ $`[5+m_0]A_{[SU(2)_L]^2\times U(1)_X}^{(i)}`$
$`A_{[U(1)_X]^3}^{(i)}`$ $`=`$ $`[5+m_0]^2A_{[SU(2)_L]^2\times U(1)_X}^{(i)}.`$ (7)
The coefficient for the mixed gravitational anomaly for each generation is zero. From Eq. (3), it follows that $`A_2=_iA_{[SU(2)_L]^2\times U(1)_X}^{(i)}=_i(1+9m_1^{(i)}+3m_2^{(i)})=0`$ can be satisfied only when all three generation contributions are included. Once $`A_2=0`$ is satisfied, all other anomaly coefficients will automatically vanish. $`A_2=0`$ can be rewritten in a familiar form as $`3_iQ_i+_i\mathrm{}_i=0`$. Thus we see that any $`U(1)`$ symmetry satisfying this condition and consistent with the $`Z_6`$ charge assignment can be a possible source of $`Z_6`$. If the $`Q_i`$ are different for different generations, quark mixings cannot be generated without additional particles. By making a shift proportional to hypercharge, we can set $`Q_i=0`$ for all $`i`$. Two obvious solutions to $`_i\mathrm{}_i=0`$ are $`\mathrm{}_i=(1,1,2)`$ and $`\mathrm{}_i=(1,1,0)`$. The latter one does not reproduce the $`Z_6`$ charge assignment while the former one does, which is our solution when $`I_R^3`$ is added to it.

A related $`B3L_\tau `$ has been discussed in Ref. . This is the same as $`BL`$ plus $`L_e+L_\mu 2L_\tau `$. In Ref , only one righthanded neutrino $`\nu _\tau ^c`$ is introduced so the seesaw mechanism applies only for one light neutrino. The other two neutrinos receive small masses from radiative correction. In our model, since there are three righthanded neutrinos, all the neutrino masses arise from the conventional seesaw mechanism.

## 4 $`\mathrm{\Delta }B=3`$ Operators and Triple Nucleon Decay

We now list the lowest dimensional (d=15) $`\mathrm{\Delta }B=3`$ effective operators which are consistent with the $`Z_6`$ symmetry. Imposing gauge invariance and Lorentz invariance, we find them to be:

$`\overline{u^c}^4\overline{d^c}^5\overline{e^c},\overline{u^c}^2\overline{d^c}^7e^c,Q\overline{u^c}^3\overline{d^c}^5\mathrm{},Q\overline{u^c}^2\overline{d^c}^6\overline{\mathrm{}},Q^2\overline{u^c}^3\overline{d^c}^4\overline{e^c},`$
$`Q^2\overline{u^c}\overline{d^c}^6e^c,Q^3\overline{u^c}^2\overline{d^c}^4\mathrm{},Q^3\overline{u^c}\overline{d^c}^5\overline{\mathrm{}},Q^4\overline{u^c}^2\overline{d^c}^3\overline{e^c},Q^4\overline{u^c}\overline{d^c}^4\nu ^c,`$
$`Q^4\overline{d^c}^5e^c,Q^5\overline{u^c}\overline{d^c}^3\mathrm{},Q^5\overline{d^c}^4\overline{\mathrm{}},Q^6\overline{u^c}\overline{d^c}^2\overline{e^c},Q^7\overline{d^c}^2\mathrm{},Q^8\overline{d^c}\overline{e^c}.`$ (8)
Here Lorentz, gauge and flavor indices are suppressed. These operators can lead to “triple nucleon decay”. The dominant processes are

$`ppp`$ $``$ $`e^++\pi ^++\pi ^+`$
$`ppn`$ $``$ $`e^++\pi ^+`$
$`pnn`$ $``$ $`e^++\pi ^0`$
$`nnn`$ $``$ $`\overline{\nu }+\pi ^0.`$ (9)
Tritium ($`{}_{}{}^{3}H`$) and Helium-3 ($`{}_{}{}^{3}He`$) are examples of three-nucleon systems in nature. These nuclei are unstable and undergo $`\beta `$-decay with relatively short lifetime. In the presence of operators of Eq. (4), $`{}_{}{}^{3}He^++\pi ^0`$ and $`{}_{}{}^{3}Hee^++\pi ^+`$ decays can occur. However, there is no stringent experimental limits arising from these nuclei. So we focus on triple-nucleon decay in Oxygen nucleus where there are experimental constraints from water detectors. To estimate the decay lifetime we need to first convert the nine-quark operators of Eq. (4) into three-nucleon operators and subsequently into Oxygen nucleus.

We choose a specific operator $`Q^5\overline{d^c}^4\overline{\mathrm{}}/\mathrm{\Lambda }^{11}`$ as an example to study the process $`pnne^++\pi ^0`$ triple nucleon decay process. This induces the effective three-nucleon operator in Oxygen nucleus

$$\frac{Q^5\overline{d^c}^4\overline{\mathrm{}}}{\mathrm{\Lambda }^{11}}\frac{\beta ^3(1+D+F)}{\sqrt{2}f_\pi \mathrm{\Lambda }^{11}}(\pi nnpe),$$
(10)
where $`\beta 0.014GeV^3`$ is the matrix element to convert three quarks into a nucleon . $`F0.47`$, $`D0.80`$ are chiral lagrangian factors, and $`f_\pi =139`$ MeV is the pion decay constant.

We now estimate the wave-function overlap factor for three nucleons inside Oxygen nucleus to find each other. This is based on a crude free fermi gas model where the nucleons are treated as free particles inside an infinite potential well. A single nucleon wave function is given by $`\psi _m(x)=\sqrt{2/r}\mathrm{sin}(m\pi x/r)`$, where $`r`$ is the size of nucleus and $`m`$ is the energy level. Incorporating isospin and Pauli exclusion principle, the highest energy level which corresponds to $`m=4`$ is found to have 2 protons and 2 neutrons. We assume the highest level has the most probability to form a Tritium-like “bound state” of three nucleons. The probability for three nucleons in Oxygen nucleus to overlap in a range the size of Tritium nucleus is

$$P\frac{4}{3}_0^{\sqrt[3]{\frac{3}{16}}}d\left(\frac{x_1}{r}\right)d\left(\frac{x_2}{r}\right)d\left(\frac{x_3}{r}\right)\left(\mathrm{sin}\left(\frac{4\pi x_1}{r}\right)\mathrm{sin}\left(\frac{4\pi x_2}{r}\right)\mathrm{sin}\left(\frac{4\pi x_3}{r}\right)\right)^20.0253,$$
(11)
where $`\sqrt[3]{\frac{3}{16}}`$ is the ratio between the radii of Tritium nucleus and Oxygen nucleus, since $`RA^{1/3}`$ (A is the atomic number). So the effective baryon number violating operator of Eq. (10) becomes

$$\frac{P\beta ^3}{\sqrt{2}f_\pi \mathrm{\Lambda }^{11}R^3}({}_{}{}^{3}H\pi e).$$
(12)
The triple nucleon decay lifetime can then be estimated to be

$$\tau \frac{16\pi f_{\pi }^{}{}_{}{}^{2}\mathrm{\Lambda }^{22}R^6}{P^2\beta ^6M_{{}_{}{}^{3}H}}.$$
(13)
By putting the current limit on proton lifetime of $`3\times 10^{33}\mathrm{yrs}`$, we obtain:

$$\mathrm{\Lambda }10^2\mathrm{GeV}.$$
(14)
Thus we see the $`Z_6`$ symmetry ensures the stability of the nucleon.

To test our crude model of nuclear transition, we have also evaluated the double nucleon decay rate within the same approach and found our results to be consistent with other more detailed evaluations .

## 5 Realistic Neutrino Masses with Low Threshold

If the threshold of new physics is low, but is still much larger than a TeV, small neutrino masses can be generated by the seesaw mechanism as discussed in Section 3 through the effective operators $`\mathrm{}\mathrm{}HH/\mathrm{\Lambda }`$. If the threshold is as low as a few TeV, the induced neutrino mass will be too large. Here we show a mechanism by which such operators can be suppressed by making use of a discrete $`Z_N`$ symmetry (with $`N`$ odd) surviving to low scale. This $`Z_N`$ has a natural embedding in $`BL`$ gauge symmetry.

Consider the following effective operators in the low energy lagrangian:

$$\mathrm{}\mathrm{}HH\frac{S^6}{\mathrm{\Lambda }^7}+\frac{S^{2N}}{\mathrm{\Lambda }^{2N4}}.$$
(15)
Here $`S`$ is a singlet field which has charge $`(1,3)`$ under $`Z_N\times Z_6`$ while $`\mathrm{}`$ has charge $`(3,2)`$. (The $`Z_6`$ charges of SM particles are as listed in Table 1.) The first term in Eq. (15) respects a $`U(1)`$ symmetry while the second term reduces this to $`Z_6\times Z_N`$. If $`S`$ develops a VEV of order $`10^2\mathrm{GeV}`$, realistic neutrino mass can arise even when $`\mathrm{\Lambda }`$ is low. For example, if $`\mathrm{\Lambda }=10\mathrm{TeV}`$ and $`S=10^2\mathrm{GeV}`$, the neutrino mass is of order $`v^2S^6/\mathrm{\Lambda }^70.4\mathrm{eV}`$, which is consistent with the mass scale suggested by the atmospheric neutrino oscillation data.

Two explicit examples of the $`Z_N`$ symmetry with $`N=5`$ and $`7`$ are shown in Table 3. These $`Z_N`$ symmetries are free from gauge anomalies. In the $`Z_5`$ example, the crossed anomaly coefficients for $`SU(3)_C`$ and $`SU(2)_L`$ are $`5N_g`$ and $`5N_g/2`$ respectively showing that $`Z_5`$ is indeed anomaly-free. For $`Z_7`$, these coefficients are $`7N_g`$ and $`7N_g/2`$, so it is also anomaly-free.

It is interesting to ask if the $`Z_N`$ can be embedded into a gauged $`U(1)`$ symmetry. A simple possibility we have found is to embed this $`Z_N`$ into the anomalous $`U(1)_A`$ symmetry of string origin with the anomalies cancelled by the Green-Schwarz mechanism . Consider $`U(1)_{BL}`$ without the right handed neutrinos but with the inclusion of vector-like fermions which have the quantum numbers of $`\mathrm{𝟓}(3)`$ and $`\overline{\mathrm{𝟓}}(2)`$ under $`SU(5)\times U(1)_A`$. This $`U(1)_A`$ is anomaly-free by virtue of the Green-Schwarz mechanism. When this $`U(1)_A`$ breaks down to $`Z_5`$, the extra particles get heavy mass and are removed from the low energy theory which is the $`Z_6\times Z_5`$ model.

Without the second term in Eq. (15), the phase of $`S`$ field will be massless upon spontaneous symmetry breaking. This Majoron field  would however acquire a mass from the second term of Eq. (15). In the $`Z_6\times Z_5`$ model, the mass of the Majoron is of order $`S^7/\mathrm{\Lambda }^6100\mathrm{keV}`$. In the $`Z_6\times Z_7`$ model, the Majoron mass is of order $`S^{11}/\mathrm{\Lambda }^{10}10\mathrm{eV}`$. Such a Majoron with a mass of either $`100\mathrm{keV}`$ or $`10\mathrm{eV}`$ is fully consistent with constraints from early universe cosmology . The interaction term $`\mathrm{}\mathrm{}HHS^6/\mathrm{\Lambda }^7`$ induces the Majoron decay $`S\nu \nu `$ with a Yukawa coupling $`Y_{S\nu \nu }=6m_\nu /S10^{11}`$. The decay rate of the Majoron can be estimated to be

$$\mathrm{\Gamma }=\frac{Y_{S\nu \nu }^2m_S}{8\pi }10^{23}m_S.$$
(16)
This corresponds to a Majoron lifetime of $`\tau 10\mathrm{sec}`$ for the $`Z_6\times Z_5`$ model and $`\tau 10^5\mathrm{sec}`$ for the $`Z_6\times Z_7`$ model. Such a Majoron can modify the big-bang nucleosynthesis processes. However the modification is not significant since the Majoron will decouple before the electro-weak phase transition. Its contribution to the expansion rate is equivalent to that of $`0.047\times 4/70.027`$ light neutrino species . This extra contribution is well within observational uncertainties.

## 6 Conclusions

We have shown in this paper that the SM lagrangian including small neutrino masses and with three generations has an anomaly-free $`Z_6`$ discrete symmetry. This $`Z_6`$ can act as gauged baryon parity and stabilize the nucleon even when the new physics threshold is below the TeV scale. The $`Z_6`$ has a natural embedding in $`I_R^3+L_i+L_j2L_k`$ gauge symmetry. We have also shown that realistic neutrino masses can be obtained through an additional anomaly-free $`Z_5`$ or $`Z_7`$ symmetry which may have its origin in the anomalous $`U(1)_A`$ symmetry.

Our framework can be readily extended to incorporate supersymmetry. The $`Z_6`$ charge assignment is still as in Table 1. The down-type Higgs $`\overline{H}`$ has a $`Z_6`$ charge of $`5`$ and the Grassmann variable $`\theta `$ has charge $`3`$ or $`0`$ . The Higgsino contribution to the $`[SU(2)_L]^2\times Z_6`$ anomaly is $`3`$ and the gaugino contribution to $`[SU(2)_L]^2\times Z_6`$ and $`[SU(3)_C]^2\times Z_6`$ anomalies are $`6`$ and $`9`$ when $`\theta `$ charge is not zero. Since these contributions are all multiples of $`6/2`$, the $`Z_6`$ symmetry remains anomaly-free in supersymmetric case. This $`Z_6`$ forbids all $`\mathrm{\Delta }B=1`$ and $`\mathrm{\Delta }B=2`$ nucleon decay as in the non-supersymmetric case.

## 7 Acknowledgement

We are grateful to Yukihiro Mimura for useful discussion. This work is supported in part by DOE Grant # DE-FG03-98ER-41076, a grant from the Research Corporation and by DOE Grant # DE-FG02-01ER-45684.


text:
# Ferroelectric ordering and electroclinic effect in chiral smectic liquid crystals

## I Introduction

Ferroelectric liquid crystals (FLC) have been studied intensively for several decades. Chiral smectic liquid crystals (LC) have unique material properties like spontaneous polarization M . Although a single chiral smectic LC molecule has a non-zero dipole moment due to symmetry considerations M1 , LC molecules in the bulk of the sample tend to form a helical structure which leads to polarization compensation. The helix can be unwound by application of a strong electric field or using thin test cells. Conventionally, alignment of LC is usually obtained by unidirectional mechanical rubbing of a thin polymer layer which coats the inner cell surfaces. Thus, the LC molecules point their long axes along the rubbing direction. This anisotropic interaction influences the ferroelectric ordering and the switching behavior which are of great practical importance for any LC device. The ferroelectrically ordered SmC\* phase is characterized by a non zero angle between the molecular long axis and the smectic layer normal. This tilt angle is the order parameter used to describe the second order SmC\*-SmA\* phase transition.

In several studies the temperature dependence of the tilt angle below the critical temperature follows a power law of 0.5 (classical behavior) M4 , whereas in other studies it follows a power law of 0.3, in strong analogy with the superfluid helium predicted by de Gennes M5 . The SmA\* phase in the vicinity of the transition point has also been the subject of intensive investigations, since in that part of the critical region the electroclinic (EC) effect can induce a tilt angle which is proportional to the electric field strength M7 . The induced tilt angle is also strongly temperature dependent and diverges as the phase transition temperature is approached from above. Another type of this phenomenon, the so-called surface EC effect, is observed in the interfacial region of the FLC cells and originates from the interaction of the sub-surface layers of chiral smectic molecules and a localized surface field ShenXue . In the vicinity of the transition point, the critical exponents of the tilt angle are in the range of 0.5 $`÷`$ 1.5 ShenXue ; M5 ; M7 ; Selinger . The surface EC effect results in the existence of sub-surface region in which the director is twisted from the rubbing axis to the bulk alignment direction. The ferroelectric properties of these twisted layers are of great importance as they play a dominant role in LC-devices.

In this paper we use linear electrooptic (EO) and second harmonic generation (SHG) techniques for investigation of ferroelectric ordering, SmC\*-SmA\* phase transitions and the EC effect in thin FLC cells. The main emphasis is made on comparative studies of the SHG reflection and transmission experimental geometries in order to figure out the roles of the sub-surface layers and the bulk in ferroelectric phase transitions. The SHG method is well-known for its sensitivity to symmetrical, structural and electronic properties of surfaces, interfaces and ultrathin films, and has been widely used along with EO for studying ferroelectric properties of chiral smectic LC shgeo . As the SHG response strongly depends on the polar state of the matter because of its unique sensitivity to the breakdown of the inversion symmetry M13 , it is a powerful instrument for probing ferroelectric phase transitions and electric-field induced effects in the vicinity of the critical region.

## II Model description

We studied linear and nonlinear quadratic responses of the FLC cells upon application of an electric field and temperature variation up to the SmC\*-SmA\* phase transition. Chiral smectic LC belong to the $`C_2`$ point group symmetry with the symmetry axes oriented in the cell plane. A schematic representation of the FLC cell and the coordinate system is shown in fig. 1. In this coordinate frame $`Oz`$ is the normal to the cell, $`XY`$ is the cell plane and the $`C_2`$ axis is parallel to the $`Oy`$ direction and coincide with the main optical axis. We assume the FLC cell to be uniaxial, with the ordinary and extraordinary refractive indices $`n_0=n_x=n_z`$, and $`n_e=n_y`$, respectively. The transmittance of an optically active plate of the thickness $`d`$ and birefringence $`\mathrm{\Delta }n`$ at the wavelength $`\lambda `$ in crossed-polarizers geometry is determined by the equation:

$$T=\mathrm{sin}^2(2\alpha )\mathrm{sin}^2\frac{\mathrm{\Delta }n\pi d}{\lambda },$$
(1)
where $`\alpha `$ is the angle between the main optical axes and polarization direction of the incident light. The linear electro-optic tensor $`r_{ijk}`$, where indices i, j, k denote the axes of the cell coordinate system, for the given symmetry has 8 non-vanishing components $`r_{xxy}`$, $`r_{yyy}`$, $`r_{zzy}`$, $`r_{yzx}`$, $`r_{yzz}`$, $`r_{xzy}`$, $`r_{xyx}`$, $`r_{xyz}`$ M16 . Application of the DC electric field E along $`Oz`$ direction leads to the rotation of the optical axes, which could be interpreted in terms of the rotation of the ellipse of refractive indexes. The equation of index ellipsoid in the presence of the electric field $`E`$ has the form

$$\frac{x^2}{n_x^2}+\frac{y^2}{n_y^2}+\frac{z^2}{n_z^2}+2yzr_{yzz}E+2xyr_{xyz}E=1.$$
(2)
The transformation of (2) to the canonical form gives the DC field-induced rotation of the undisturbed system of coordinates based on the main optical axes. Considering the rotation only around $`Oz`$-axis, the influence of the electric field can be described by the tilt of the main optical axis in the cell plane with the dependence of the turn angle $`\mathrm{\Delta }\alpha `$ on the electric field as:

$$\mathrm{\Delta }\alpha =\frac{1}{2}\mathrm{arctan}\frac{2r_{xyz}E}{\frac{1}{n_o^2}+\frac{1}{n_e^2}}.$$
(3)
The nonlinear quadratic susceptibility tensor of the FLC cell with a $`C_2`$ symmetry is given by following non-vanishing components $`\chi _{ijk}`$, $`\chi _{yxx}`$, $`\chi _{yyy}`$, $`\chi _{yzz}`$, ,$`\chi _{yzx}`$, $`\chi _{xyz}`$, $`\chi _{xxy}`$, $`\chi _{zyz}`$, $`\chi _{zxy}`$ M17 . The contribution to the SHG intensity of each quadratic susceptibility component depends on the angle of incidence of the fundamental radiation, on the azimuthal position of the sample and polarizations of the input and SHG light. In transmission geometry at normal incidence only three components participate in SHG - $`\chi _{yxx}`$, $`\chi _{xxy}`$ and $`\chi _{yyy}`$, and for pp -geometry the SHG intensity can be expressed by:

$$I_{pp}^{2\omega }\mathrm{cos}^2\theta (\chi _{yyy}\mathrm{cos}^2\theta +(2\chi _{xxy}+\chi _{yxx})\mathrm{sin}^2\theta )^2,$$
(4)
where the azimuthal angle $`\theta `$ is the angle between the direction of the p-polarization and $`C_2`$ axes of the cell symmetry. By fitting the anisotropy dependences of the SHG intensity in all polarization combinations, the corresponding components of the quadratic susceptibility can be extracted.

According to a well-known model, in the smectic planar layers FLC molecules can precess on the surface of the smectic cone cones18 . Ferroelectric switching is attributed to the interaction between molecular dipoles and external electric field and results in rotation of the molecules and, correspondingly, of the symmetry axes, within the half of the smectic cone. If the molecular director positions, corresponding to the saturating electric fields of opposite values, are in the cell plane and on the opposite sides of the smectic cone, then (3) can be used for the explanation of linear and nonlinear switching, as it shows the in-plane rotation of the main optical axis.

The angle between the molecular long axis and the normal to the smectic layers, i.e. the apex angle of the cone, is the order parameter of the SmC\*-SmA\* phase transition, vanishing in the SmA\* phase. It can also have a finite magnitude in the SmA\* phase and diverge in the vicinity of the critical temperature due to the EC effect. The temperature-induced changes of this angle lead to the changes in $`\alpha `$ and $`\theta `$ angles which are determined by the orientation of the main axis. We show below that the temperature and electric field dependences of the main axis orientation can be deduced from linear and nonlinear optical experiments, in order to describe the ferroelectric switching, phase transitions and the EC effect in FLC cells.

## III Experiment

The principle arrangement of the experimental setup is described elsewhere M15 . Briefly, an OPO laser system is used as a source of the fundamental radiation, with the output wavelength of 537 nm, repetition rate 10 Hz and pulse duration 4 ns. An appropriate set of color filters is used for linear and nonlinear effects measurements. A photodiode or a PMT tube and gated electronics are used as a registration system for linear or nonlinear experiments. The temperature is varied from 20 to 55 $`{}_{}{}^{o}C`$, and automatically controlled by a digital thermocouple thermometer with the accuracy of $`1^oC`$. An electric field in the range of -15 to +15 MV/m is applied along $`OZ`$ axis through the ITO electrodes. Commercial FLC cells (E.H.C. Co., Tokyo) of a nominal thickness of 2 $`\mu m`$ are used. These cells have an unidirectional rubbed polyimide surface layer. The samples are prepared by a slow cooling of the mixture M15 , capillary filled in the isotropic phase. The critical temperature of the SmC\*-SmA\* phase transition for the studied mixture is about $`42^oC`$.

To determine the in-plane components of the nonlinear susceptibility, the SHG intensity vs. azimuthal angle is measured at normal incidence in transmission geometry for different polarization combinations - pp, ss, ps, sp, where the first letter denotes polarization of the incident light, and the second one - polarization of the light transmitted through the analyzer, placed after the sample. The observed anisotropy is shown in fig. 2a, b and correlates well with the dependence (4), obtained for a $`C_2`$ symmetry structure. No electric field is applied.

Ferroelectric switching of the cells is studied by measuring the anisotropy dependences of the SHG intensity at different values of the external DC electric field applied along the normal to the cell plane. Fig. 3a shows the SHG intensity in transmission geometry for pp-polarizations and for two applied DC fields of opposite sign. The dependences are shifted relative to each other, the phase difference is about 30o. In fig. 3b the dependences of the SHG intensities vs. electric field are shown for different anisotropy positions of the cell. The anisotropy position determines the contrast value of the switching dependences and the ratio between SHG responses corresponding to opposite biases. Switching in reflection geometry also exhibits a shift of the anisotropy dependences for opposite electric fields with a phase difference of about 20o.

Linear transmittance anisotropy is measured in crossed polarizers geometry and is shown in fig. 2c. Application of the electric field leads to a shift of the linear transmittance anisotropy, as in the nonlinear case. In both cases, we observed a shift of about 30o.

Temperature dependences of the SHG response are taken in reflection geometry at incident angle about 45o and in transmission at normal incidence. The cases of the presence and absence of external DC field are studied. Fig. 4a shows temperature dependences of the SHG intensity in pp transmission geometry for biased (filled circles and squares) and unbiased (open circles) cells. In the absence of the electric field the temperature dependence of the SHG intensity approaches a constant value at critical temperature in accordance with a power law at $`T<T_c`$, and shows no manifestation of the surface EC effect at $`T>T_c`$ in the SmA\* phase. For biased LC samples, a critical behavior of the SHG intensity due to the EC effect is obtained in the SmA\* phase in the vicinity of the phase transition, which is seen from the divergent behavior of the SHG intensity at $`T>T_c`$, similar to the critical behavior of the tilt angle.

## IV Discussion

FLC cell structure usually contains layers with different direction of spontaneous polarization - twisted or helicoidally wounded layers. This modulation of the space orientation of molecular dipoles is governed by the intermolecular forces and substrate influence. We suppose that in the absence of the DC electric field the cell has a net dipole moment and axis of $`C_2`$ symmetry laying in the cell plane. Then it is possible to approximate anisotropy dependences of SHG intensity in all possible polarization combinations in transmission at normal incidence by formulas analogous to (4) for pp case. In fig. 2 calculated curves are shown for all geometries. The approximation gives the ratios between in-plane components of nonlinear susceptibility:

$$\chi _{xxy}:\chi _{yyy}:\chi _{yxx}=1:23.9:3.3$$
(5)
The linear transmittance anisotropy has four equal peaks and no dark extinction, which is typical for twisted structures. Assuming the main optical axis to lay in the cell plane, approximation of the linear anisotropy dependence by (1) gives a birefringence of $`0.041`$, which corresponds to the characteristic values for smectic LC measured by other authors shgeo . Comparison of nonlinear and linear anisotropy dependences indicates that $`C_2`$ symmetry axis, corresponding to the maximum of the SHG anisotropy in pp-geometry, and the main optical axis, corresponding to the minimum of the linear transmittance anisotropy in sp-geometry, are parallel.

Temperature dependences of the SHG intensity for the applied electric field of different polarities reveal a hysteresis free SmC\*-SmA\* second order phase transition. A typical critical behavior in the vicinity of the SmC\*-SmA\* phase transition is obtained. A non-zero SHG signal in the SmA\* phase denotes that some part of the cell is still oriented and have a polar order with an in-plane component of the nonlinear polarization. The temperature dependence of the SHG response for reflection geometry shows qualitatively the same behavior as in transmission geometry, which is connected with the SmC\*-SmA\* phase transition in the sub-surface layer of the cell. The SHG temperature dependence can be explained by the interference of the field and temperature dependent and independent contributions to nonlinear polarization, responsible for the SHG signal. Then, we can write for the nonlinear polarization

$$\stackrel{}{𝐏}_{NL}^{2\omega }=\stackrel{}{𝐏}_{\text{surf}}^{2\omega }(E,T)+\stackrel{}{𝐏}_{\text{bulk}}^{2\omega }(E,T)+\stackrel{}{𝐏}_{\text{const}}^{2\omega },$$
(6)
where $`\stackrel{}{𝐏}_{\text{surf}}^{2\omega }(E,T)`$ and $`\stackrel{}{𝐏}_{\text{bulk}}^{2\omega }(E,T)`$ are electric field $`E`$ and temperature $`T`$ dependent contributions to the nonlinear polarization at double frequency from the sub-surface layers and the bulk, respectively, and $`\stackrel{}{𝐏}_{\text{const}}^{2\omega }`$ is a field and temperature independent component. For the transmitted SHG the correlation length is about 5 microns, while for reflection it is about 0.1 micron M12 , which indicates that in reflection a thin sub-surface layer participates in SHG. As SHG temperature dependences of sub-surface layers and the bulk have the same qualitative character, then either $`\stackrel{}{P}_{\text{bulk}}^{2\omega }(E,T)`$ and $`\stackrel{}{P}_{\text{surf}}^{2\omega }(E,T)`$ from (6) should be the same functions of $`E`$ and $`T`$, or the bulk contribution to nonlinear polarization should vanish.

Explanation of temperature dependences demands the existence of an electric field and temperature independent in-plane contribution $`\stackrel{}{P}_{\text{const}}^{2\omega }`$ to the nonlinear polarization. Orientation of the molecules in the bulk layers of the cell can be electric-field independent in some experimental geometries depending on the layer packing, but temperature dependence is obligatory because of the presence of the SmC\*-SmA\* phase transition. Then the non-switching layer should be in the sub-surface region. We can assume that a ”frozen” sub-surface layer, with a thickness smaller than the correlation length for SHG in reflection geometry, exists. It is strongly stabilized by the surface coupling, and does not respond to any external electric field or temperature variations, as well as the layer closest to the substrate. The latter has only z-component of the nonlinear polarization, as the dipole moments of its molecules are directed along the normal to the substrate plane, and then it does not contribute to the SHG at normal incidence. But the rest of the ”frozen” region contains layers with in-plane components of spontaneous polarization due to the twist near the surface, and thus participates in the SHG at any angle of incidence.

Let us figure out the interconnection between the temperature dependences of SHG and the order parameter of the SmC\*-SmA\* phase transition. Usually the order parameter of this phase transition is the tilt angle, as spontaneous polarization is compensated because of the helical structure. The thickness of the cells studied in this paper is smaller than the helix pitch, besides, anisotropic interaction with the substrate produces an additional orientational order of the FLC molecules, so that the cells have spontaneous polarization $`\stackrel{}{P}_{\text{sp}}`$ in SmC\* phase. $`\stackrel{}{P}_{\text{sp}}`$ decreases in SmA\* phase and is a linear function of the tilt angle $`\theta (E,T)`$. So we can write for the SHG intensity:

$$\sqrt{I^{2\omega }}\stackrel{}{P}_{NL}^{2\omega }\stackrel{}{P}_{\text{sp}}\theta (E,T),$$
(7)
assuming that $`\sqrt{I^{2\omega }(E,T)}`$ has the same critical dependence on temperature as the tilt angle. Then, according to (6) and (7), the temperature dependent contribution $`\sqrt{I^{2\omega }(E,T)}`$ in the absence of the electric field and for $`T<T_c`$ can be approximated by the following expression:

$$\sqrt{I^{2\omega }(E=0)}P_{const}^{2\omega }P_0^{2\omega }(1\frac{T}{T_c})^\beta ,$$
(8)
where $`P_{const}^{2\omega }`$ and $`P_0^{2\omega }`$ are proportional to isotropic and anisotropic contributions of nonlinear polarizations, $`T_c`$ is the temperature of the phase transition, and $`\beta `$ is the critical exponent. Approximation of the curve is shown in fig. 4a by a solid line, with $`\beta =0.31`$, in correspondence to the theory of phase transition in a superfluid helium.

To compare the temperature behavior of the nonlinear polarization in the vicinity of the phase transition for reflection and transmission cases, we introduce the contrast of the dependences as

$$K=\frac{\sqrt{I^{2\omega }(+E)}\sqrt{I^{2\omega }(E)}}{\sqrt{I^{2\omega }(+E)}+\sqrt{I^{2\omega }(E)}},$$
(9)
where $`I^{2\omega }(\pm E)`$ are the SHG intensities for positive and negative voltages applied to the cell. In fig. 4b the contrast vs. temperature is shown for reflection and transmission geometries. The contrast determined in such a way shows a behavior similar to that of the order parameter near the phase transition. Introducing SHG intensity dependences on temperature in the presence of the electric field at $`T<T_c`$ in the way, analogous to (8),

$$\sqrt{I^{2\omega }(\pm E)}P_{const}^{2\omega }\pm P_0^{2\omega }(\pm E)(1\frac{T}{T_c})^\beta ,$$
(10)
and supposing $`P_0^{2\omega }(E)=P_0^{2\omega }(E)`$, we can express the contrast K as:

$$K\frac{P_0^{2\omega }(E)}{P_{const}^{2\omega }}(1\frac{T}{T_c})^\beta .$$
(11)
Approximation of experimental dependences by (11) in the SmC\* phase gives $`\beta =0.31`$ (fig. 4b), as in the case of the absence of the field. The amplitude of the contrast (11) depends on the ratio $`P_0^{2\omega }(E)/P_{const}^{2\omega }`$ of field dependent and independent contributions to the nonlinear polarization. In transmission this ratio can be bigger due to the larger amount of layers participating in ferroelectric switching, besides, the field independent contribution of the boundary layer leads to the contrast decrease in reflection geometry. These considerations explain why the contrast of temperature dependences and electric field-induced shift of anisotropy dependences in transmission are bigger than in reflection geometry, as obtained in the experiment.

Temperature dependence of the contrast in the region of the EC effect, at $`T>T_c`$, can be written in the same way, considering the temperature dependence $`\sqrt{I^{2\omega }(E,T)}`$ similar to the temperature dependence of the tilt angle above the critical temperature:

$$\sqrt{I^{2\omega }(\pm E)}P_{const}^{2\omega }\pm P_0^{2\omega }(\pm E)(1\frac{T}{T_c})^\gamma ,$$
(12)
where $`\gamma `$ is the critical exponent, which gives under the same considerations the contrast for $`T>T_c`$:

$$K\frac{P_0^{2\omega }(E)}{P_{const}^{2\omega }}(1\frac{T}{T_c})^\gamma ,$$
(13)
Approximation of experimental data (fig. 4b) gives $`\gamma =1.4`$. Thus we observe that both critical exponents are in agreement with the superfluid helium model.

Application of a DC electric field to the FLC leads to a shift of anisotropy dependences of the linear and nonlinear responses, as shown in fig. 3a. Similar pictures are obtained in the reflection case. As the character of anisotropy dependences is the same in the presence of the field, our suggestion of the switching mechanism based on the symmetry axis rotation in the cell plane seems to be rather realistic. The shift originates from the rotation of the molecules in the switching layers, the resulting nonlinear polarization of the cell is the sum of unchanged quadratic polarization from the ”frozen” layer and from field-dependent region with the new position of the symmetry axes. Thus for approximating the SHG anisotropy dependence in the presence of the electric field, we divide the nonlinear polarization into two parts, and (4) can be rewritten as:

$$I_{pp}^{2\omega }(gP_{NL,pp}^{2\omega }(\theta +d\theta )+(1g)P_{NL,pp}^{2\omega }(\theta ))^2,$$
(14)
where $`g`$ is the effective thickness of switching layer relatively to the non-switching one and $`d\theta `$ is a field-induced change of the azimuthal angle. The value of $`g`$ is estimated to be 0.5 which indicates that the magnitude of the field-independent contribution is comparable to the field-dependent one, allowing effective mutual interference in accordance with (6). Symmetry axis shift $`d\theta `$ between opposite values of the voltage is about 52o. For electric field dependences of the SHG intensity, shown in fig. 3b, the independent variable is the electric field and in (14) the azimuthal angle $`\theta `$ is fixed by the anisotropy position. Then $`d\theta `$ becomes a function of the electric field, chosen, according to (3), as a saturation function:

$$d\theta (E)=b\mathrm{arctan}(c(EdE)),$$
(15)
where $`b`$ and $`c`$ are constants, depending on the strength of interaction between field and angle, and $`dE`$ is a constant determined by the switching history, as ferroelectric switching of SHG intensity exhibits characteristic field hysteresis. Experimental dependences in fig. 3b are taken for two anisotropy positions with the difference in $`\theta 130`$ degrees. Fitting one of the curves by (14) with (15) allows to extract $`b=0.5`$ and $`c=0.63`$. Varying then only $`dE`$ and $`\theta `$ during approximation of the second curve gives the difference of the anisotropy position of about 126 degrees, in good agreement with experimental conditions. Values of $`dE`$ in both cases lay in the region of typical hysteresis width and do not exceed 1 V/micron. The value of $`b`$ corresponds to the tilt of the symmetry axis in the cell plane, discussed above (3).

In conclusion, the electroclinic effect, the SmC\*-SmA\* phase transition and ferroelectric switching have been studied in thin planar cells of ferroelectric chiral liquid crystals by means of electrooptic and second harmonic generation techniques. The analysis of the temperature dependences of nonlinear quadratic response in the critical region leads to the assumption of a strong surface coupling existence, resulting in the stabilizing of several ”frozen” subsurface twisted layers, independent on electric field and temperature. Critical exponents in the vicinity of the SmC\*-SmA\* phase transition show a behavior corresponding to the superfluid helium theory. Ferroelectric switching has been observed in linear and nonlinear anisotropic responses and has been explained by the model, assuming in-plane rotation of symmetry and main optical axes of the cells in the presence of external electric field.

###### Acknowledgements.

We thank Prof. G. Heppke’s group for synthesis of the samples used in our work. The work was supported by INTAS grant 2002-113/F1b, President grant ”Leading Russian Scientific Schools” 1604.2003.2 and EU network grant ERB FMRXCT 980209 (SILC). S.S acknowledges funding from European Regional Development Fund and from the Ministerio de Ciencia y Tecnologia through the Ramon y Cajal program.


text:
# Bose-Einstein Condensation in Random Directed Networks

## I Introduction

Recently there has been much interest in random growing networks review ; review1 ; review2 , both from the point of view of theoretical modelling, as well as the empirical study of real networks. There is considerable evidence that while traditional Erdős-Renyi random graphs have a Poisson degree distribution bela , real graphs and random growing graphs have a power-law degree distribution.

Of particular interest are directed networks, which can be used to model systems in which directed flow takes place. Such networks include the world wide web 1 , the phone-call graph 2  and the citation graph 3 .

As a complement to computer simulations and exact solutions of simplified systems, thermodynamic formulations have been used to study a number of complex co-operative systems including granular media, econophysics, breaking phenomena and many others. In a recent paper bb1 , Bianconi and Barabási (BB) mapped the different behaviour of random growing networks with fitness onto the thermodynamically distinct phases of a free Bose gas. The fitness model predicts that, in the large network limit, the fittest node will have the most links. This is called the fit-get-rich (FGR) phase. Unlike the scale-free model (SFM), in which the degree distribution of the network is power-law, the FGR behaviour (phase in the thermodynamic language) has nodes with a very large degree which dominate the network. Another phase revealed in bb1  is the Bose-Einstein condensation (BEC), when the fittest node acquires a finite fraction of the total degree. In contrast to BEC, in the FGR phase the richest node is not an absolute winner, since its share of the links (i.e. the ratio of its degree and the total degree of the network) decays to zero for large system sizes, whereas in BEC the winner maintains its share irrespective of system size.

In this case, the fact that this ratio is constant corresponds to the extensivity property of a Bose gas, when the gas keeps a finite fraction of its particles in the ground state. In this paper we examine the phenomenon of Bose-Einstein condensation, which was considered for an undirected graph in bb ; bb1  on the random directed growing network introduced in er .

To do this we work with the model introduced in er , introducing a dependence between the vertex fitness and an energy. As each vertex has two fitnesses (for in-degree and out-degree) then it is necessary to assign two different energy levels to represent each vertex. This is detailed in the first sub-section, together with the description of edges in terms of particles assigned to two energy levels.

The necessary coexistence of two sub-systems forming the network and the formulation of an equilibrium condition in a formal way is discussed in detail in sections 3 and 4.

The occurrence of Bose-Einstein condensation, and its interpretation in a physical framework makes it possible to describe the directed network in terms of canonical concepts of statistical physics, such as a phase diagram and first and second order phase transitions. These are introduced in the following sections.

In the conclusions we emphasize the usefulness of this formulation in providing a clear and concise interpretation of different phenomena in networks. So, the possibility of coexistence of phases in which the in-degree distribution exhibits a clear winner whereas the out-degree shows scale free behaviour, or the simultaneous existence of two different winners (bipolarity) and the conditions for its existence are predicted in thermodynamic language.

## II The Model

In BB bb1 , a correspondence between fitness $`\eta `$ and energy $`ϵ`$ given by

$$\eta =e^{\beta \epsilon }$$
(1)
was introduced. In our model the vertices have, in general, different fitnesses $`a`$ and $`b`$ for the addition of an in- or out-degree. As in er , here the model consists of the addition of bare vertices (i.e. without edges, but with fitness $`a`$ for in-degree and $`b`$ for out-degree) to the network with probability $`p`$, and the creation of directed edges between vertices with probability $`q=1p`$.

A kinetic equation describing the process of directed networks must include the kinetics of both in-degree and out-degree simultaneously since they coexist and influence each other. Then, the kinetic equation for the mean number of vertices with in-degree $`i`$ and out-degree $`j`$ is, as in er ,

$$\frac{N_{ij}(a,b)}{t}=\frac{qa}{M_1}[iN_{i1,j}(i+1)N_{ij}]+\frac{qb}{M_2}[jN_{i,j1}(j+1)N_{ij}]+p\delta _{i0}\delta _{j0}f(a,b).$$
(2)
where the fitnesses $`(a,b)`$ are chosen from the fitness distribution $`f(a,b)`$. The normalisation factors $`M_1`$ and $`M_2`$ are given by

$$M_1=\underset{ijab}{}(i+1)aN_{ij}(a,b)$$
(3)
and

$$M_2=\underset{ijab}{}(j+1)bN_{ij}(a,b).$$
(4)
The first term in the first square brackets of Eq.(2) represents the increase in $`N_{ij}`$ when vertices with in-degree $`i1`$ and out-degree $`j`$ gain an incoming edge and the second term represents the corresponding loss. The second square brackets contain the analogous terms for outgoing edges and the last term ensures the continuous addition of new vertices with fitnesses $`a,b`$.

The translation of this problem into an energetic formulation is straightforward. We associate the addition of a vertex with fitnesses $`(a,b)`$ to the network with the creation of two energy levels; one representing the fitness of the incoming edge and the other for the outgoing edge. This means that the edge can be mapped into two separate isolated sub-systems. Then the creation of a directed edge corresponds to the creation of two particles, one in each sub-system, simultaneously. The particles are created in the energy levels corresponding to incoming fitness of the receiving vertex and the outgoing fitness of the emitter vertex, as Fig. 1 illustrates.

Then, using the energy variables, rather that the fitness, we have

$$a=e^{\beta _1\epsilon _1},$$
(5)
$$b=e^{\beta _2\epsilon _2}$$
(6)
and

$$\frac{N_{ij}(\epsilon _1,\epsilon _2)}{t}=\frac{qe^{\beta _1\epsilon _1}}{M_1}[iN_{i1,j}(i+1)N_{ij}]+\frac{qe^{\beta _2\epsilon _2}}{M_2}[jN_{i,j1}(j+1)N_{ij}]+p\delta _{i0}\delta _{j0}f(\epsilon _1,\epsilon _2).$$
(7)
and the normalisations are re-written as

$$M_1=\underset{ijab}{}(i+1)e^{\beta _1\epsilon _1}N_{ij}(\epsilon _1,\epsilon _2)$$
(8)
and

$$M_2=\underset{ijab}{}(j+1)e^{\beta _2\epsilon _2}N_{ij}(\epsilon _1,\epsilon _2).$$
(9)
As we distinguish between incoming and outgoing edges, and their respective fitnesses, when going to the energy representation the edges and their fitnesses must also be distinguished as belonging to different isolated sub-systems of the same system. Let us denote the sub-systems as 1 and 2 So, the creation of a vertex in the network implies the simultaneous creation of an energy level in each sub-system, i.e. $`\epsilon _{1\text{ }}`$and $`\epsilon _{2\text{ }}`$.

The creation of an edge joining two vertices implies the simultaneous creation of one particle in each sub-system in the energy levels corresponding to the fitnesses of the vertices gaining an in- and out-degree, as Fig.1 illustrates.

It should be noted that we are doing a simultaneous description of two isolated sub-systems which compose the whole system (the network). These sub-systems do not exchange energy or particles, but are correlated in the sense that the creation of an energy level (a particle) in one of the sub-systems implies the simultaneous phenomenon in the other.

By definition of the model, the total number of particles increases linearly with time, so

$$\underset{i,j,\epsilon _1,\epsilon _2}{}N_{ij}(\epsilon _1,\epsilon _2)=pt.$$
(10)
Let us define the reduced moments $`m_1`$ and $`m_2`$ by

$$M_1(t)=m_1t,$$
(11)
$$M_2(t)=m_2t$$
(12)
and introduce $`n_{ij}`$ as

$$N_{ij}(t)=n_{ij}t.$$
(13)
As shall be seen these magnitudes will be useful to calculate the characteristics of the network. Then, following the procedure in er , we can show that

$$m_1=\frac{q}{m_1}\underset{i,j,\epsilon _1,\epsilon _2}{}e^{2\beta _1\epsilon _1}(i+1)n_{ij}+pe^{\beta _1\epsilon _1}f(\epsilon _1,\epsilon _2)𝑑\epsilon _1𝑑\epsilon _2$$
(14)
and

$$m_2=\frac{q}{m_2}\underset{i,j,\epsilon _1,\epsilon _2}{}e^{2\beta _2\epsilon _2}(j+1)n_{ij}+pe^{\beta _2\epsilon _2}f(\epsilon _1,\epsilon _2)𝑑\epsilon _1𝑑\epsilon _2$$
(15)
which determines $`m_1`$ and $`m_2`$ as functions of $`p`$, $`\beta _1`$ and $`\beta _2`$.

### II.1 The coexistence of two systems

In this representation, as we already pointed out, the addition of a vertex means the creation of two energy levels, each corresponding to its fitness. When an edge joining two vertices is created, we add one particle in each corresponding level.

The distribution of in- and out-degree then corresponds to the distribution of particles in systems 1 and 2

$$g_i=\underset{j}{}n_{ij}=p\frac{\mathrm{\Gamma }(i+1)\mathrm{\Gamma }(1+e^{\beta (\epsilon _1\mu _1)})}{\mathrm{\Gamma }(i+2+e^{\beta (\epsilon _1\mu _1)})}e^{\beta (\epsilon _1\mu _1)}f(\epsilon _1,\epsilon _2)$$
(16)
is the number of vertices with in-degree $`i`$, and

$$h_j=\underset{i}{}n_{ij}=p\frac{\mathrm{\Gamma }(j+1)\mathrm{\Gamma }(1+e^{\beta (\epsilon _2\mu _2)})}{\mathrm{\Gamma }(j+2+e^{\beta (\epsilon _2\mu _2)})}e^{\beta (\epsilon _2\mu _2)}f(\epsilon _1,\epsilon _2)$$
(17)
is the number of vertices with out-degree $`j`$. For large degree their asymptotic expressions are

$$g_ii^{(1+e^{\beta (\epsilon _1\mu _1)})}$$
(18)
and

$$h_jj^{(1+e^{\beta (\epsilon _2\mu _2)})}.$$
(19)
where we have introduced

$$\mu _k=\frac{m_k}{q}$$
(20)
for $`k=1`$ and $`k=2`$ as chemical potentials. Unlike BB, the chemical potential is introduced here as an exact expression, not as an asymptotic one.

### II.2 “Equilibrium” condition

The equations for $`m_1`$ and $`m_2`$ can be transformed by introducing the generating function

$$g(x,y,\epsilon _1,\epsilon _2)=\underset{i,j}{}x^iy^jn_{ij}(\epsilon _1,\epsilon _2)$$
(21)
where

$$g(1,1,\epsilon _1,\epsilon _2)=pf(\epsilon _1,\epsilon _2)$$
(22)
with

$$\frac{g}{x}(1,1,\epsilon _1,\epsilon _2)=\frac{pe^{\beta _1\epsilon _1}}{e^{\beta _1\mu _1}[1e^{\beta _1(\epsilon _1\mu _1)}]}$$
(23)
and

$$\frac{g}{y}(1,1,\epsilon _1,\epsilon _2)=\frac{pe^{\beta _2\epsilon _2}}{e^{\beta _2\mu _2}[1e^{\beta _2(\epsilon _2\mu _2)}]}.$$
(24)
Then the equations for $`m_k`$ $`(k=1,2)`$ become

$$\frac{p}{q}\frac{d\epsilon _1d\epsilon _2f(\epsilon _1,\epsilon _2)}{e^{\beta _k(\epsilon _k\mu _k)}1}=1.$$
(25)
This means that

$$\frac{d\epsilon _1d\epsilon _2f(\epsilon _1,\epsilon _2)}{e^{\beta _1(\epsilon _1\mu _1)}1}=\frac{d\epsilon _1d\epsilon _2f(\epsilon _1,\epsilon _2)}{e^{\beta _2(\epsilon _2\mu _2)}1}.$$
(26)
Unlike BB, this is a new condition generated by the model. Let us denote it as a generalized equilibrium condition. A *strong equilibrium* condition implies that

$$\epsilon _1=\epsilon _2$$
(27)
and

$$\mu _1=\mu _2.$$
(28)
This condition will be discussed later. In this formalism the occupation probability is, then,

$$n_k(\epsilon _k,\mu _k)=\frac{p}{q}\frac{1}{e^{\beta _k(\epsilon _k\mu _k)}1},$$
(29)
for $`k=1,2`$, representing the probability for a given edge to belong to a given vertex (probability of a particle to belong to an energy level). If, now, the fitness distribution is separable then

$$f(\epsilon _1,\epsilon _2)=f_1(\epsilon _1)f_2(\epsilon _2)$$
(30)
and

$$\frac{p}{q}\frac{d\epsilon _kf(\epsilon _k)}{e^{\beta _k(\epsilon _k\mu _k)}1}=1$$
(31)
separately for each system $`k=1,2`$. In each case, $`\mu _k`$ is the solution of this equation, which depends on the density of states $`f_k(\epsilon )`$.

### II.3 Bose-Einstein condensation

To illustrate, let us start with a particular case, when $`f(\epsilon )=C\epsilon ^\theta `$ with $`0<\epsilon <\epsilon _{\mathrm{max}}`$. The normalisation $`C=\frac{\theta +1}{\epsilon _{\mathrm{max}}^{\theta +1}}`$ as in BB.

Then the integral in Eq.(25),

$$I(\beta ,\mu )=\frac{p(\theta +1)}{q\epsilon _{\mathrm{max}}^{\theta +1}}\frac{dϵϵ^\theta }{e^{\beta (\epsilon \mu )}1}=1.$$
(32)
But this integral reaches its maximum value when $`\mu =0`$. Then if for a given value of $`\beta `$ we have

$$\frac{d\epsilon f(\epsilon )}{e^{\beta \epsilon }1}<1,$$
(33)
then we have a Bose-Einstein Condensation (BEC). If we calculate the temperature for BEC in this particular case, we obtain

$$T_{BE}=\frac{1}{\epsilon _{\mathrm{max}}}\left\{\frac{p}{q}\mathrm{\Gamma }(\theta +2)\zeta (\theta +1)\right\}^{\frac{1}{\theta +1}}.$$
(34)
Unlike BB,the factor $`\frac{p}{q}`$ appears here. If $`p=q`$ we recover BB result. Even in this simple example, $`T_{BE}`$ may be different for the sub-systems 1 and 2, if the value of $`\theta `$ is different. In a more general form, we can state the condition for Bose Einstein condensation. Then, as

$$I(\beta _k,\mu _k)=\frac{p}{q}\frac{d\epsilon _1d\epsilon _2f(\epsilon _1,\epsilon _2)}{e^{\beta _k(\epsilon _k\mu _k)}1}=1$$
(35)
for $`k=1,2`$., if $`I(\beta _k,\mu _k)<1`$ then BEC occurs. The conditions for BEC depend on $`f(\epsilon _1,\epsilon _2).`$

### II.4 Phase diagram

Below $`T_{BE}`$ there is BEC. Above it the phase is fit-get-richer as in bb1 , but here both phases must be considered for both systems. Then, a phase diagram can be drawn with the temperatures $`T_1`$ =$`\frac{1}{\beta _1}`$ and $`T_2`$ =$`\frac{1}{\beta _2}`$ characterizing each sub-system (Fig. 2).

Imposing a given dependence of fitness with energy so that

$$\frac{d\epsilon _1d\epsilon _2f(\epsilon _1,\epsilon _2)}{e^{\beta _1(\epsilon _1\mu _1)}1}\frac{d\epsilon _1d\epsilon _2f(\epsilon _1,\epsilon _2)}{e^{\beta _2(\epsilon _2\mu _2)}1}$$
(36)
means that the equilibrium is broken, (i.e. in a computer simulation we would impose temperatures that do not satisfy Eq. 26 ). Then there can be a winner in one of the phases. One of the vertices accumulates, say, most of the in-degrees, coming from the other vertices. If both phases are in BEC there is a vertex with most of the incoming and another with most of the outgoing edges. Then both vertices will be highly correlated.

### II.5 Out of equilibrium first order phase transition

Though a more complicated situations are possible, let us illustrate the case $`T_{1BE}=T_{2BE}`$. The point A represents that both sub-systems are in FGR phase. If, for example, we make the simulations with the parameters given by A1 then will be an FGR and a BEC phase coexisting, but not in equilibrium since one of them satisfies Eq. 35 but not the other. There is a coexistence of phases out of equilibrium. We might call this an out of equilibrium first order phase transition. The motion along AA1A2…leads to the region in which BEC occurs for both phases. Here is hard to speak about equilibrium unless A belongs to the line $`T_1=T_2`$. All we know in BEC region is that none of the sub-systems satisfies Eq. 35

The concept of equilibrium in this model is purely formal since the sub-systems are not in contact, but it is a comfortable tool to make a map of directed networks into Bose-Einstein statistics.

### II.6 The scale free model

If one of the systems is such that $`f(\epsilon )=\delta (\epsilon )`$, then Eq. (31) gives

$$\frac{p}{q}\frac{1}{e^{\beta \mu }1}=1$$
(37)
then

$$g_ii^{(1+e^{\beta \mu })}i^{(2+\frac{p}{q})}$$
(38)
with a corresponding expression for $`h_j`$. If $`p=q`$ the scale free model is recovered with a degree distribution $`P_kk^3`$. Thus the model with directed edges gives the possibility of a scale free phase for both in the in-degree and the out-degree.

### II.7 Second order phase transition

If in Fig. 3 we move in such a way as to cross into the BEC region through the point B then BEC occurs simultaneously for both phases. No coexistence of phases will occur, the system goes from FGR to BEC or vice-versa instantly. This is a second order phase transition since the system as a whole changes its phase. Then the point B is the critical point for the transition.

BEC is itself a second order phase transition since it effects the system as a whole, so that when sub-system 1 or 2 experiences BEC, the sub-system itself suffers a second order phase transition, but our system contains both sub-systems 1 and 2 so that when they are in different phases we consider our whole system as composed of two phases.

## III Discussion

We have studied Bose-Einstein condensation in a random growing network model. In this model, with a formulation considering incoming and outgoing edges, a formal analogy can be adopted to describe the change in the behaviour of the directed edges in terms of Bose statistics.

In this formulation, as in that of BB, the temperature plays the role of a dummy variable. Real systems are characterised by the functional form of the distribution function of fitness and there is no $`\beta `$ to consider. It only emerges when we translate from the language of fitness to that of energies. But the formulation of networks in terms of quantum statistics and the introduction of temperature is an elegant and simple description of the behaviour of networks under different conditions. In this sense, the temperature plays its role in the simulations of the networks, determining the strength of the dependence between the fitness and its associated energy, and this dependence is reflected in the behaviour of the network by determining its position in the phase diagram.

Though phase diagrams and phase transitions are not exclusive to thermodynamic formulations (see, e.g. Dorogovtsev and Mendes dm ), such a formulation is a very powerful tool to describe different phases of random directed networks, where the presence of directed edges make the description in terms of fitness difficult. Considering this network as a thermodynamic system, makes it simple to understand, provided we appropriately interpret the thermodynamic parameters.

## IV Acknowledgements

We would like to thank The Royal Society, London and the London Mathematical Society for their financial support. One of us (O. S) is greatly indebted to the staff of the Department of Mathematical Sciences at Brunel University for their warm hospitality.

## V Figure Captions

1. Correspondence between fitness and energy for directed networks.

2. Phase diagram.

3. First and second order phase transitions.


text:
# Chapter 1 A New Estimate of the Extragalactic Gamma-Ray Background from EGRET Data

## Chapter 1 A New Estimate of the Extragalactic Gamma-Ray Background from EGRET Data

### Abstract

We use the GALPROP model for cosmic-ray (CR) propagation to obtain our best prediction of the Galactic component of gamma rays, and show that away from the Galactic plane it gives an accurate prediction of the observed EGRET intensities. On this basis we re-evaluate the isotropic extragalactic gamma-ray background (EGRB). We find that for some energies previous work underestimated the Galactic contribution and hence overestimated the background. Our new EGRB spectrum shows a positive curvature similar to that expected for models of the extragalactic emission based on the blazar population.

### 1. Introduction

The GALPROP model for CR propagation produces explicit predictions for the angular distribution of diffuse Galactic $`\gamma `$-rays. In a companion paper  we describe a comparison of a model for Galactic diffuse continuum $`\gamma `$-rays with EGRET data. In this paper we use this model to determine the EGRB.

### 2. Data and method

Given the success of the model  in reproducing the $`\gamma `$-ray sky, we can use it to make improved estimates of the EGRB. Since the model is nevertheless not exact the best approach is to fit the observed intensities with a free scaling factor; in this way the EGRB is determined as the intercept, thus removing any residual uncertainty in the absolute level of the Galactic components. (This is the same method as in , the difference lies in the model). To reduce the effects of Galactic structure, point sources etc. the fits are made excluding the plane, using the region $`360^{}<l<0^{},10^{}<|b|<90^{}`$ (“region G”); ideally both IC and gas-related components would be left free but they are difficult to separate statistically at high latitudes, so we make a linear fit to the total IC+$`\pi ^0`$-decay+bremsstrahlung, with the scaling factor and EGRB as parameters. The fit and errors are based on a simple $`\chi ^2`$ analysis, with ($`l,b`$) bins $`360^{}\times 2^{}`$ to obtain sufficient statistics (at least 10 counts per bin were required). For comparison we also made fits to the entire sky; in this case IC and gas-related contributions are easily separated; the two fit regions then give some indication of the model-dependent systematic error in our EGRB estimates.

The scaling factors determined for the region G fits reflect the deviations from the model and are shown in Table 1. They are typically between 0.8 and 1.2 which is satisfactory. (Only statistical errors are shown in Table 1.) For 30–50 MeV, 1–2 and 2–4 GeV the scaling factors deviate further from 1 reflecting the discrepancy in the spectrum so that these are the least reliable ranges of our EGRB determination. The EGRB is however not very sensitive to the scaling factor.

Table 2 summarizes the fitted EGRB values. The two fitted regions (region G and all sky) give consistent results, indicating that there is no large systematic effect; it shows a model-dependent systematic uncertainty of 5–30%, comparable to the formal statistical errors. This is comparable to the $`15`$% systematic uncertainty on EGRET data so we adopt 30% for our total error estimate.

Fig. 1 shows the extragalactic X- and $`\gamma `$-ray background, using the compilation by  but with our new EGRET values, and also updated COMPTEL results ( and references therein). Our estimates lie significantly below those of , in each energy range. The positive curvature in our EGRB spectrum is interesting and not unexpected  but in view of the systematic uncertainties should not be taken too literally; a similar, less pronounced effect is present in the Sreekmar spectrum. A power-law is a poor fit to our spectrum. Although the 50 MeV – 2 GeV range can be fit satisfactorily by a power law, it is clearly inconsistent with the points above 2 GeV.

The reason for the difference between our spectrum and that given in  is the improved modelling of high-latitude $`\gamma `$-rays based on inverse Compton emission from the halo. An indication of this effect is already apparent in Table 1 of  which requires scaling factors up to 1.8 above 1 GeV. Detailed comparisons of our model with the EGRET data will be given in a forthcoming publication.

I.V.M. acknowledges partial support from a NASA Astrophysics Theory Program grant.

### 3. References

1. Salamon M.H., Stecker F.W. 1998, ApJ 493, 547

2. Sreekumar P. et. al. 1998, ApJ 494, 523

3. Strong A.W., Moskalenko I.V., Reimer O. 2000, ApJ 537, 763

4. Strong A.W., Moskalenko I.V., Reimer O. 2003, these Proceedings

5. Weidenspointner G. et al. 2000, in AIP Conference Proceedings 510, 467