text:
# Two-photon interference of multimode two-photon pairs with an unbalanced interferometer

## Abstract

Two-photon interference of multimode two-photon pairs produced by an optical parametric oscillator has been observed for the first time with an unbalanced interferometer. The time correlation between the multimode two photons has a multi-peaked structure. This property of the multimode two-photon state induces two-photon interference depending on delay time. The nonclassicality of this interference is also discussed.

Quantum interference is one of the most interesting phenomena in quantum physics. Since the observation of nonclassical effects in the interference of two photons by Ghosh and Mandel Ghosh , several types of quantum interference experiments have been demonstrated using correlated two-photon pairs generated by spontaneous parametric down-conversion (SPDC) Hong ; Ou1989 ; Ou1990 ; OuFranson ; Chiao ; Rarity ; Brendel ; ShihMach . In Refs. Ghosh ; Hong ; Ou1989 ; Ou1990 ; OuFranson ; Chiao ; Rarity ; Brendel ; ShihMach , higher visibility of two-photon interference than in the classical case is discussed as a typical nonclassical effect. Another feature of quantum interference is characterized by the shorter period of interference than in the classical case. Fonseca et al. have observed a twice narrower interference pattern than single-photon interference with correlated two-photon pairs generated by SPDC Fonseca . Quantum lithography has been proposed by Boto et al. in order to surpass classical diffraction limit utilizing this feature of quantum interference Boto . D’Angelo et al. have reported a proof-of-principle quantum lithography ShihLithography . This feature of quantum interference has also been confirmed using a conventional Mach-Zehnder interferometer Edamatsu .

In these quantum interference experiments, SPDC process has been used to prepare correlated two photons. Ou et al. have succeeded in generating correlated two photons with new property by using an optical parametric oscillator (OPO) OuPRL ; OuPRA . The bandwidth of the two photons is narrow and the correlation time is long ($``$10 ns). This property of the two photons produced by an OPO has enabled to directly observe their correlation function by coincidence counting. These narrow-band two-photon state have been used to observe nonclassical photon statistics OuNonclassical . Goto et al. have recently reported the observation of another type of correlated two photons, that is, multimode two photons produced by an OPO Goto . The correlation function of the multimode two-photon pairs has a multi-peaked structure. In this Brief Report, we report the first observation of two-photon interference of the multimode two-photon pairs with an unbalanced interferometer. Our experiment is shown in Fig. 1. The output beam from an OPO is incident at one of the input ports of an unbalanced interferometer. The correlation function of one of the outputs of the interferometer is observed with two photodetectors and a coincidence counter. First, we discuss briefly what happens in this experiment. Next, we explain our experimental setup and show our experimental results. Finally, we discuss our results with theoretical calculation of the correlation function of the output from the interferometer. The nonclassical feature of the two-photon interference is also discussed.

The distance between multimode two photons produced by an OPO is $`n\tau _rc`$, where $`\tau _r`$ is the round-trip time of the OPO cavity, $`c`$ is the speed of light in the vacuum, and $`n`$ is an integer Goto . We assume that the propagation time difference, $`T`$, between the short and long paths in the interferometer is nearly equal to $`\tau _r/2`$. There are two cases, Case1 and Case2: in Case 1, both the photons are reflected or transmitted at the first beam splitter of the interferometer; in Case 2, one of the two photons is reflected and the other is transmitted there. The distance between two photons in the output of the interferometer is $`n\tau _rc`$ in Case 1 and $`(n+1/2)\tau _rc`$ in Case 2. Therefore, the two-photon pairs in Case 1 and Case 2 will induce the peaks of coincidence counts at delay times $`n\tau _r`$ and $`(n+1/2)\tau _r`$, respectively. This enables to distinguish Case 1 and Case 2 through delay time. The height of the peaks of coincidence counts at delay times $`(n+1/2)\tau _r`$ will be constant with regard to the path-length difference of the interferometer because two photons in Case 2 do not interfere with each other. On the other hand, in Case 1, the two-photon pairs provide two-photon interference because we can not say which path a two-photon pair propagate on. Therefore, the height of the peaks of coincidence counts at delay times $`n\tau _r`$ will change with regard to the path-length difference of the interferometer. Thus, it is expected that two-photon interference depending on delay time will be observed.

The schematic of the experimental setup is shown in Fig. 2. The differences between this setup and that used in our previous work Goto  are an unbalanced interferometer between the OPO and detectors, and a locking beam for the phase lock of the interferometer. The light source is a single-mode cw Ti:Sapphire laser of wavelength 860 nm. The round-trip length of the OPO is set long (560 mm) in order to time resolve the oscillatory structure in the correlation function. The output beam from the OPO is incident at one of the input ports of the interferometer. The path-length difference of the interferometer is set at about 29 cm, which gives $`T0.97`$ ns. The phase difference of the interferometer is locked by a servo-control system. To prevent the locking beam for the interferometer from making noise, the beam propagates in the interferometer in the opposite direction to the signal beam from the OPO. Furthermore, the polarization of the locking beam is perpendicular to that of the signal in order to remove the beam by a polarization beam splitter (see Fig. 2). One of the two outputs of the interferometer is split into two with a 50/50 beam splitter. The two beams are coupled to optical fibers and detected with avalanche photodiodes (APD, EG&G SPCM-AQR-14). The coincidence counts of the signals from the two APDs are measured with a time-to-amplitude converter (TAC, ORTEC 567) and a multichannel analyzer (MCA, NAIG E-562). We measured the coincidence counts at $`\theta =(\pi /8)\times j`$ ($`j=0,1,\mathrm{},8`$), where $`\theta `$ is the phase difference of the interferometer defined as the intensity of the output of the interferometer is proportional to $`(1+\mathrm{cos}\theta )`$ when classical light of wavelength 860 nm is incident to the interferometer. The experimental results at $`\theta =(\pi /8)\times j`$ ($`j=0,1,\mathrm{},8`$) are shown in Fig. 3(a)-(i).

From now, we discuss our experimental results. First, we derive the correlation function of the output from the interferometer. The output operator, $`\stackrel{~}{a}(t)`$, from one of the output ports of the interferometer is expressed as

$`\stackrel{~}{a}(t)=`$ $`{\displaystyle \frac{\stackrel{~}{a}_{out}(tT_S)+\stackrel{~}{c}_{in}(tT_S)}{2}}`$
$`{\displaystyle \frac{\stackrel{~}{c}_{in}(tT_L)\stackrel{~}{a}_{out}(tT_L)}{2}}.`$ (1)
Here $`\stackrel{~}{a}(t)`$ denotes a Fourier transform of a field operator, $`a(\mathrm{\Omega })`$, of frequency $`\omega _0+\mathrm{\Omega }`$ ($`\omega _0`$ is the degenerate frequancy of the OPO). That is, it is defined as

$`\stackrel{~}{a}(t)={\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle 𝑑\mathrm{\Omega }a(\mathrm{\Omega })e^{i(\omega _0+\mathrm{\Omega })t}}.`$ (2)
$`a_{out}`$ is the output operator of the OPO far below threshold Goto  and $`c_{in}`$ is an annihilation operator of the vacuum entering the interferometer from the other of the input ports (see Fig. 1). $`cT_S`$ and $`cT_L`$ are the short and long path lengths of the interferometer, respectively. The intensity correlation function is derived as follows OuPRA ; Goto :

$`\mathrm{\Gamma }(\tau )`$ $`=\stackrel{~}{a}^{}(t)\stackrel{~}{a}^{}(t+\tau )\stackrel{~}{a}(t+\tau )\stackrel{~}{a}(t)`$
$`=\left[{\displaystyle \frac{2\sqrt{\mathrm{\Gamma }_0(\tau )}\mathrm{cos}\theta +\sqrt{\mathrm{\Gamma }_0(\tau T)}+\sqrt{\mathrm{\Gamma }_0(\tau +T)}}{4}}\right]^2+\left|{\displaystyle \frac{2\sqrt{\delta \mathrm{\Gamma }_0(0)}+e^{i\theta }\sqrt{\delta \mathrm{\Gamma }_0(T)}+e^{i\theta }\sqrt{\delta \mathrm{\Gamma }_0(T)}}{4}}\right|^2`$
$`+\left|{\displaystyle \frac{2\sqrt{\delta \mathrm{\Gamma }_0(\tau )}+e^{i\theta }\sqrt{\delta \mathrm{\Gamma }_0(\tau T)}+e^{i\theta }\sqrt{\delta \mathrm{\Gamma }_0(\tau +T)}}{4}}\right|^2,`$ (3)
with

$`\mathrm{\Gamma }_0(\tau )`$ $`=|ϵ|^2\left({\displaystyle \frac{F}{F_0}}\right)^2e^{\mathrm{\Omega }_c|\tau |}{\displaystyle \frac{\mathrm{sin}^2[(2N+1)\mathrm{\Omega }_F\tau /2]}{\mathrm{sin}^2(\mathrm{\Omega }_F\tau /2)}},`$ (4)
$`\delta `$ $`={\displaystyle \frac{4|ϵ|^2}{\mathrm{\Omega }_c^2}}.`$ (5)
Here $`ϵ`$ is the single-pass parametric amplitude gain; $`F`$ and $`F_0`$ are the finesse of the OPO with and without loss, respectively; $`\mathrm{\Omega }_c`$ and $`\mathrm{\Omega }_F`$ are the bandwidth and free spectral range of the OPO, respectively; $`2N+1`$ is the number of the longitudinal modes in the OPO output. $`\mathrm{\Gamma }_0(\tau )`$ is a multi-peaked function of delay time $`\tau `$. The width of the peaks is about $`\tau _r/(2N+1)`$. We assume that $`T`$ is nearly equal to $`\tau _r/2`$ but the difference between $`T`$ and $`\tau _r/2`$ is longer than the width of the peaks. This assumption is feasible in our experiment. This results in the following approximations:

$`\mathrm{\Gamma }_0(\tau )\mathrm{\Gamma }_0(\tau \pm T)`$ $`0,\mathrm{\Gamma }_0(\tau T)\mathrm{\Gamma }_0(\tau +T)0,\mathrm{\Gamma }_0(\pm T)0.`$ (6)
In addition, $`\delta `$ is much smaller than one when the OPO is operated far below threshold. Therefore, Eq. (3) can be approximated as follows:

$`\mathrm{\Gamma }(\tau ){\displaystyle \frac{\mathrm{\Gamma }_0(\tau )\mathrm{cos}^2\theta }{4}}+{\displaystyle \frac{\mathrm{\Gamma }_0(\tau T)+\mathrm{\Gamma }_0(\tau +T)}{16}}+{\displaystyle \frac{\delta \mathrm{\Gamma }_0(0)}{4}}.`$ (7)
The first term in the right-hand side of Eq. (7) corresponds to two-photon interference in Case 1. This is a periodic function of $`\theta `$ with a period $`\pi `$, which is a half of that of classical interference. It means that two-photon interference in Case 1 is nonclassical. The second term in the right-hand side of Eq. (7), which is constant with regard to $`\theta `$, corresponds to coincidence counts in Case 2. These two terms are due to correlated two photons. The last term in the right-hand side of Eq. (7) corresponds to the contribution from higher photon-number states than two. This term is not negligible when the pump power of the OPO is relatively high. As discussed in Ref. Goto , the coincidence rate measured in experiments is an average of the correlation function over the resolving time, $`T_R`$, of detectors. According to Ref. Goto , the coincidence counts measured in this experiment will become

$`\mathrm{\Gamma }_c(\tau )=C_1\left[4\mathrm{\Gamma }_c^{(0)}(\tau )\mathrm{cos}^2\theta +\mathrm{\Gamma }_c^{(0)}(\tau T)+\mathrm{\Gamma }_c^{(0)}(\tau +T)\right]+C_2,`$ (8)
with

$`\mathrm{\Gamma }_c^{(0)}(\tau )=e^{\mathrm{\Omega }_c|\tau \tau _0|}{\displaystyle \underset{n}{}}\left(1+{\displaystyle \frac{2|\tau n\tau _r\tau _0|\mathrm{ln}2}{T_R}}\right)\mathrm{exp}\left({\displaystyle \frac{2|\tau n\tau _r\tau _0|\mathrm{ln}2}{T_R}}\right).`$ (9)
Here $`C_1`$ and $`C_2`$ are constants and $`\tau _0`$ is an electric delay. The lines in Fig. 3 are fits to Eq. (8). Fitting parameters are two constants, $`C_1`$ and $`C_2`$, and the phase difference $`\theta `$. Constant parameters are set as follows: $`\tau _0=47.5`$ ns; $`\tau _r=2.07`$ ns; $`T_R=280`$ ns; $`\mathrm{\Omega }_c/(2\pi )=11`$ MHz. The range of the data used for the fitting is from 18 ns to 77 ns, while the range plotted in Fig. 3 is from 41 ns to 54 ns. The fitting in Fig. 3 is good. The term, $`C_2`$, independent of the delay time, $`\tau `$, is mainly due to the last term in the right-hand side of Eq. (7). In our experiment, $`C_2`$ is comparable to $`C_1`$. It means that the pump power of the OPO is relatively high and the contribution from higher photon-number states than two is not negligible.

The phase determined from the fitting is plotted in Fig. 4 against the phase locked experimentally. The inclination of the line in Fig. 4 is unity. The error bars are estimated from the fluctuation of the phase. The deviations of the circles from the line are probably due to the fluctuation of the phase difference, which is shown by error bars in Fig. 4, and due to an imperfect visibility, which makes larger the deviations around $`\theta =0`$, $`\pi /2`$, and $`\pi `$. Taking these points into consideration, the correspondence between theory and experiment seems not to be bad. From these fitting results, it is concluded that our experimental results can be explained by Eq. (8), which is derived from Eq. (7). Since Eq. (7) includes a term corresponding to nonclassical interference, it is also concluded that nonclassicality has appeared in this experimental results.

In conclusion, we have observed two-photon interference of multimode two-photon pairs produced by an OPO for the first time with an unbalanced interferometer. This two-photon interference has been dependent on delay time. The experimental results have been explained theoretically. Nonclassical feature of the two-photon interference has been also discussed.


format: text

text:
# Droplet Spreading on Heterogeneous Surfaces using a Three-Dimensional Lattice Boltzmann Model

## 1 Introduction

Wetting processes, such as the spreading of a droplet over a surface, have attracted the attention of scientists for a long time . A great deal is understood about the wetting behaviour of equilibrium droplets. However less is known about the dynamics of these systems, a problem of considerable industrial relevance with the advent of ink-jet printing. The droplets involved in printing typically have length scales of microns. Experimental work on such mesoscopic droplets is difficult and expensive because of the length- and time-scales involved. Therefore there is a need for numerical modelling both to investigate the physics and to help design and interpret the experiments.

Lattice Boltzmann models are a class of numerical techniques ideally suited to probing the behaviour of fluids on mesoscopic length scales . Several lattice Boltzmann algorithms for a liquid-gas system have been reported in the literature . They solve the Navier-Stokes equations of fluid flow but also input thermodynamic information, typically either as a free energy or as effective microscopic interactions. They have proved successful in modelling such diverse problems as fluid flows in complex geometries , two-phase models , hydrodynamic phase ordering  and sediment transport in a fluid .

Here we show that it is possible to use a lattice Boltzmann approach to model the spreading of mesoscale droplets and, in particular, to illustrate how a droplet spreads on a substrate comprising of hydrophilic and hydrophobic stripes.

We consider a one-component, two-phase fluid and use the free energy model originally described by Swift et al.  with a correction to ensure Galilean invariance . The advantage of this approach for the wetting problem is that it allows us to tune equilibrium thermodynamic properties such as the surface tension or static contact angle to agree with analytic predictions. Thus it is rather easy to control the wetting properties of the substrate. Three dimensional simulations of spreading on smooth and rough substrates have previously been reported in  using a different lattice Boltzmann algorithm.

The paper is organised as follows. First we summarise the main features of the lattice Boltzmann approach. The model is validated by showing the consistency of the measured equilibrium contact angle with that predicted by Young’s law and by measuring the base radius of the spreading droplet as a function of time obtaining, as expected, a power law growth. We show that when the reduced base radius is plotted as a function of reduced time the data fall on a universal curve for several values of surface tension and viscosity.

We then consider spreading on a heterogeneous substrate consisting of alternate hydrophobic and hydrophilic stripes. We find that the spreading velocity is anisotropic and that the final droplet shape reflects the wetting properties of the underlying substrate. Finally, a conclusion suggests extensions to the work presented here.

## 2 Simulating spreading

### 2.1 The lattice Boltzmann model

The lattice Boltzmann approach solves the Navier-Stokes equations by following the evolution of partial distribution functions $`f_i`$ on a regular, $`d`$-dimensional lattice formed of sites $`𝐫`$. The label $`i`$ denotes velocity directions and runs between $`0`$ and $`z`$. $`DdQz+1`$ is a standard lattice topology classification. The $`D3Q15`$ lattice topology we use here has the following velocity vectors $`𝐯_i`$: $`(0,0,0)`$, $`(\pm 1,\pm 1,\pm 1)`$, $`(\pm 1,0,0)`$, $`(0,\pm 1,0)`$, $`(0,0,\pm 1)`$ in lattice units.

The lattice Boltzmann dynamics are given by

$$f_i(𝐫+\mathrm{\Delta }t𝐯_i,t+\mathrm{\Delta }t)=f_i(𝐫,t)+\frac{1}{\tau }\left(f_i^{eq}(𝐫,t)f_i(𝐫,t)\right)$$
(1)
where $`\mathrm{\Delta }t`$ is the time step of the simulation, $`\tau `$ the relaxation time and $`f_i^{eq}`$ the equilibrium distribution function which is a function of the density $`n=_{i=0}^zf_i`$ and the fluid velocity $`𝐮`$ defined through the relation $`n𝐮=_{i=0}^zf_i𝐯_i`$.

The relaxation time tunes the kinematic viscosity as

$$\nu =\frac{\mathrm{\Delta }𝐫^2}{\mathrm{\Delta }t}\frac{C_4}{C_2}(\tau \frac{1}{2})$$
(2)
where $`\mathrm{\Delta }𝐫`$ is the lattice spacing and $`C_2`$ and $`C_4`$ are coefficients related to the topology of the lattice. These are equal to $`3`$ and $`1`$ respectively when one considers a $`D3Q15`$ lattice (see  for more details).

It can be shown  that equation (1) reproduces the Navier-Stokes equations of a non-ideal gas if the local equilibrium functions are chosen as

$`f_i^{eq}`$ $`=`$ $`A_\sigma +B_\sigma u_\alpha v_{i\alpha }+C_\sigma 𝐮^2+D_\sigma u_\alpha u_\beta v_{i\alpha }v_{i\beta }+G_{\sigma \alpha \beta }v_{i\alpha }v_{i\beta },i>0,`$
$`f_0^{eq}`$ $`=`$ $`n{\displaystyle \underset{i=1}{\overset{z}{}}}f_i^{eq}`$ (3)
where Einstein notation is understood for the Cartesian labels $`\alpha `$ and $`\beta `$ (i.e. $`v_{i\alpha }u_\alpha =_\alpha v_{i\alpha }u_\alpha `$) and where $`\sigma `$ labels velocities of different magnitude.

The coefficients $`A_\sigma `$, $`B_\sigma `$, $`C_\sigma `$, $`D_\sigma `$ and $`G_\sigma `$ are chosen so as to satisfy the relations

$`{\displaystyle \underset{i}{}}f_i^{eq}`$ $`=`$ $`n,`$
$`{\displaystyle \underset{i}{}}f_i^{eq}v_{i\alpha }`$ $`=`$ $`nu_\alpha ,`$
$`{\displaystyle \underset{i}{}}f_i^{eq}v_{i\alpha }v_{i\beta }`$ $`=`$ $`P_{\alpha \beta }+nu_\alpha u_\beta +\nu \left(u_\alpha _\beta n+u_\beta _\alpha n+u_\gamma _\gamma n\delta _{\alpha \beta }\right),`$
$`{\displaystyle \underset{i}{}}f_i^{eq}v_{i\alpha }v_{i\beta }v_{i\gamma }`$ $`=`$ $`{\displaystyle \frac{c^2n}{3}}(u_\alpha \delta _{\beta \gamma }+u_\beta \delta _{\alpha \gamma }+u_\gamma \delta _{\alpha \beta })`$ (4)
where $`P_{\alpha \beta }`$ is the pressure tensor, $`c`$ is defined to be $`\mathrm{\Delta }𝐫/\mathrm{\Delta }t`$ and the last term of the third expression in equation (4) is included to ensure Galilean invariance.

Considering a $`D3Q15`$ lattice and a square-gradient approximation to the interface free energy ($`\kappa (_\alpha n)^2/2`$, a possible choice of the coefficients is

$`A_\sigma `$ $`=`$ $`{\displaystyle \frac{w_\sigma }{c^2}}\left(p_b{\displaystyle \frac{\kappa }{2}}(_\alpha n)^2\kappa n_{\alpha \alpha }n+\nu u_\alpha _\alpha n\right),`$
$`B_\sigma `$ $`=`$ $`{\displaystyle \frac{w_\sigma n}{c^2}},C_\sigma ={\displaystyle \frac{w_\sigma n}{2c^2}},D_\sigma ={\displaystyle \frac{3w_\sigma n}{2c^4}},`$
$`G_{1\gamma \gamma }`$ $`=`$ $`{\displaystyle \frac{1}{2c^4}}\left(\kappa (_\gamma n)^2+2\nu u_\gamma _\gamma n\right),G_{2\gamma \gamma }=0,`$
$`G_{2\gamma \delta }`$ $`=`$ $`{\displaystyle \frac{1}{16c^4}}\left(\kappa (_\gamma n)(_\delta n)+\nu (u_\gamma _\delta n+u_\delta _\gamma n)\right)`$ (5)
where $`w_1=1/3`$, $`w_2=1/24`$, $`\kappa `$ is a parameter related to the surface tension and $`p_b`$ is the pressure in the bulk which is defined below. One can show  that the pressure tensor can be written as

$$P_{\alpha \beta }=\left(p_b\frac{\kappa }{2}(_\gamma n)^2\kappa n_{\gamma \gamma }n\right)\delta _{\alpha \beta }+\kappa (_\alpha n)(_\beta n).$$
(6)
### 2.2 Wetting boundary conditions

In this paper, we will focus our attention on flat substrates normal to the $`z`$ direction. The derivatives in that direction should then be handled in such a way that the wetting properties of the substrate can be controlled. A boundary condition can be established using the Cahn model . He proposed adding an additional surface free energy $`\mathrm{\Psi }_c(n_s)=\varphi _0\varphi _1n_s+\mathrm{}`$ at the solid surface where $`n_s`$ is the density at the surface. Neglecting the second order terms of $`\mathrm{\Psi }_c(n)`$ and minimizing $`\mathrm{\Psi }_b+\mathrm{\Psi }_c`$ (where $`\mathrm{\Psi }_b`$ is the free energy in the bulk), a boundary condition valid at $`z=0`$ emerges

$$_zn=\frac{\varphi _1}{\kappa }.$$
(7)
Equation (7) is imposed on the substrate sites to implement the Cahn model in the lattice Boltzmann approach. Details are given in .

The Cahn model can be used to relate $`\varphi _1`$ to $`\theta `$ the contact angle defined as the angle between the tangent plane to the droplet and the substrate. Within the Cahn model the surface tension at the interfaces is given by

$`\sigma _{lg}`$ $`=`$ $`{\displaystyle _{n_g}^{n_l}}\sqrt{2\kappa W(n,T)}𝑑n`$
$`\sigma _{sl,sg}`$ $`=`$ $`\left|{\displaystyle _{n_s}^{n_{g,l}}}\sqrt{2\kappa W(n,T)}𝑑n\right|+\varphi _0\varphi _1n_s`$ (8)
where $`W(n,T)`$ is the excess free energy, $`\sigma _{lg}`$, $`\sigma _{sg}`$, $`\sigma _{sl}`$ are the surface tensions at the liquid-gas, solid-gas and solid-liquid interface respectively and $`n_s`$, $`n_l`$, $`n_g`$ are the densities at the substrate, of the liquid phase and of the gas phase respectively. Young’s law  gives a relation between the static contact angle and the surface tensions of the three phases

$$\sigma _{lg}\mathrm{cos}\theta =\sigma _{sg}\sigma _{sl}.$$
(9)
A convenient choice of bulk pressure is

$$p_b=p_c(\nu _p+1)^2(3\nu _p^22\nu _p+12\beta \tau _p)$$
(10)
where $`\nu _p=(nn_c)/n_c`$, $`\tau _p=(T_cT)/T_c`$ and $`p_c=1/8`$, $`n_c=3.5`$ and $`T_c=4/7`$ are the critical pressure, density and temperature respectively and $`\beta `$ is a constant typically equal to $`0.1`$. The excess free energy then becomes

$$W(n,T)=p_c(\nu _p^2\beta \tau _p)^2.$$
(11)
Inserting equation (11) into relation (8) and using equation (9), gives the relation between $`\varphi _1`$ and $`\theta `$

$$\varphi _1=2\beta \tau _p\sqrt{2p_c\kappa }\mathrm{sign}(\theta \frac{\pi }{2})\sqrt{\mathrm{cos}\frac{\alpha }{3}\left(1\mathrm{cos}\frac{\alpha }{3}\right)}$$
(12)
where $`\alpha =\mathrm{cos}^1(\mathrm{sin}^2\theta )`$ and the function $`\mathrm{sign}`$ returns the sign of its argument.

We impose a no-slip boundary condition on the velocity. Because a collision takes place on the boundary the usual bounce-back condition must be extended to ensure mass conservation (see  for a wider discussion). This is done by a suitable choice of the rest field, $`f_0`$, to correctly balance the mass of the system.

## 3 Spreading on a homogeneous surface

We consider a $`80\times 80\times 40`$ lattice on which a spherical drop of radius $`R_0=16`$ just touches a flat surface at $`z=0`$. Unless otherwise specified the temperature is $`T=0.4`$ which leads to two phases of density $`n_l=4.128`$ and $`n_g=2.913`$. Fig. 1 shows how the droplet evolves in time to reach an equilibrium shape with contact angle $`60^{}`$.

We first present a check on the accuracy of the equilibrium properties of the model. Fig. 2 reports a comparison between two methods of measuring the contact angle. $`\theta ^y`$ is the contact angle obtained from equation (9) with the surface tensions measured at equilibrium. $`\theta ^g`$ is the contact angle measured from the profile of the simulated droplet once equilibrium is reached. The agreement is good. Small errors results from the difficulty of a direct measurement of the contact angle on a discrete lattice.

The shape of the area formed by the contact of a droplet with a homogeneous substrate is a disk. Its radius $`R_c`$ is a quantity which is rather simple to measure and has consequently attracted the attention of many scientists, see  and references therein. The time evolution of $`R_c`$ has been found to follow a power law $`R_c=mt^{n/2}`$. The exponent $`n`$ has been widely reported in the literature but with no consistent result. Marmur  in his review lists exponents between $`0.06`$ and $`0.6`$. The value of $`m`$ appears to be related to the droplet volume.

Fig. 3 shows the time evolution of $`R_c`$ for different values of the viscosity and the surface tension. The curves correspond to a value $`n=0.56`$ which is within the range reported in the literature. The power law is independent of the surface tension and the viscosity.

Indeed if the evolution curves are plotted as a function of the dimensionless time

$$tt^{}=\frac{\sigma _{lg}}{\eta R_0}t$$
(13)
the data collapses onto a single curve as shown in fig. 3(b). Experimental data taken from  shows similar behaviour.

## 4 Spreading on heterogeneous surfaces

Almost any surface will contain physical and chemical inhomogeneities which will affect the spreading of a mesoscopic droplet. It has recently become feasible to fabricate surfaces with well-defined chemical properties on micron length scales and it is becoming possible to perform well-controlled experiments which probe the behaviour of mesoscopic droplets on chemically and physically heterogeneous substrates. Thus it is particularly interesting to develop techniques to model the effect of these surfaces on the spreading properties of a droplet.

One of the simplest heterogeneous surfaces can be formed by alternating stripes of materials with different wetting properties. The static properties of droplets on such substrates have been discussed . However less attention has been paid to the dynamics of the spreading.

In this section we consider heterogeneous surfaces formed by alternating hydrophilic and hydrophobic stripes. They are characterised by widths $`w_{phi}`$,$`w_{pho}`$ and contact angles $`\theta _{phi}`$, $`\theta _{pho}`$ respectively. Fig. 4 presents the behaviour of a three-dimensional droplet spreading on such a surface with $`\theta _{phi}=50^{}`$, $`\theta _{pho}=110^{}`$, $`w_{phi}=6`$, $`w_{pho}=5`$. The droplet has an initial radius $`R_0=18`$.

It is apparent from the figure that the behaviour of the droplet depends on whether it is on a hydrophobic or a hydrophilic stripe. The equilibrium shape of the contact line shown in fig. 4(b) reflects the pattern of the underlying substrate which is comparable to that found in laboratory experiments .

The time evolution of the contact line is also shown in fig. 4(b). Note that its velocity decreases smoothly in the $`y`$-direction parallel to the stripes but not in the $`x`$-direction where it moves faster on the hydrophilic than on the hydrophobic stripes. Note also that the droplet remains symmetric about an axis perpendicular to the stripes but that the shape becomes asymmetric about an axis parallel to the stripes, depending upon the initial position of the center of the droplet.

Observation of the movement of the contact line in the $`x`$-direction shows that in a hydrophilic region the contact angle tends to decrease and the velocity of the contact line increase. When the contact line reaches the boundary its progress is stopped and the contact angle increases until it is large enough to cross the hydrophobic stripe.

It has been proposed that an equilibrium droplet on such a surface has a spatially averaged contact angle following Cassie’s law

$$\mathrm{cos}\theta =p_{phi}\mathrm{cos}\theta _{phi}+p_{pho}\mathrm{cos}\theta _{pho}$$
(14)
where $`p_{phi}`$ and $`p_{pho}`$ are the proportion of the substrate area which are hydrophilic or hydrophobic respectively. However, this relation is not universally accepted. In particular it has been argued that there should be a dependence on the relative size of the droplet and the surface stripes . Fig. 4(c)-(d) show characteristic angles for the droplet considered here. Their average is $`76.5^{}`$ which is close to the one predicted by Cassie’s law, $`78.7^{}`$.

## 5 Conclusion

We have used a three-dimensional lattice Boltzmann algorithm to model the spreading of a mesoscopic droplet. By incorporating the Cahn theory of wetting into the simulation we obtain a way of easily tuning the contact angle of the droplet on the substrate. This gives us the ability to simulate spreading on both homogeneous and heterogeneous surfaces.

The approach provides a well-controlled way of investigating the dependence of the spreading on such properties as the droplet volume, contact angles and the substrate geometry. Further work is in progress to systematically determine how these parameters affect the velocity and shapes of the spreading droplets. It would also be interesting to investigate the effect of physical inhomogeneities on the spreading and to consider a droplet spreading on a porous surface.

A particular aim of the work will be to compare the results to forthcoming experiments on substrates which have chemical patterning on mesoscopic length scale. This will allow increased understanding of both the experimental results and the model assumptions. For example we assume, as is the standard practice, no-slip boundary conditions on the velocity. These may not be appropriate on short length scales near a contact line. Moreover the liquid-gas density difference in lattice Boltzmann models is very small compared to real fluids and it is important to undertake further work to assess the effect of this on the modelling results.

## Acknowledgments

We thank D. Bucknall, J. Leopoldes and S. Willkins for helpful discussions. Supercomputing resources were provided by the Oxford Supercomputing Centre. AD acknowledges the support of the EC IMAGE-IN project GR1D-CT-2002-00663.


text:
# Effects of dark-matter density profiles on bulge formation in very late-type galaxies

## 1 Introduction

In our previous work (Fu, Huang & Deng 2003a, (FHD03)), we managed to construct a set of models for simulating the dynamical evolution of super star clusters (SSCs) embedded in dark matter (DM) halo with central cusp, motivated by observations in nuclear regions of galaxies. Among these observations we would mention the discovery of the presence of a black hole with intermediate mass (IMBH) in a young SSC located at some distance from the nucleus of M82 (Matsumoto et al. 2000). While another possible IMBH has been found in an evolved SSC located in the nucleus of M33 (Gebhardt et al 2001). These discoveries strongly imply the role of dynamical evolution of SSCs in the physical process causing the $`M_{\mathrm{BH}}M_{\mathrm{bulge}}`$ relation (see, Fu, Huang & Deng 2003b), as well as in bulge formation in late-type galaxies.

The simulations in FHD03 have yielded bulges that are similar in many aspects to the observational ones. In particular, the derived surface density profiles can be well fitted by an exponential structure with nuclear cusps, which is consistent with Hubble Space Telescope (HST) observations (Carollo 1999). The progress in observations (e.g. Marchesini et al. 2002; and the references therein) shows, however, that the cuspy-core DM profiles, such as the NFW profile, are not adequate for a large fraction of dwarf galaxies which are dominated by the DM halos, and suggested that the soft-core density profiles, the Burkert profile for example, is more suitable to these galaxies. It would then be very much instructive to perform simulations for the dynamical evolution of SSCs with both the cuspy- and the soft-core DM density profiles for the bulge formation, with the hope of constraining the bulge-formation model or the properties of the dark matter in late-type galaxies. That is what we would present in this paper.

## 2 Models

In view of what we described in introduction, we would then model the sinking and tidal stripping of SSCs in DM halos depicted with the NFW or the Burkert density profiles in this paper so as to derive something important as compared with observations. The part of the models in FHD03 are adjusted accordingly and are described in the following subsections. The readers are referred to FHD03 for detailed description about the other part of the models which remain unchanged in this work, including dynamical friction and tidal stripping.

### 2.1 SSC model

Recently, Parodi & Binggeli (2003) studied the spatial distribution of bright star-forming complexes in a volume-limited sample of 72 very late-type galaxies located within 10 Mpc. They derived the projected radial number distribution $`R\mathrm{exp}(R/R_l)`$, with median scale length $`R_l0.6kpc`$ and median number of bright star-forming lumps 17, respectively. Taking this spatial distribution with its median values as our SSCs’ radial distribution, rather than 1/R function inferred from 3 galaxies in FHD03, we generate a set of SSCs over the whole galaxy instead of setting both inner and outer boundary for its distribution. As in FHD03, these SSCs are assigned with local circular speed. Both the initial mass function of SSCs and the SSC itself remain unchanged as those in FHD03.

### 2.2 Background

Since very-late type galaxies are DM dominated at all radii (Marchesini et al., 2002; Persic, Salucci & Stel, 1996), we only take, as in FHD03, a DM halo as the initial background for the bulge-formation process. In light of considering the effects of DM halo models on the bulge formation, we adopt the NFW (Navarro, Frenk & White, 1997) and the Burkert (1995) profiles as the cuspy- and the soft-core density profiles, respectively. The mass of the dark matter halo is taken to be $`M_{200}=10^{11}M_{}`$.

The NFW profile can be written as (Binney & Merrifield, 1998)

$$\rho _{NFW}(r)=\frac{M_{0h}}{r(a_h+r)^2},$$
(1)
where $`r`$ is the distance from the halo center, and, $`a_h`$ and $`M_{0h}`$ are two parameters. With a $`\mathrm{\Lambda }`$CDM cosmological model (Jing, 2000), the NFW halo is scaled to $`M_{200}=10^{11}M_{}`$, which gives

$$a_h=9.1kpc,M_{0h}=5.3\times 10^{10}M_{}.$$
(2)
The Burkert profile writes

$$\rho _{Burkert}(r)=\frac{\rho _0r_0^3}{(r+r_0)(r^2+r_0^2)},$$
(3)
where the central density $`\rho _0`$ is taken to be $`0.05M_{}/pc^3`$ following Marchesini et al. (2002), and the scale radius $`r_0`$ of 4 kpc is computed from the definition of $`M_{200}`$.

## 3 Results and Discussions

### 3.1 Properties of formed bulges and formation rates

Fig 1 shows the evolution of surface density profiles of simulated individual galaxies obtained with the NFW and the Burkert density profiles, respectively. Obviously, the displayed profiles in Fig 1 fall into two categories, one with high-density central component and the other without. The latter ones basically correspond to the cases where no SSC contributes its mass to the bulge area, or to the cases of being at the pre-stage of bulge-forming. On the contrary, the former ones correspond to the cases where there is some mass contribution from SSCs down to less than 1 pc from the halo center, and the profiles clearly show that bulges are well-formed.

The mean surface density profiles of the formed bulges derived with the NFW or the Burkert profiles at 3 Gyr are illustrated in Fig 2, together with the curves depicting the fitting to the model (Carollo & Stiavelli 1998; FHD03),

$$\sigma (R)=\sigma _0\mathrm{exp}(1.678\frac{R}{R_e})+\sigma _1(1+\frac{R_c}{R})^\gamma \mathrm{exp}(\frac{R}{R_s}).$$
(4)
As obtained in FHD03, the formed bulges, either with the NFW or the Burkert density profiles, do share the same characteristics with the HST observational results (Carollo 1999), i.e. the general presence of central cusps on top of the exponential bulges.

The half-mass radii, $`R_e`$, of the formed exponential bulges for individual galaxies, derived from the fitting with equa (4), range over 100pc to 1000pc. Fig 3a,b show the histograms of the scale radii $`R_e`$ of individual bulges formed at 3 Gyr. Either with the NFW or the Burkert density profile, the scale radii $`R_e`$ distributed over a few hundred parsecs are in accord with observed results by the HST (Carollo 1999).

The above analyses indicate that neither the NFW nor the Burkert profile make a difference to the formed bulges in view of their structure. Nevertheless, the fraction of simulated galaxies with no bulges formed turn out to be quite different as Fig 1 illustrated. This situation is clearly manifested in Fig 4. The bulge formation rates are quite low in simulated galaxies derived with the Burkert density profile, about 30% and 70% at the age of 1 and 4 Gyr, respectively. On the contrary, the formation rates derived with the NFW profile are above 97% and about 100% at 1 and 4 Gyr, respectively.

Observations on searching for bulges in very late-type galaxies are at the very beginning indeed (see, e.g. Böker et al. 2003; Carollo 1999). These researches provide the bulge detection rates of about 50%, compatible with the bulge formation rates derived from the simulations with the Burkert profile. Though it would be natural to conclude that the DM density profiles suitable to very late-type dwarfs are soft-core ones, we should be cautious at the moment. For one thing, these surveys are based on samples of small number of sources, 19 galaxies at the most. On the other hand, Matthews & Gallagher (1997) found that bulges of extremely late-type spirals are often very weak.

### 3.2 The mass relation of bulges to their nuclear cusps

Though Carollo et al. (2002) discovered that each exponential bulge in their sample hosts a stellar nucleus in its center, they claimed no necessarily causal connection between these two components. In our simulations (FHD03, and this work), however, every formed bulge is always paired with its central, stellar cusp. The stripped stars from SSCs contribute to the bulge formation, while the remaining, undissociated SSCs form the nuclear cusps. Needless to say, to confirm the causal connection between nuclear cusps and the exponential bulges would shed much light on the evolution of galaxies.

Recent work by Balcells et al. (2003) seems to confirm the existence of the causal connection between the paired bulges and nuclear cusps, where a tight luminosity correlation in H-band has been reported. In order to see if there is a corresponding correlation between the masses of our simulated bulges and their paired cusps, we have measured the scale sizes of nuclear cusps quantified by their FWHMs, following Carollo et al. (2002). Considering the minimum angular size of 0.″03 that WFPC2 on board HST can resolve and the distance ranges of 12 and 37 Mpc for observed galaxies, we have averaged the simulated nuclear cusps over central radius of 2 or 5 pc so that the measured FWHMs of simulated nuclear cusps are comparable to the observed ones. Fig 3c,d show the FWHMs thus obtained of the simulated nuclear cusps, which conforms with observed results (Carollo et al. 2002), ranging from a few pc up to 20 pc.

The derived mass relation between simulated bulges and their paired nuclear cusps is illustrated in Fig 5, where the correlations are obviously reached either for the NFW or the Burkert profile at 3 Gyr. The consistency between the observed and simulated mass correlation provides further evidence favorable to the proposed scenario (FHD03 and this work) for bulge formation in very late-type galaxies.

Besides, it would be worth noticing what Fig 5 reveals that the formed bulges and nuclear cusps are more heavier in case of the NFW density profile, with median values of about $`2\times 10^7M_{}`$ and $`1.5\times 10^6M_{}`$ respectively at 3 Gyr. On the other hand, the relevant values derived in the Burkert profile are about $`8\times 10^6M_{}`$ and $`6\times 10^5M_{}`$ respectively. This difference is mainly made by the different local densities and density variation of DM within the circumnuclear region in the cuspy- and soft-core profiles. Further observations are surely needed to confirm this kind of disparity.

### 3.3 Formation of nuclear star clusters without bulges

Now it is believed that the occurrence of compact star clusters at nuclei of very late-type spirals is a common phenomenon (e.g. Böker et al. 2002). A fraction of these clusters are located in bulge-less galaxies, with ages roughly ranging from $`10^7`$ to $`10^9yr`$ (e.g. Walcher et al. 2003).

As we have shown (FHD03 and this work), our proposed scenario provides a promising mechanism for the formation of massive, nuclear star clusters in galaxies having detected bulges, with masses compatible with observed values ranging very crudely from a few times $`10^5`$ to $`10^7M_{}`$ (e.g. Matthews & Gallagher 2002).

The existence of nuclear star cluster in bulgeless galaxies, however, remains a mystery. This phenomenon is closely related to several hot topics, especially the properties of DM halos in galaxies. We have then made records of simulations for spatial distances of sinking SSCs during 1Myr and 3 Gyr, either in case of the NFW or Burkert profile, to see if some massive SSCs could survive tidal stripping so as to sink into center of simulated galaxies before bulge formed.

The snapshots of this kind are illustrated in Fig 6 for DM of NFW and Burkert profile, respectively. It is very intriguing to find from Fig 6 that massive, nuclear star clusters do form in NFW profile at 10Myr and in Burkert profile around 100Myr with no discernible bulges. Indeed, this situation has already been manifested in Fig 1 at corresponding ages. After 10Myr, the probability to have nuclear star cluster with no bulge formed becomes low for DM depicted with NFW density profile, while it is low before around 10Myr in case of Burkert profile. After 1 Gyr, this probabilities are quite low in both cases, especially in case of NFW profile due to very high bulge formation rate, above 97%.

It seems then the scenario we proposed could be a promising mechanism to form nuclear star clusters in bulge-less galaxies. Before drawing a conclusion about the suitability of cuspy-core density profiles for very late-type dwarfs, however, the models presented in FHD03 do need inclusion of more physical processes, e.g. the reaction of DM halo to the sinking of SSCs. This process might have strong effect on changing the DM density profiles (El-Zant et al. 2001).

This work is supported by NKBRSF G19990754, and NSFC.


text:
# Critical Behavior of the Ferromagnetic Ising Model on a Sierpiński Carpet: Monte Carlo Renormalization Group Study

## 1 Introduction

Since Mandelbrot attracted people’s attention to the fractals in the 60’s and 70’s , scientists began to take auto-similarity and scaling invariance as fundamental rules of nature. In the description of second order phase transitions, Widom’s homogeneity hypothesis  and Kadanoff’s real space renormalization group scheme  are based on the invariance of the physical behavior under any change of scale at the criticality where the correlation length is divergent. The validity of these hypothesis and scheme has been verified in different systems with an integer dimension which possess a translational symmetry . As a matter of fact, a system with translational symmetry is auto-similar and can be considered as a particular case of fractal. Fractals are natural candidates to represent systems with non-integer dimensions. It turns out to be of fundamental relevance to know if the hypothesis and the scheme mentionned above still work on a general fractal system, where the translational symmetry is lost and replaced by scale invariance.

The critical behavior of the Ising model on fractal lattices has been firstly studied by Gefen and his co-workers in the 80’s . Their works were based on the Migdal-Kadanoff bond-moving renormalization method and showed that the topological features of the fractals plays an important role in the determination of the critical behavior. As a main result, they found that a phase transition at non-zero temperature can occur only if the fractal has an infinite ramification order. Later, Bonnier et al. used an alternative decimation method  and high-temperature expansions  to study the phase transition on various fractals, namely Sierpiński carpets. They also found that the critical exponents depend on the geometrical properties of the fractal. In the above studies, the applied methods are approximative and the results are not always consistent. Moreover, in these theoretical analyses, the spins were placed at the corners of the occupied squares of the Sierpiński carpets; in these cases, the number of spins doesn’t follow a power law of the lattice size and the Hausdorff dimension should not be expected to enter in the description of the critical phenomena on such systems. Recently, due to the progress in the simulation methods and the growth of the computer power, the critical behavior of the Ising model on fractals of Hausdorff dimension $`d_f`$ between 1 and 3 has been studied numerically in a much more precise way . The results showed that the finite-size scaling (FSS) analysis works in the case of fractals, although the convergence towards the thermodynamical limit can be very slow when $`d_f<2`$, and that the hyperscaling law $`d_f=2\beta /\nu +\gamma /\nu `$ is verified. Moreover, discrepancies with the predictions of the $`ϵ`$-expansions  were observed. The universality of phase transitions on fractals is said to be weak . It is worth noting that the spins were placed at the center of the occupied squares of the Sierpiński carpets in these studies; consequently, the number of the spins increases as a power law of the lattice size with an exponent equal to the Hausdorff dimension $`d_f`$. In this case, $`d_f`$ takes the place of the space dimension in Widom’s homogeneity hypothesis because of the conservation of the total free energy under a change of length scale; the generalization of this hypothesis in the case of fractals is quite straightforward .

Monte Carlo renormalization group (MCRG) has been shown to be a powerful tool in the study of critical phenomena . It combines Monte Carlo (MC) simulation techniques with the analysis of the real space renormalization group (RSRG). It has been, so far, applied to systems with translational symmetry. A RSRG analysis involves an infinite number of couplings when dealing with the Hamiltonian. Theoreticians usually truncate the number of couplings to a finite one in order to make the calculations tractable. It, therefore, makes RSRG an approximate method and the accuracy of the results is linked to the number of the couplings considered. Even if only few couplings are considered, the calculations associated with the process of decimation remain tedious. With the aid of MC simulation, a renormalization-group calculation with a larger number of couplings becomes possible and the results are expected to be more reliable.However, some crucial problems are encountered in a MCRG study:

1. The critical slowing down:
Hypothetically, the accuracy of the thermodynamical averages can be improved by performing a larger number of MC steps to reduce the statistical errors. In fact, due to the divergence of the correlation length at the criticality, local MC algorithms (for e.g., the Metropolis one) suffer from critical slowing down, which hinders the computation of accurate thermodynamical averages. Fortunately, the use of cluster algorithms, for e.g., the Swendsen-Wang algorithm  and the Wolff algorithm , enables to improve significantly the efficiency of the simulations.
2. The finite size effects:
RSRG method involves an infinite size lattice, whereas computer simulations can only be performed on finite size lattices. One should perform the simulation on several large systems and then extrapolate the results to the infinite one.

The MCRG method, according to which one plays the game of changing the length scale, should be a natural approach to explore the critical behavior of phase transitions on hierarchical lattices. The purpose of this paper is to study the critical behavior of the Ising model on a Sierpiński carpet using the MCRG method. As far as we know, it is the first time that the MCRG method is used in the case of fractal lattices. The paper is organized as follows : In Sec. II, we present the model and recall briefly the MCRG method. We explain how the simulation is setup in Sec. III. Sec. IV is devoted to the calculation of the critical temperature $`T_c`$, the temperature eigen-exponent $`y_t`$, and the magnetic eigen-exponent $`y_h`$. A FSS analysis is set out in Sec. V for a consistency check.

## 2 Sierpiński carpet and MCRG method

A 3 by 3 square lattice with its central subsquare removed is chosen as “generating cell”. The Sierpiński carpet we deal with is constructed iteratively from this generating cell:

1. We take the generating cell as the carpet at the first iteration step and denote it $`SC(3,1,1)`$.
2. The carpet at the $`(k+1)`$-th iteration step $`SC(3,1,k+1)`$ is constructed by enlarging three times the size of the carpet at the $`k`$-th iteration step and replacing each enlarged un-removed subsquare by the whole generating cell. This carpet becomes a “true” fractal in the mathematical sense when $`k`$ tends to infinity.

The size of the lattice $`SC(3,1,k)`$ is equal to $`L_k=3^k`$. A spin is set at each center of an un-removed subsquare (so-called occupied square). The number of spins on $`SC(3,1,k)`$ is therefore equal to $`N_k=8^k`$. We notice that $`N_k`$ is equal to $`L_{k}^{}{}_{}{}^{d_f}`$ where $`d_f=\mathrm{log}8/\mathrm{log}31.892789`$ is the Hausdorff dimension. The critical behavior of the ferromagnetic Ising model on this Sierpiński carpet has been recently studied independently by P. Monceau et al.  and J. Carmona et al.  by means of MC simulation. Thus, estimations of the critical temperature are available and the comparison between the two methods will be possible.

Let us briefly recall the MCRG method.The geometrical structure of a Sierpiński carpet remains invariant during a RSRG process involving a change in the length scale from 1 to $`b`$, provided that $`b`$ is equal to an integer power of the size of its generating cell. In our case, we set $`b`$ equal to 3 and increase the length scale by this factor at each renormalization step. The majority rule is used to decimate the spin blocks: a new spin is assigned to a given block according to the sign of the summation of the spin states of the 8 occupied sites in this block; if the sum is zero, a spin state $`+1`$ or $`1`$ is assigned to the block with the same probability. The reduced Hamiltonian after $`n`$ renormalization steps ($`n`$ is a positive integer) reads as

$$^{(n)}(\{K_\alpha ^{(n)}\}_{\alpha =1}^{\mathrm{}})=\underset{\alpha =1}{\overset{\mathrm{}}{}}K_\alpha ^{(n)}S_\alpha ^{(n)},$$
(1)
where $`\{K_\alpha ^{(n)}\}_{\alpha =1}^{\mathrm{}}`$ is the coupling set after $`n`$ renormalization steps and $`S_\alpha ^{(n)}`$’s are the conjugate lattice sums of the spin products on the reduced network. The matrix $`𝒯`$, which describes the flow of the couplings from the $`n`$-th renormalization step to the $`(n+1)`$-th one, is defined by

$$𝒯_{\alpha \beta }^{(n+1,n)}=\frac{K_\alpha ^{(n+1)}}{K_\beta ^{(n)}}$$
(2)
and satisfies the following relation:

$$\left[S\right]_{\gamma \beta }^{(n+1,n)}=\underset{\alpha =1}{\overset{\mathrm{}}{}}\left[S\right]_{\gamma \alpha }^{(n+1,n+1)}𝒯_{\alpha \beta }^{(n+1,n)},$$
(3)
where $`\left[S\right]_{\alpha \beta }^{(n,m)}S_\alpha ^{(n)}/K_\beta ^{(m)}`$ can be calculated by $`S_\alpha ^{(n)}S_\beta ^{(m)}S_\alpha ^{(n)}S_\beta ^{(m)}`$. One can hence calculate the matrix $`𝒯`$ by inverting the above relation. Two critical exponents enable to describe the static scaling properties of the system: the temperature eigen-exponent $`y_t`$ and the magnetic one $`y_h`$. They describe respectively the scaling behavior of the reduced temperature $`t=T/T_c1`$ (where $`T_c`$ is the critical temperature) and the external magnetic field $`h`$ near the criticality under a change of the length unit from $`1`$ to $`b`$:

$`t`$ $``$ $`t^{}=b^{y_t}t,`$
$`h`$ $``$ $`h^{}=b^{y_h}h.`$
$`b^{y_t}`$ and $`b^{y_h}`$ are associated with the two relevant directions of the renormalization flows and can be obtained by finding the largest eigenvalues of the $`𝒯`$ matrix in the even- and odd-coupling subspaces, respectively.

## 3 Simulation setup

Up to 4-spin couplings have been considered in our study and the interaction range has been restricted within a $`3\times 3`$ square. There are 25 even couplings and 11 odd couplings as shown in Fig.1(a) and Fig.1(b), respectively. The symmetry number obtained by rotating and reflecting a given coupling is indicated in the third column of the figure. These couplings are listed in decreasing order according to the importance factor proposed by Blöte et al. : $`F=(2^{s/2}\overline{r})^1`$, where $`\overline{r}`$ is the average distance between the $`s`$ spins. The value of $`F`$ is given in the right column of the figure. The simulations have been carried out at two different temperatures, $`T=1.4813`$ obtained by Carmona et al.  and $`T=1.4795`$ obtained by Monceau et al. . The size of the starting Sierpiński carpet varies between $`3^4=81`$ and $`3^8=6561`$. The simulation procedure is organized as follows:

1. An equilibrium spin configuration on the starting Sierpiński carpet $`SC(3,1,k)`$ is generated by the Wolff algorithm  at the simulation temperature.
2. The lattice sums conjugated to the 25 even couplings and to the 11 odd ones are calculated.
3. A decimation is done, according to the majority rule, by dividing the size of the Sierpiński carpet by 3 and the lattice sums associated with the renormalized carpet are calculated.
4. The step 3 is repeated until the size of the renormalized Sierpiński carpet is smaller than 9.

A data sample of the lattice sums is built up from $`10^5`$ steps of the above simulation, with periodic boundary conditions. At each simulation temperature, for the five different sizes of the starting Sierpiński carpet, 10 independent data samples are collected. The errors bars are estimated from a standard statistical analysis of these data.

## 4 MCRG Results

### 4.1 Critical temperature $`T_c`$

The size of the lattice after $`n+p`$ renormalization steps from the starting Sierpiński carpet of size $`L_{k+p}`$ is equal to the size of the lattice after $`n`$ renormalization steps from the starting Sierpiński carpet of size $`L_k`$. At the fixed point $`𝐊^{}`$, the thermodynamical averages of the lattice sums calculated from renormalized carpets with the same size should be independent of the size of the starting Sierpiński carpet. Let us assume that the initial coupling vector $`𝐊^{(\mathrm{𝟎})}`$ is close to $`𝐊^{}`$. The difference of the lattice sums on the reduced lattices of equal size, derived from the small displacement $`\delta 𝐊=𝐊^{}𝐊^{(\mathrm{𝟎})}`$, should satisfy the following relation:

$$\left(S_\alpha ^{(n+p)}_{L_{k+p}}S_\alpha ^{(n)}_{L_k}\right)|_{𝐊^{(\mathrm{𝟎})}}\underset{\beta }{}\left(\left[S_{L_{k+p}}\right]_{\alpha \beta }^{(n+p,0)}\left[S_{L_k}\right]_{\alpha \beta }^{(n,0)}\right)|_{𝐊^{(\mathrm{𝟎})}}\delta K_\beta ,$$
(4)
where the suffix $`L_{k+p}`$ or $`L_k`$ indicates that the physical quantity is obtained from the starting Sierpiński carpet of size $`L_{k+p}`$ or $`L_k`$. In the even coupling subspace, the nearest neighborhood coupling $`K_1`$ represents the inverse of the temperature. If we neglect the contribution from the directions other than the nearest neighborhood coupling and set $`\alpha =1`$ in Eq.(4), the critical temperature can be calculated by

$`T_c^{[k+p,k]}(n)`$ $`=`$ $`\left(K_1^{}\right)^1=\left[K_1^{(0)}+\delta K_1\right]^1`$ (5)
$``$ $`\left[K_1^{(0)}{\displaystyle \frac{\left(S_1^{(n+p)}_{L_{k+p}}S_1^{(n)}_{L_k}\right)|_{K_1^{(0)}}}{\left(\left[S_{L_{k+p}}\right]_{\mathrm{1\; 1}}^{(n+p,0)}\left[S_{L_k}\right]_{\mathrm{1\; 1}}^{(n,0)}\right)|_{K_1^{(0)}}}}\right]^1,`$
where the symbol $`T_c^{[k+p,k]}(n)`$ indicates that the temperature is obtained from two starting carpets of size $`3^{k+p}`$ and $`3^k`$ after $`n+p`$ and $`n`$ renormalization steps respectively. Fig. 2 shows the evolution of $`T_c^{[k+p,k]}(n)`$ calculated from pairs of the starting Sierpiński carpets at the simulation temperatures $`T_{sim}=\left[K_1^{(0)}\right]^1=1.4813`$ and $`T_{sim}=\left[K_1^{(0)}\right]^1=1.4795`$. These results deserve the following comments :

i) At the two simulation temperatures, $`T_c^{[k+p,k]}(n)`$ decreases as $`n`$ increases and tends to converge to some value.

ii) If we fix $`p=1`$, the value of $`T_c^{[k+p,k]}(n)`$ decreases as $`k`$ increases. We can see that $`T_c^{[5,4]}(n)>T_c^{[6,5]}(n)>T_c^{[7,6]}(n)>T_c^{[8,7]}(n)`$ and the curves $`T_c^{[7,6]}(n)`$ and $`T_c^{[8,7]}(n)`$ are very close to each other.

iii) $`T_c^{[k+p,k]}(n)`$ decreases as $`k`$ decreases if the size of the largest starting Sierpiński carpet $`3^{k+p}`$ is fixed. One can observe that $`T_c^{[8,7]}(n)>T_c^{[8,6]}(n)>T_c^{[8,5]}(n)>T_c^{[8,4]}(n)`$ where the size of the largest starting carpet is fixed at 6561. These curves tend to converge to some temperature while $`n`$ increases.

The best estimation for the critical temperature can be obtained from the 6-data points $`T_c^{[8,7]}(n)`$ for $`n=0,1,\mathrm{},5`$, by performing the following three-parameter $`(T_c,A_0,n_0)`$ fit:

$$T_c^{[k+p,k]}(n)=T_c+A_03^{n/n_0}.$$
(6)
The results of the fits at the two simulation temperatures lead to :

at $`T_{sim}=1.4813`$, $`(T_c,A_0,n_0)=(1.47959(15),\mathrm{\hspace{0.17em}0.02718}(28),\mathrm{\hspace{0.17em}0.839}(25))`$
with the reliability $`R^2=0.99968`$ and
at $`T_{sim}=1.4795`$, $`(T_c,A_0,n_0)=(1.47933(27),\mathrm{\hspace{0.17em}0.02827}(49),\mathrm{\hspace{0.17em}0.868}(42))`$
$`\text{with the reliability }R^2=0.99911.`$
The best value of the critical temperature can be estimated by taking the average of the results of the two fits; we hence have $`T_c=1.47946(16)`$.

In the following subsections, the eigen-exponents will be calculated from simulations carried out at $`T_{sim}=1.4795`$, since it lies closer to the best estimate of the critical temperature $`T_c=1.47946(16)`$. The simulations carried out at $`T_{sim}=1.4813`$ will provide a consistency check. It is well known that the exponents provided by the MCRG method are not very sensitive to small deviations from the critical temperature.

### 4.2 Temperature eigen-exponent $`y_t`$

We work in the even-coupling subspace and consider in turn the first 1 to 25 couplings shown in Fig. 1(a). We calculate the $`𝒯`$ matrix in the subspace at each renormalization level by inverting Eq.(3) with the help of the Gauss-Jordan elimination method ($`\mathrm{\S }2.1`$ in Ref.). The dimension of $`𝒯`$ is $`N_{coupling}\times N_{coupling}`$ where $`N_{coupling}`$ varies between 1 and 25. The largest eigenvalue $`\lambda _t`$ of $`𝒯`$ is obtained by applying the eigenvalue searching method for a real non-symmetric square matrix given in $`\mathrm{\S }11.5`$ and $`\mathrm{\S }11.6`$ of Ref.. The temperature eigen-exponent is then calculated by $`y_t=\mathrm{log}_3\lambda _t`$. Fig. 3 shows the evolution of the eigen-exponent $`y_t`$ as a function of the number of even couplings $`N_{coupling}`$, obtained from starting Sierpiński carpets of different sizes and at different renormalization levels. One can see that $`y_t`$ tends to converge towards stable values when $`N_{coupling}`$ increases, within the statistical errors, excepted at the highest renormalization levels where the sizes of the reduced lattices are $`3^3`$ and $`3^2`$. We estimate $`y_t`$ in the infinite couplings limit by taking the average of the values of $`y_t`$ obtained from 16 to 25 even couplings. The average value at each renormalization level is reported in Table 1 where we have disregarded the unreliable results at the highest renormalization levels. We find that the finite-size effect is not significant because at the same renormalization level, the values of $`y_t`$ obtained from the starting Sierpiński carpets of different sizes are about the same. A similar situation has been observed in the study of the 3D Ising model on regular lattices . In order to extrapolate the value of $`y_t`$ on the starting Sierpiński carpet of infinite size, we plot in the Fig. 4 the value of $`\lambda _t`$ obtained at level $`n(n+1)`$ from the starting carpet $`SC(3,1,n+4)`$ with respect to $`n`$. We find that, up to the renormalization levels which were performed, $`\lambda _t`$ seems not to be convergent. This behavior is quite different from the one observed in the case of the 3D Ising model on regular lattices  where $`\lambda _t`$ tends to converge to some value while $`n`$ increases. We, hence, cannot estimate the convergence by performing the three-parameter fit ($`\lambda _t^{},a_t,\omega _t`$) proposed in Ref.:

$$\lambda _t=\lambda _t^{}+a_t3^{\omega _tn}.$$
(7)
The lack of convergence in the case of the Sierpiński carpet may be due to a small value of the scaling correction exponent $`\omega _t`$ . At the present stage, we are only able to provide an upper bound for the Sierpiński carpet of infinite size: $`\lambda _t<1.781(10)`$ or $`y_t<0.5254(51)`$. The same situation is also observed at $`T_{sim}=1.4813`$. The difficulty in estimating the value of $`y_t`$ has already been observed in the study based on FSS analysis . According to Widom’s homogeneity hypothesis, $`y_t`$ is equal to the inverse of the correlation-length exponent $`\nu `$. Monceau et al.  have shown that the maxima values of the logarithmic derivatives $`\mathrm{\Phi }_i`$ (for $`i=1,2`$), defined by

$$\mathrm{\Phi }_i=E\frac{E|M|^i}{|M|^i}$$
(8)
where $`E`$ is the total energy and $`M`$ the total magnetic moment, are affected by scaling corrections. Instead of scaling properly as a power law, the maxima values of $`\mathrm{\Phi }_i`$ (for $`i=1,2`$) appear to exhibit a slight concavity in a log-log plot with respect to the lattice size. Monceau et al. could only provide an upper bound for $`y_t`$ from their MC study: $`1/\nu <1/1.5650.639`$. Our estimation is consistent with their result.

### 4.3 Magnetic eigen-exponent $`y_h`$

We now focus on the odd-coupling subspace. We consider in turn the first $`N_{coupling}`$ odd couplings shown in Fig.1(b) where $`N_{coupling}=1,2,\mathrm{},11`$. The $`𝒯`$ matrix in the odd-coupling subspace is obtained from Eq.(3) by considering only the lattice sums conjugated to the odd couplings. We denote $`\lambda _h`$ the largest eigenvalue of this matrix, and the magnetic eigen-exponent is calculated by $`y_h=\mathrm{log}_3\lambda _h`$. Fig. 5 shows the evolution of $`y_h`$ as a funtion of the number of odd couplings at each renormalization level obtained from the starting Sierpiński carpet $`SC(3,1,k)`$ with $`k=4,5,\mathrm{},8`$. We can see that $`y_h`$ tends to be stable as $`N_{coupling}`$ increases. Larger fluctuations are observed at the highest renormalization level where $`y_h`$ is extracted from the reduced carpets of sizes $`3^3`$ and $`3^2`$; we will disregard these values. We estimate $`y_h`$ in the infinite couplings limit by taking the average value of them when 7 to 11 odd couplings are considered. The results are reported in Table 2. The size effect of the starting Sierpiński carpet on the evolution of $`y_h`$ with respect to the renormalization level is more significant than it is on the evolution of $`y_t`$. However, the value of $`y_h`$ at level $`n(n+1)`$ from the starting carpet $`SC(3,1,n+4)`$ tends to converge to some value as $`n`$ increases (see Fig. 6). We hence perform the three-parameter fit suggested by Baillie et al. :

$$\lambda _h=\lambda _h^{}+a_h3^{\omega _hn}.$$
(9)
The fit yields $`(\lambda _h^{},a_h,\omega _h)=(7.38298(66),0.0737(12),\mathrm{\hspace{0.17em}1.86}(12))`$ with a reliability equal to $`R^2=0.99949`$. According to these results, the magnetic eigen-exponent $`y_h`$ takes the value $`1.81973(9)`$ and the associated scaling correction exponent $`\omega _h`$ the value $`1.86(12)`$. $`y_h`$ is slightly larger than the one obtained by Monceau et al. . They found that $`\gamma /\nu =1.732(4)`$ and $`\beta /\nu =0.075(10)`$. According to the relations: $`\gamma /\nu =2y_hd_f`$ and $`\beta /\nu =d_fy_h`$ , their results correspond to $`y_h=1.812(2)`$ and $`y_h=1.818(10)`$, respectively. No significant discrepancy can be brought out from these results, since the relative difference between the exponent provided respectively by the MCRG method and MC simulations remains smaller than $`0.5`$ percent.

## 5 Results from FSS analysis

The lattice sums associated with the first odd and the first even couplings represent respectively the total magnetic moment $`M`$ and minus the total energy $`E`$. We can, hence, perform a FSS analysis for a consistency check. According to this analysis, the thermodynamical average of the total magnetic moment at $`T_c`$ should scale as $`|M|L^{d_f\beta /\nu }=L^{y_h}`$ and the derivative of the Binder’s cumulant $`U=1M^4/(3M^2^2)`$ at the critical temperature $`T_c`$ scales as

$$\frac{dU}{d\beta _B}=(1U)(\mathrm{\Phi }_42\mathrm{\Phi }_2)L^{1/\nu }=L^{y_t},$$
(10)
where $`L`$ is the lattice size and $`\beta _B=(k_BT)^1`$ is the inverse of temperature. In Fig. 7, we plot the values of $`|M|`$ and $`dU/d\beta _B`$ at $`T=1.4795`$ for the 5 different lattice sizes from $`3^4`$ to $`3^8`$ in logarithmic coordinates. We find that the values of $`|M|`$ line up along straight lines, with the lattice size covering several order of magnitude. The slope provided by a least-square fit from the $`5`$ above points, is $`1.8198(11)`$, with a reliability $`R^2=0.99994`$. It turns out that this value is consistent with the value of $`y_h`$ obtained from MCRG method in the subsection 4.3. The behavior of $`dU/d\beta _B`$ is not affected by scaling corrections, (see Fig.7), although the $`\mathrm{\Phi }_i`$ (for $`i=1,2`$) exhibit a slight concavity as a function of the lattice size in a log-log plot. The reason why $`dU/d\beta _B`$ is not affected by scaling corrections may be a kind of “magic” cancellation of these scaling correction effects in the difference between $`\mathrm{\Phi }_4`$ and $`2\mathrm{\Phi }_2`$ . The slope measured from $`dU/d\beta _B`$ is $`0.449(6)`$ with the fitting reliability $`R^2=0.99948`$. This value is consistent with the result obtained in the subsection 4.2 where $`y_t<0.5254(51)`$.

## 6 Conclusion

The MCRG method has been shown to provide reliable values of the critical temperature $`T_c`$ and the magnetic eigen-exponent $`y_h`$, for the Ising model in a case where the underlying network has a fractal structure. The scaling correction exponent $`\omega _h`$ associated with $`y_h`$ can also be calculated. Moreover, difficulties are encountered in the estimation of the temperature eigen-exponent $`y_t`$ where, unlike regular lattices, the convergence of the renormalization flow in this direction seems to be very slow; nevertheless, an upper bound for $`y_t`$ can be provided. A similar slowness in the convergence of the logarithmic derivatives of the magnetization $`\mathrm{\Phi }_i`$ to the thermodynamical limit has already been brought out from MC simulations in the case of the same Sierpiński fractal  . In reference $`11`$, scaling corrections have been shown to affect strongly the behavior of the $`\mathrm{\Phi }_i`$’s and only an upper bound for $`y_t`$ has been calculated from their finite size behavior. The difficulties in calculating $`y_t`$ for fractal structures arise in the two methods, but emerge differently. At last, the results provided by the MCRG method have been shown to be consistent with the one obtained by Monceau et al.  and with an additional FSS analysis of the MCRG results.

## Acknowledgements

A part of the numerical simulations has been carried out in the Institut de Développement et des Ressources en Informatique Scientifique (IDRIS), supported by the Centre National de la Recherche Scientifique (project number $`021186`$). We acknowledge the scientific committee and the staff of the center. We are also grateful to the Centre de Calcul Recherche (CCR) of the University Paris VII–Denis Diderot, where the rest of the simulations have been done.

## Figure captions

1. (a) 25 even couplings and (b) 11 odd couplings, considered in our MCRG study. They are ordered according to the Blöte’s importance factor $`F`$ whose value is given in the right colomn. The symmetry number is given in the next to last right colomn.
2. Evolution of $`T_c^{[k+p,k]}(n)`$ with respect to the different pairs of the starting Sierpiński carpets $`SC(3,1,k+p)`$ and $`SC(3,1,k)`$ (denoted by $`[k+p,k]`$) obtained at (a) $`T_{sim}=1.4813`$, (b) $`T_{sim}=1.4795`$.
3. Variation of $`y_t`$ with respect to the number of even couplings $`N_{coupling}`$ at different renormalization levels from the starting Sierpiński carpet $`SC(3,1,k)`$ with $`k=4,5,\mathrm{},8`$.
4. $`\lambda _t`$ obtained at level $`n(n+1)`$ from the starting carpet $`SC(3,1,n+4)`$.
5. Variation of $`y_h`$ with respect to the number of odd couplings $`N_{coupling}`$ at different renormalization levels from the starting Sierpiński carpet $`SC(3,1,k)`$ with $`k=4,5,\mathrm{},8`$.
6. $`\lambda _h`$ obtained at level $`n(n+1)`$ from the starting carpet $`SC(3,1,n+4)`$.
7. $`|M|`$ and $`dU/d\beta _B`$ at $`T=1.4795`$ on SC(3,1,k) with $`k=4,5,\mathrm{},8`$ in logarithmic coordinates.


text:
# Optimization as a result of the interplay between dynamics and structure

## 1 Introduction

In recent times the possibility of using the tools of statistical physics to analyze the rich dynamical behaviors observed in social, technological, biological and economical systems has attracted a lot of attention from the physics community. So far, one of the main contributions to these fields has been the analysis of simple models that capture the basic features of the investigated phenomena. The goal is to identify their relevant ingredients as well as the essential mechanisms governing their dynamics with the hope that this information will help us to understand the physical behavior of real complex systems. A great part of this effort has been devoted to the characterization of real networks, identifying their main features, and understanding how they arise since in nature the coupling and interaction between units always has a clear motivation . This approach provides a valuable tool to analyze how the underlying network affects the dynamics of interest. This is a natural step forward along this line of investigation.

For example, studies of an Ising model defined on a small world network  showed that for any finite disorder strength the system undergoes a crossover to a mean-field like region for small temperatures in the thermodynamic limit. Also, the analysis of a biological evolution model defined on a small world  has shown that even for a small number of quenched long-range jumps in the genotypic space the results are indistinguishable from those obtained by assuming all mutations equiprobable. The dramatic change in some static quantities such as the mean distance between units has been argued as the key mechanism to explain these phenomena. However, there is another set of works pointing out that other elements are necessary in order to observe changes when the underlying structure is modified, something which cannot be directly related to known static properties of the network. Along this line, numerical simulations on a non-equilibrium directed Ising model on a small world  showed a non trivial phase diagram. This system exhibits a line of continuous phase transitions below a critical threshold of disorder, above which the transition becomes first order. In agreement to this model other systems display transitions between different regimes but only when a finite number of shortcuts is introduced. In a model for the spreading of rumors  this critical value separates two regimes. In the first regime the rumor is bound to a finite neighborhood of an initially infected site while on the other one it reaches a finite fraction of the population. A similar phenomenon was observed in a model for the spreading of diseases  where a transition from an endemic situation to an oscillatory one was found.

These are just a few modern examples where significant changes in the properties of models are observed when their underlying structure is modified. Our work fits in this intriguing area. However, we will not focus our attention only on what happens when a structure is replaced by another. There are systems whose behavior can be described in terms of macroscopic observables, which give an idea about how well a system performs in a given scenario. Precisely, in this work the interplay between the dynamics and the underlying network is going to be exploited in order to reach an optimal behavior. Keeping this idea in mind we will show that a clear separation between both ingredients is not possible. In this way given a particular set of features which characterize the system, we wonder for the binomial dynamics/structure that optimizes its performance, measured in terms of the appropriate magnitudes. As we will see the answer is not trivial, even for the simple model of diffusion of innovations  that will be under consideration in this paper.

The outline of this work is as follows. First we will define the model, and describe its main features. Then we will analyze the influence of different underlying networks on the dynamical properties of the system. Using this information, we will study the interplay between structure and dynamics. Finally, we will show that both elements need to be taken into account in order to achieve an optimal regime.

## 2 The model

Our starting point is a model of diffusion of innovations by imitation recently proposed in . In spite of its apparent simplicity it displays a very rich behavior presenting self-organization, subcritical, critical and supercritical regimes and deserves further analysis. The model was originally presented in a computational economics context , describes how an external perturbation that affects individuals propagates to the rest of the population by a mechanism of imitation. Such a perturbation, modeled in terms of random updates, is independent of the underlying network and may affect any agent in the system. When an agent changes her state such information will be available only to the agents which are in contact with her. These agents should decide if it is beneficial for them to imitate, or remain in their present state. In the case that any of these agents decides to imitate the information will also propagate to her neighbors. The process ends when everybody is satisfied with their current situation. In this way, the features that characterize the members of the population are diffused as a wave of changes throughout the system, eventually reaching all agents. Clearly, the main features of the whole process will depend on how we model the mechanism of imitation as well as the way in which the agents are connected.

To be more precise, let us explicitly enumerate the dynamical rules governing the behavior of the individuals. Each site (or agent) $`i`$ is characterized by a real variable $`a_i`$. In a general way we can consider this quantity as a characteristic of a given individual that other agents might want to imitate. When an agent has adopted a new characteristic, her neighbors become aware of the change and balance their interest (quantified as $`a_ia_j`$) with their resistance to change $`C`$ to decide if they would like to imitate this change. In this way $`C`$ controls the mechanism of imitation. This parameter is constant and the same for all the agents in the system.

The dynamics can be summarized as follows :

1. The system is asynchronously updated, that is, at each time step a randomly selected agent $`a_i`$ updates her state

$$a_ia_i+\mathrm{\Delta }_i,$$
(1)
where $`\mathrm{\Delta }_i`$ is a random variable exponentially distributed with mean $`\lambda `$. This driving process accounts for the external pressure that may lead an individual to spontaneously adopt a new characteristic.
2. All agents $`j\mathrm{\Gamma }(i)`$, where $`\mathrm{\Gamma }(i)`$ is the set of neighbors of agent $`i`$, decide whether they want to upgrade or not, according to the following rule

$$a_ia_jCa_j=a_i,$$
(2)
3. If any $`j\mathrm{\Gamma }(i)`$ has decided to imitate, we also let her neighbors decide whether they want to imitate this behavior or not. In this way the information of an update may spread beyond the first neighbors of the originally perturbated site. This procedure is repeated until no one else wants to change, concluding an avalanche of imitation events. In this way we have assumed that the time scale of the imitation process is much shorter than the one corresponding to the external driving.

Regardless of the particular connectivity pattern between agents, there are some common trends that helps us to understand intuitively the collective behavior of the system. Hence, it is easy to understand the role of the parameter $`C`$ and see how different dynamical regimes come out in a natural way. For small $`C`$ there is almost no resistance to change, and the information spreads easily through the system. Avalanches involve a large number of units and therefore agents are continually updating their state, although by small amounts. On the other hand, for large $`C`$, since the resistance to change is high, the number of units involved in the diffusion process is rather small which implies a localized region of propagation for most of the updates. As a consequence, distant agents can present very different characteristics, and only in few occasions they are able to find out about their differences. When this happens the changes (variation of $`a_i`$) adopted will be large. For intermediate values of $`C`$ the behavior is rather complex and the most interesting from a collective standpoint.

A desirable situation would be one in which most of the time the agents adopt changes that lead to large advances. In this context the following question arises: what is the resistance to change that permits the agents to reach a given average level with a minimum number of upgrades? By measuring the mean rate of advance $`\rho `$ it has been shown that there exists a unique value of $`C`$ for which this optimal regime can be achieved . This quantity is simply defined as the ratio between the total advancement over the total number of upgrades, that is

$$\rho =\underset{T\mathrm{}}{lim}\rho (T)=\underset{T\mathrm{}}{lim}\frac{_{t=0}^TA(t)}{_{t=0}^Ts(t)}$$
(3)
where $`s(t)`$ is the number of sites that have changed (i.e. the avalanche size) at time $`t`$, and $`A(t)=_i\left(a_i(t)a_i(t1)\right)`$ is the advance triggered by the avalanche. The optimal growth regime is characterized by the presence of a peak, $`\rho _{max}`$, in the mean rate of advance for a given value of the parameter $`C`$. This peak scales with system size as a power law $`\rho _{max}N^\alpha `$. This property is very convenient in certain frameworks. For example, from an economic point of view the model displays scale effects, i.e., large economies grow faster.

So far, we have considered general aspects of the model. Since the main goal of this paper is to analyze the interplay between underlying structure and dynamics, let us put some attention in the influence of the connectivity pattern on the collective properties of the system and, in particular, in macroscopic observables such as $`\rho `$. There are two simple cases which are systematically analyzed: a fully connected network and a regular lattice. In the globally coupled case the information referent to any change elicited in an arbitrary position of the network is immediately available to every other agent. As a consequence, the state of all the agents $`a_i`$ is bound in a gap of width $`C`$. This limits the advance of any agent to a maximum of $`C`$. On the other hand, when the system is defined on a 1D ring this limitation only applies to the nearest neighbors. If the information spreads beyond these neighbors, the advance achieved can exceed $`C`$. In this extreme cases different qualitative behaviors are observed, which can be quantified by the exponent $`\alpha `$. When one considers the dynamics of the model on a ring $`\alpha =0.20(2)`$ while mean field calculations and numerical simulations of the dynamics of the model on a fully connected network show that $`\alpha =0.50`$ .

Summarizing, we have presented a model that presents an optimal regime which can be characterized by the mean rate of advance $`\rho `$. The scaling properties of this magnitude can be used to quantitatively distinguish among different dynamical behaviors. These dynamical behaviors have been observed in two particular underlying networks. In the next section we will consider more general structures and analyze their effects on the dynamical properties of the model.

## 3 Scaling analysis

The two particular cases considered so far are usually chosen either for their numerical simplicity, as in the ring, or because they allow for simple mean field calculations as in the fully connected network. Clearly, the structure of both cases is far from the much more complex pattern of interactions observed in realistic systems. However, they appear as paradigms of two opposite generic situations either a scenario where the propagation of information can be constrained to a local neighborhood, or one where it may reach the whole population in just one step. One wonders which will be the dominant features in a more general situation. In this section we will consider this issue. In particular, we will study the relation between the salient features present in a more general structure and the dynamical properties of the system.

It is well known that, starting from a ring, the random addition of a few links produces changes in the properties of the network, such as a rapid drop in the average distance between nodes, maintaining the local structure . It has been shown that this feature can be related to significant changes in the dynamical properties in some systems . In order to study what effects may be present in our model, we will consider the dynamics on a ring lattice with $`N`$ vertices and $`k=2`$ edges per vertex, adding a new link at random with probability $`p`$ per edge . In general, when one modifies the underlying structure, a quantitative variation of the numerical values of the dynamical magnitudes that characterize the system is observed. However, we will focus our attention on whether the scaling properties are modified, since they are an indication of qualitative changes in these properties. In order to quantify the changes in the scaling behavior of the system we have computed the exponent $`\alpha `$ for different values of $`p`$, as shown in Fig. 1.

The main conclusion that can be extracted from the figure is that the addition of new links does not modify substantially the scaling properties of the system. For small $`p`$ we observe a small decrease in the value of $`\alpha `$, effect that seems to be correlated with the rapid drop in the average distance between agents. As $`p`$ increases to $`p=1`$ a slight growth of $`\alpha `$ is reported. It is important to understand why the qualitative behavior of the system is similar to the one characteristic of the ring. In this case, an avalanche propagates through steps in which the information reaches the nearest neighbors of a modified site. To generate a large event or simply to reach agents far away from the initial updated unit a large number of steps is required. In this sense, the information propagates by a local process. The addition of new links allows the information to reach agents through short cuts, but it does not change the mean mechanism of diffusion, i.e. in order to proceed further, an avalanche still requires a large number of steps, it is still dominated by a local process. One cannot under-stress the fact that for the imitation strategy described in this paper, in contrast to what is observed in other models, the characteristic behavior of the ring is dominant even when a more general structure, such as a small world network, is considered.

The local character of the diffusion process typical of low connectivity networks is lost when considering densely connected networks. In this situation the information about the state of a given agent is available to any other agent in just a few steps. The small world construction is far away from this limit, even when $`p=1`$ the mean number of links per site, i. e. the connectivity of the system, has increased very slightly. In order to reach a highly connected state we added links randomly to the ring up to a fraction $`f`$ of the total possible links in the system. This means adding $`f[N(N1)/2]N`$ connections to the ring. This recipe allows us to interpolate between the ring ($`f=0`$) and the fully connected model ($`f=1`$). For finite values of $`f`$ the system corresponds to a ring with a superimposed random graph. For this construction the connectivity of the system will be $`fN`$, and thus, for a fixed value of $`f`$, the connectivity will increase as $`N`$ grows. If this plays a significant role in the dynamical behavior of the system, we expect that it will be reflected in an important quantitative change in $`\alpha `$. In Fig. 2 we present the behavior of the peak of the mean rate of advance $`\rho _{max}`$ as a function of system size $`N`$ when $`f=0.005`$. For low values of $`N`$ the behavior resembles the one observed in the small world case. As $`N`$ grows there is a crossover to the fully connected behavior and $`\rho _{max}N^{0.5}`$.

To complete the analysis we have performed a finite size study of the crossover between both regimes. The results are presented in Fig. 3. As one can expect, and in agreement with the previous discussion, the jump from the ring-like scaling behavior to the highly connected network scaling behavior becomes sharper as N grows being the crossover located at $`f\frac{1}{N}`$. Numerical simulations support this result suggesting that in the thermodynamic limit a transition to $`\alpha =0.5`$ should be present for any non-vanishing value of $`f`$, i.e. provided the connectivity of the system grows with the system size. Therefore, by using this construction we have found that there are two dominant behaviors, one for networks with a number of links scaling with the size of the system, $`O(N)`$, and another where the connectivity grows with the size of the system, i.e., the number of links go as $`O(N^2)`$. It will be very interesting to analyze the situation for other complex situations such as scale free networks. This is currently under study and will be published elsewhere.

We have remarked the fact that the behaviors observed in the ring and the fully connected network were generic regarding the way the information propagates through the system. We found that these dynamical behaviors are dominant in more complex structures. In a small world network the dynamics is dominated by local processes. In order to reach the qualitative behavior observed in a fully connected network a higher connectivity is necessary. In fact, we have shown that by introducing a mechanism which takes this effect into account we were able to reach this qualitative behavior. The identification of the dominant behaviors of our model in more general networks is an interesting issue in itself, but it is also an essential step in exploiting the interplay between dynamics and structure to reach an optimal regime. Precisely, to consider both elements in a unified framework is the main purpose of this work. The following section is devoted to the investigation of our model under this new exciting point of view.

## 4 Identification of an optimal behavior

Up to now, we have analyzed different mechanisms for the spreading of the information. On one hand, $`C`$ appears as a parameter that can be tuned in order to reach an optimal regime for a fixed underlying network. On the other hand, we have also seen that diverse structures lead to different ways in which the information spreads. This line of reasoning naturally leads to the question of which is the structure that gives the optimal regime for a fixed value of $`C`$. For small $`C`$ there is almost no resistance to change and the information easily spreads. In fact we expect that if $`C0`$ then $`\rho `$ will also decrease independently of the underlying network. On the other extreme, for $`C\mathrm{}`$ the behavior of the system will resemble a random deposition process, and again the behavior of $`\rho `$ is expected to decrease when considering any general structure. To analyze what happens for intermediate values of $`C`$ we should proceed with care as we will see immediately.

In order to do this we consider a system with a certain distribution of couplings between agents and a fixed $`C`$. In general, the value of $`\rho `$ will be smaller than $`\rho _{max}`$. Now, suppose that we can modify the connectivity pattern in order to reach the optimal behavior. What should we do? Is it better to increase the number of links and tend to a densely connected network or, on the contrary, we should go towards a regular structure with a small number of connections per site? The question cannot be answered properly without looking at the parameter that control the dynamics, i.e., the strategy to follow depends precisely on $`C`$. To optimize the behavior of the system one cannot split the problem in two independent parts; dynamics and the underlying structure must be considered as a whole.

Let us analyze some features associated to the mechanism of imitation and consider the difference between the highest ($`a_{max}`$) and the lowest ($`a_{min}`$) characteristic values in the system. For sufficiently high connectivity this gap will be limited to $`a_{max}a_{min}C`$ and, as a consequence, for these systems the value of $`\rho `$ cannot exceed $`C`$. On the other hand, when the propagation of the information is constrained to advance through local processes, a more heterogeneous profile of characteristic values can be formed, allowing for a larger gap between $`a_{max}`$ and $`a_{min}`$. In this case, avalanche events consisting of a large number of steps will produce a large advance, allowing the value of $`\rho `$ to exceed $`C`$. In fact, by using the probability distribution of avalanches $`P(s)`$ and advances $`P(H)`$ obtained using numerical simulations in  this can be easily verified analitically. In this situation a sparsely connected network will necessarily have a greater $`\rho `$ than a highly connected one. For increasing values of $`C`$, avalanche events will easily get blocked in a few steps. In this context, large advances in a sparsely connected network will become very rare. More frequent advances will be observed in a highly connected network since many agents are permitted to find out about an update in any step. This situation offers the possibility for a higher $`\rho `$ to be observed in highly connected structures.

Following this analysis we have considered a general system of fixed size $`N`$. The evolution of the mean rate of advance $`\rho `$ vs $`C`$ has been studied for two different situations: the ring and another structure where the number of links (measured in terms of $`f`$) is large enough to observe fully connected behavior. The results are illustrated in Fig. 4. Note that as $`f`$ is varied the qualitative shape of the $`\rho `$ curve is similar. However, the position of the peak corresponds to different values of $`C`$. For increasing values of $`f`$ the peak corresponds to larger values of $`C`$, eventually reaching the curve corresponding to the globally coupled case. It is important to stress that, as $`C`$ is varied, the optimal network may change from a highly connected to a sparsely connected one. This behavior is clearly reflected in the inset, where we present the behavior of $`\rho `$ vs. $`f`$ for two different values of $`C`$. For $`C=2`$ a decrease in $`\rho `$ is observed as $`f`$ grows. The addition of links is harmful for the system. On the other hand, for $`C=4.0`$ the opposite behavior is observed, and $`\rho `$ increases its value as $`f`$ grows. Clearly, in this case, the addition of links is beneficial.

These results show that in order to optimize the behavior of the system a non trivial combination of dynamical rules and underlying structure should be considered. When the interplay between both allows for $`\rho /C>1`$, a sparsely connected structure performs better than a highly connected one. On the other hand, when $`\rho /C<1`$ the opposite is true. Note that when $`\rho /C1`$ the behavior should be independent of the underlying network. In fact, Fig 4 shows that when $`\rho /C1`$ both curves intersect. Numerical simulations for different $`N`$ and $`f`$ also show these qualitative behaviors.

## 5 Conclusions

The search for the essential mechanisms that govern the dynamics of real complex systems has led recent efforts to highlight the role of the complex underlying structures present in these systems. A large amount of work has been devoted to the characterization and identification of their main features and the role they play on the dynamical properties. This scenario suggests that a more general framework, devised to address these issues, should consider the interplay between the dynamics and the underlying structure as a single ingredient.

The main aim of this work has been to take a first step into this direction. With this picture in mind we analyzed a model of diffusion by imitation which presents a very interesting property, the possibility of defining an optimal regime by tuning a parameter. This regime can be quantified by a macroscopic observable. In order to be able to consider the interplay between dynamics and structure we have studied the behavior of the model on general networks and characterized the dominant behaviors in these structures. An analysis of how the optimal regime can be achieved has led to the conclusion that the dynamical rules and the underlying structure cannot be considered as separate ingredients.

The rich dynamical behaviors observed in the more general framework makes a detailed analysis of the model necessary. In particular, we are focusing our attention on the behavior when $`f0`$, in order to characterize a possible transition. Since the connectivity plays an important role in the properties of the system, a natural extension to more general structures, such as scale-free networks, is under consideration. Finally, a very interesting new point of view which we are also analyzing concerns the behavior of the model on dynamical networks. Although we cannot underestimate the results that will follow from these studies we believe that the main results presented in this work will not be affected in their generality.

## 6 Acknowledgements

We acknowledge financial support from MCYT, grant number BFM2000-0626, and also from European Commission, Fet Open Project COSIN IST-2001-33555. P. M. Gleiser acknowledges financial support from Fundación Antorchas.