text:
# A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks

## 1 INTRODUCTION

Dwarf novae are close binary systems that consist of a white dwarf and a Roche-lobe-filling secondary star. The matter overflowing from the secondary star forms an accretion disk around the white dwarf. Dwarf nova systems show repetitive outbursts in which the luminosity of the disks increase rapidly. The duration of the outbursts ranges from a few to 20 days, and the recurrence time is about 20 – 300 days. The disk instability model (Osaki 1974), in which the evolution of the disk is regulated by changes in the rate of angular momentum transport, is generally accepted as the explanation of the outbursts (Cannizzo 1993; and references therein). Since ordinary molecular viscosity is too low to explain the evolutionary timescale of the disks, an anomalous stress associated with turbulence and characterized by the $`\alpha `$ parameter (Shakura & Sunyaev 1973), is thought to be the source of angular momentum transport. Detailed comparisons between observed light curves of dwarf novae and the theoretical disk models suggest that the viscous parameter $`\alpha `$ must vary between the hot (outburst) and cold (quiescent) state, and the amplitude of $`\alpha _{\mathrm{hot}}`$ and $`\alpha _{\mathrm{cold}}`$ is typically of the order of 0.1 and 0.01, respectively. However, while the origin of the anomalous stress has become better understood in recent years (see below), the reasons for the difference between $`\alpha _{\mathrm{hot}}`$ and $`\alpha _{\mathrm{cold}}`$ are still unclear.

The magnetorotational instability (MRI; Balbus & Hawley 1991) is the most promising source of the anomalous stress. MHD turbulence driven by the MRI can transport angular momentum by the Maxwell (magnetic) stress. During the hot state, the disk gas is fully ionized so that the ideal MHD approximation is appropriate. Numerical simulations that adopt ideal MHD have shown that the Maxwell stress caused by the MRI can give $`\alpha 0.1`$ to 0.01, where the saturation level is determined by the geometry and strength of the magnetic field (e.g., Hawley, Gammie, & Balbus 1995; 1996). Therefore, MHD turbulence can be the primary mechanism of angular momentum transport at least during outbursts, that is it can account for $`\alpha _{\mathrm{hot}}`$.

During the cold state, on the other hand, the gas is only weakly ionized. Gammie & Menou (1998) have pointed out that ohmic dissipation could modify the nature of the MHD turbulence at quiescence, and this could make the difference between $`\alpha _{\mathrm{hot}}`$ and $`\alpha _{\mathrm{cold}}`$. Calculations of the ionization fraction in dwarf nova disks reveal that the Hall effect as well as ohmic dissipation must be considered if the temperature is $`T2000`$ K (Sano & Stone 2002a). At low temperatures, ohmic dissipation can suppress the growth of the MRI (Jin 1996; Sano & Miyama 1999). If turbulence dies away completely as a result of ohmic dissipation, the Maxwell stress cannot provide the required anomalous stress during quiescence (Menou 2000). Thus, it is important to calculate numerically the amplitude of the Maxwell stress at the cold state.

In this paper, we examine the behavior of MHD turbulence during the transition from an outburst to quiescence using local three-dimensional MHD simulations. These simulations include the dominant non-ideal MHD effects as determined from a self-consistent calculation of the ionization state in the disk. Since we adopt a local approximation to study the turbulent stresses (usually referred to as a one-zone model, e.g., Mineshige & Osaki 1983; Cannizzo & Wheeler 1984), our calculations do not include global effects (such as spiral shocks in the disk, Sawada, Matsuda, & Hachisu 1986) which may contribute to angular momentum transport.

The plan of this paper is as follows. In § 2, we calculate the ionization fraction in dwarf nova disks by solving the Saha equation to estimate the magnitude of nonideal MHD effects at quiescence. Our numerical method is described in § 3. The results of numerical simulations which examine the dependence of the stress on the geometry of magnetic field and the decay timescale of turbulence are presented in § 4. In § 5, we discuss the activity of MHD turbulence and the source of angular momentum transport during quiescence.

## 2 IONIZATION STATE OF THE DISKS

During outburst, the temperature of the disk is above $`10^4`$ K, so that the gas is fully ionized and the ideal MHD approximation holds. In this case, the disk will be dominated by MHD turbulence driven by the MRI, with the plasma beta $`\beta 10^2`$$`10^4`$ (e.g., Hawley, Gammie, & Balbus 1995; 1996), so that the Alfvén speed is typically $`v_\mathrm{A}10^4`$$`10^5`$ cm s-1. During the transition to quiescence, the temperature drops to a few thousands K. The ionization fraction decreases dramatically, and thus nonideal MHD effects must be included. The induction equation is then

$$\frac{𝑩}{t}=\mathbf{}\times \left[𝒗\times 𝑩\frac{4\pi \eta 𝑱}{c}\frac{𝑱\times 𝑩}{en_e}+\frac{\left(𝑱\times 𝑩\right)\times 𝑩}{c\gamma _i\rho _i\rho }\right],$$
(1)
where $`𝒗`$ is the neutral velocity, $`\eta `$ is the magnetic diffusivity, $`𝑱`$ is the current density, $`c`$ is the speed of light, $`e`$ is the electron charge, $`n_e`$ is the number density of electrons, $`\gamma _i`$ is the drag coefficient, and $`\rho _i`$ and $`\rho `$ are the mass densities of ions and neutrals, respectively. Each term in the brackets of equation (1) from left to right means the inductive term ($`I`$), ohmic dissipation ($`O`$), Hall effect ($`H`$), and ambipolar diffusion ($`A`$).

The importance of the nonideal effects can be measured by the relative ratio of each term in equation (1) to the inductive term. These ratios are determined mainly by the ionization fraction. At a typical density of dwarf nova disks ($`n_n10^{18}`$ cm-3, where $`n_n`$ is the number density of neutrals), ohmic dissipation and Hall effect are very important and affect the evolution of the MRI (Sano & Stone 2002a). The size of the ohmic dissipation is given by the magnetic Reynolds number;

$$Re_M=\frac{v_\mathrm{A}^2}{\eta \mathrm{\Omega }},$$
(2)
where $`\mathrm{\Omega }=(GM/r^3)^{1/2}`$ is the angular velocity, $`G`$ is the gravitational constant, $`M`$ is the mass of the central object, and $`r`$ is the distance from the central star. The magnetic diffusivity is given by

$$\eta =\frac{c^2m_e\nu _e}{4\pi e^2n_e}=234\left(\frac{n_e}{n_n}\right)^1T^{1/2}\mathrm{cm}^2\mathrm{s}^1,$$
(3)
where $`\nu _e=n_n\sigma v_e`$ is an effective collision frequency for electrons. Here we take $`\sigma v_e=8.28\times 10^{10}T^{1/2}`$ (Draine, Roberge, & Dalgarno 1983). The magnetic Reynolds number $`Re_M`$ is the ratio of the inductive term to the ohmic dissipation ($`Re_M=I/O`$; Sano & Stone 2002a), and a smaller $`Re_M`$ represents a more resistive case. The Hall effect on the MRI is characterized by the Hall parameter (Balbus & Terquem 2001);

$$X=\frac{cB\mathrm{\Omega }}{2\pi en_ev_\mathrm{A}^2},$$
(4)
which is proportional to the ratio between the Hall and the inductive term ($`X=2H/I`$; Sano & Stone 2002a). If the Hall parameter $`X`$ is large, the Hall term is important.

Both the magnetic Reynolds number and Hall parameter can be written as a function of the number density of electrons and the temperature;

$$Re_M=68\left(\frac{v_\mathrm{A}}{10^5\mathrm{cm}\mathrm{s}^1}\right)^2\left(\frac{n_e/n_n}{10^6}\right)\left(\frac{T}{3000\mathrm{K}}\right)^{1/2}\left(\frac{M}{M_{}}\right)^{1/2}\left(\frac{r}{10^{10}\mathrm{cm}}\right)^{3/2},$$
(5)
and

$$X=0.0059\left(\frac{v_\mathrm{A}}{10^5\mathrm{cm}\mathrm{s}^1}\right)^1\left(\frac{n_e/n_n}{10^6}\right)^1\left(\frac{n_n}{10^{18}\mathrm{cm}^3}\right)^{1/2}\left(\frac{M}{M_{}}\right)^{1/2}\left(\frac{r}{10^{10}\mathrm{cm}}\right)^{3/2}.$$
(6)
Interestingly, the product $`Re_MX`$ is independent of the ionization fraction;

$$Re_MX=0.40\left(\frac{v_\mathrm{A}}{10^5\mathrm{cm}\mathrm{s}^1}\right)\left(\frac{n_n}{10^{18}\mathrm{cm}^3}\right)^{1/2}\left(\frac{T}{3000\mathrm{K}}\right)^{1/2}.$$
(7)
The number density of electrons is a function of the temperature, because the free electrons in dwarf nova disks are mainly generated by thermal ionization of Al and K. Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks shows the temperature dependence of the electron abundance $`n_e/n_n`$. We have solved the Saha equation assuming Solar abundance for the case of $`n_n=10^{18}`$ cm-3, which is a typical density of dwarf nova disks. The ionization fraction decreases as the temperature decreases. As seen from the figure, the fraction is very sensitive to the temperature when $`T`$ 2000 K. For example, the abundance is $`n_e/n_n3.5\times 10^8`$ at $`T=2000`$ K, but it is reduced by 2 orders of magnitude ($`n_e/n_n3.7\times 10^{10}`$) at $`T=1500`$ K.

Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks also shows the temperature dependence of the parameters $`Re_M`$ and $`X`$, which are given by equations (5) and (6), respectively. Here we assume $`v_\mathrm{A}=10^5`$ cm s-1, $`M=M_{}`$, and $`r=10^{10}`$ cm. When $`Re_M1`$ or $`X1`$, these nonideal MHD terms can modify the linear growth and nonlinear evolution of the MRI111The definition of the magnetic Reynolds number given by equation (2) is different from that used in the previous work ($`Re_M^{}=c_sH/\eta `$; Gammie & Menou 1998; Menou 2000), and thus the critical value is also different (see Sano & Stone 2002a,b). (Sano, Inutsuka, & Miyama 1998; Wardle 1999; Fleming, Stone, & Hawley 2000; Balbus & Terquem 2000; Sano & Stone 2002a,b). Then, nonideal MHD effects are clearly inefficient if the temperature is higher than 3000 K ($`Re_M1`$ and $`X1`$), so that ideal MHD may be valid. However, both the ohmic dissipation and Hall effect become important when $`T`$ 2000 K ($`Re_M1`$ and $`X1`$). Therefore, these nonideal MHD effects must be included in the analysis when the temperature of the disk becomes below 2000 K. At this range of the temperature, the electron abundance and number density is $`n_e/n_n10^8`$ and $`n_e10^{12}\mathrm{cm}^3`$, respectively.

For a typical condition of dwarf nova disks, the relative ratio of the ambipolar diffusion to the inductive term is given by

$$\frac{A}{I}=\frac{\mathrm{\Omega }}{\gamma _i\rho _i}=6.1\times 10^6\left(\frac{n_i}{10^{12}\mathrm{cm}^3}\right)^1\left(\frac{M}{M_{}}\right)^{1/2}\left(\frac{r}{10^{10}\mathrm{cm}}\right)^{3/2},$$
(8)
where we assume that $`\gamma _i=\sigma v_i/(m_i+m_n)`$, $`\sigma v_i=1.9\times 10^9`$ cm3 (Draine et al. 1983), the ion and neutral particle mass are $`m_i=30m_\mathrm{H}`$ and $`m_n=1.27m_\mathrm{H}`$, and $`m_\mathrm{H}`$ is the mass of a hydrogen atom. As mentioned in Gammie and Menou (1998), this ratio given by equation (8) is so small that the ambipolar diffusion does not affect the evolution of the MRI (Blaes & Balbus 1994; Brandenburg et al. 1995; Hawley & Stone 1998). Therefore, we ignore the ambipolar diffusion in the numerical simulations in this paper.

## 3 NUMERICAL METHOD

We study the evolution of the MRI and MHD turbulence during the transition to quiescence using three-dimensional MHD simulations. The shearing box approximation (Hawley et al. 1995) is adopted to investigate the local evolution of the instability. The equations of adiabatic MHD with the ratio of specific heats $`\gamma =5/3`$ are solved in a Cartesian frame of reference ($`x`$, $`y`$, $`z`$) corotating with the disk at the angular velocity $`\mathrm{\Omega }`$, where $`x`$ is the radial, $`y`$ is the azimuthal, and $`z`$ is the vertical direction. The induction equation includes terms that describe ohmic dissipation and the Hall effect. The details of the equations and our numerical scheme are described in Sano & Stone (2002b). Vertical gravity is ignored, so that the calculations begin with a uniform density, $`\rho =\rho _0`$, and a uniform pressure, $`P=P_0`$. The magnetic field structure of dwarf nova disks is still highly uncertain. If the white dwarf has a strong magnetic field, then the dipole field may penetrate the disk vertically. On the other hand, a toroidal field amplified by the differential rotation may be dominant in the disk. In this case, the net flux of the vertical field could be negligible. Nonlinear evolution of the MRI started with a toroidal field is quite similar to that with a zero net flux $`B_z`$ (Sano & Stone 2002b). Therefore, in this paper, we consider two kinds of initial field geometries; a uniform vertical field, $`B_z=B_0`$, and a zero net flux vertical field, $`B_z=B_0\mathrm{sin}(2\pi x/H)`$, where $`H=(2/\gamma )^{1/2}c_{s0}/\mathrm{\Omega }`$ is the scale height of the disk and $`c_{s0}`$ is the initial sound speed. We choose normalizations with $`\rho _0=1`$, $`H=1`$, and $`\mathrm{\Omega }=10^3`$. The magnetic field strength is characterized by the plasma beta $`\beta _0=P_0/(B_0^2/8\pi )=3200`$. This is about 100 G in dwarf nova disks, which corresponds to a dipole field of a white dwarf with a moderately strong magnetic field ($`10^6`$ G) at a radius $`r=10^{10}`$ cm. The size of the simulation box is $`H\times 4H\times H`$ for all the models calculated in this paper. The initial ratio of the MRI wavelength, $`2\pi v_{\mathrm{A0}}/\mathrm{\Omega }`$, to the box height is about 0.11. A standard grid resolution of $`32\times 128\times 32`$ uniform zones is used for most of our calculations.

To set up a disk in outburst, the induction equation is integrated for 50 orbits assuming ideal MHD. The simulations are then restarted with nonideal MHD effects included in the induction equation, and the models are evolved for a further 50 orbits. Although the transition should proceed at the thermal timescale of the disk, our technique assumes an instantaneous change of the temperature at 50 orbits, which is thereafter kept constant. The quiescent disk is characterized by the non-dimensional parameters $`Re_M`$ and $`X`$. We have computed 4 different models with $`(Re_M,X)=(100,0.01)`$, (10,0.1), (1,1), and (0.1,10), which correspond to the temperature $`T=2800`$, 2200, 1800, and 1600 K, respectively (see Fig. 1).

As seen from equations (5) and (6), the parameters $`Re_M`$ and $`X`$ are functions of the magnetic field strength. Because the vertical component of the field is essential for the evolution of the MRI (Sano & Stone 2002b), we define these parameters ($`Re_{M,50}`$ and $`X_{50}`$) using the volume-averaged vertical field strength $`B_z^2^{1/2}`$ and squared Alfvén speed $`v_{\mathrm{A}z}^2`$ at 50 orbits, i.e., $`Re_{M,50}=v_{\mathrm{A}z}^2/\eta \mathrm{\Omega }`$ and $`X_{50}=cB_z^2^{1/2}\mathrm{\Omega }/2\pi en_ev_{\mathrm{A}z}^2`$. Then the size of $`\eta `$ and $`n_e`$ after 50 orbits are obtained from $`Re_{M,50}`$ and $`X_{50}`$, respectively. For simplicity, the electron abundance $`n_e/n_n`$ is assumed to be constant in time and uniform in space throughout the calculation. Since the pressure has little effect on the evolution of the MRI, we find it does not matter whether or not the pressure is changed to be consistent with the reduced temperature of the disk after 50 orbits. We have not changed the pressure in the models reported here.

## 4 NUMERICAL RESULTS

### 4.1 Models without a Net Flux of the Magnetic Field

First, models started with a zero net flux $`B_z`$ are studied in this subsection. All the calculated models are listed in Table 1. To obtain MHD turbulence in the hot state, model S is calculated for 50 orbits using the ideal MHD equations. The models of quiescent disk are labeled SA, SB, SC, and SD. For each model, ohmic dissipation and the Hall term are included after 50 orbits. The size of the parameters $`Re_{M,50}`$ and $`X_{50}`$ are listed in the Table. The magnetic diffusivity and the electron abundance at quiescence are calculated using the magnetic field strength of model S at 50 orbits. The turbulent stress is the sum of the Maxwell stress, $`w_M=B_xB_y/4\pi `$, and the Reynolds stress, $`w_R=\rho v_x\delta v_y`$, where $`\delta v_y`$ is the perturbed azimuthal velocity from the Keplerian flow. The $`\alpha `$ parameter of Shakura & Sunyaev (1973) is given by $`\alpha =(w_M+w_R)/P_0`$. The time- and volume-averaged stresses ($`w_M`$ and $`w_R`$) are listed in Table 1. The time averages are taken from 75 to 100 orbits.

Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks shows the time evolution of the volume-averaged Maxwell stress $`w_M/P_0`$, where time is measured in orbits $`t_{\mathrm{rot}}=2\pi /\mathrm{\Omega }`$. The ideal MHD model (S) is shown by the dashed curve. The saturation level of the Maxwell stress in model S is of the order of 0.01 throughout the evolution. The solid curves after 50 orbits are for the models that include the nonideal MHD effects. From top to bottom, they are models SA ($`T=2800`$ K), SB (2200 K), SC (1800 K), and SD (1600 K).

The suppression of the MRI is mainly caused by ohmic dissipation (Sano & Stone 2002b). When ohmic dissipation is inefficient ($`Re_{M,50}>100`$), the amplitude of the MHD turbulence maintains a level similar to the ideal MHD case. At $`Re_{M,50}=10`$, the Maxwell stress is reduced by a factor of 4, but still has a significant amplitude ($`w_M/P_0=4.7\times 10^3`$). However, if $`Re_{M,50}1`$ or $`T1800`$ K, the MHD turbulence dies away very quickly. As the temperature decreases, the amplitude of the magnetic stress decreases dramatically. For the model SD, the magnetic field and stress are almost vanished at the end of the calculation, i.e., 50 orbits after the transition. The duration of quiescence in dwarf novae is typically more than a few tens of days, which corresponds to over 1000 orbits. Thus, the magnetic field in the disks can disappear locally during quiescence if the temperature is less than about 1500 K.

Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks shows the saturation level of the Maxwell $`w_M`$ and Reynolds stress $`w_R`$ as a function of the magnetic Reynolds number $`Re_{M,50}`$. The stresses are normalized in terms of the initial pressure $`P_0`$, so that they have the same dimension as the $`\alpha `$ parameter. The upper and lower arrows indicate the Maxwell and Reynolds stress for the ideal MHD model ($`Re_M\mathrm{}`$). As shown in Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks, active MHD turbulence can be sustain by the MRI if the magnetic Reynolds number is much larger than unity. Angular momentum transport is dominated by the Maxwell stress, which is about 4 times larger than the Reynolds stress.

When $`Re_{M,50}1`$, the MRI is suppressed by the ohmic dissipation and MHD turbulence dies away immediately. The saturation level of the Maxwell stress is quite small ($`w_M/P_010^6`$) compared with the less resistive runs (S, SA, and SB). This dependence on the magnetic Reynolds number is consistent with the results of previous work (Fleming et al 2000; Sano & Stone 2002b). While the Maxwell stress decreases as the temperature decreases, it is found that the Reynolds stress has an amplitude which is comparable to or larger than the magnetic stress. Since we do not include an explicit kinematic viscosity in our calculations, and since the numerical viscosity of our algorithm is much lower than the explicit magnetic diffusivity added to model quiescence, epicyclic oscillations remain even after the magnetic field disappears. Although the amplitude of the Reynolds stress is very small ($`\alpha w_R/P_010^6`$), these motions, which are remnants of the turbulence during outburst, could give non zero accretion stress during quiescence. However, such motions may be swamped by global effects which cannot be studied in our local simulations.

The reduction of the Maxwell stress due to the nonideal MHD effects is independent of the numerical resolution. We have performed the same simulations as shown by Figures A Local One-Zone Model of MHD Turbulence in Dwarf Nova DisksA Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks with double the resolution ($`64\times 256\times 64`$). The saturation level of the stress in the high-resolution runs are listed in Table 1. For the ideal MHD case, the Maxwell stress in the high-resolution run is in fact slightly smaller than that in the standard run. For less resistive models (SA and SB), the Maxwell stress is always larger than the Reynolds stress and the saturation level of the stress is almost unrelated with the magnetic Reynolds number. When $`Re_M1`$, on the other hand, the ohmic dissipation suppresses the MRI and MHD turbulence. The Maxwell stress is $`w_M/P_0=1.0\times 10^6`$ and $`1.0\times 10^{10}`$ for model SC and SD, respectively. For these cases, angular momentum transport is dominated by the Reynolds stress with a similar amplitude to the standard-resolution runs ($`\alpha 10^6`$).

### 4.2 Models with a Net Flux of the Vertical Field

Next we consider the effect of the initial field geometry. Models started with a uniform $`B_z`$ are listed in Table 2, along with the model parameters and saturation level of the Maxwell and Reynolds stress. With this field geometry, the saturation level of $`\alpha `$ in the ideal MHD run (Z) is of the order of 0.1. This is an order of magnitude larger than that in the zero net flux $`B_z`$ model (S). For models with a net flux of vertical field, a two-channel flow forms repeatedly at the nonlinear phase of the MRI (Sano & Inutsuka 2001). The channel flow is destroyed through magnetic reconnection triggered by the parasitic instability (Goodman & Xu 1994). Since the magnetic energy and stress increase exponentially when the channel flow is growing, the time evolution of the magnetic stress shows frequent spike-shaped excursions. The growth of the two-channel flow can be seen only in the uniform $`B_z`$ runs, and this causes the relatively larger magnetic stress compared with the zero net flux $`B_z`$ runs (Sano & Stone 2002b).

For models of a quiescent disk (ZA, ZB, ZC, and ZD), we choose the same sets of $`Re_{M,50}`$ and $`X_{50}`$ as for the zero net flux $`B_z`$ runs. In Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks, the Maxwell and Reynolds stress are depicted by the filled and open circles, respectively. The dependence of the stresses on the size of the nonideal effects is qualitatively quite similar to that for the zero net flux $`B_z`$ runs. When the magnetic Reynolds number is larger than 10, the stress is almost independent of $`Re_{M,50}`$ and $`X_{50}`$. The difference in the Maxwell stress between the ideal MHD (Z) and the $`Re_{M,50}=10`$ (ZB) runs is only 10 %. The Maxwell stress always dominates the Reynolds stress by a factor of 5.

If $`Re_{M,50}1`$, on the other hand, the Maxwell stress becomes inefficient ($`w_M/P_010^6`$) as a result of the suppression of the MRI. For models ZC and ZD, the spatial fluctuations in the magnetic field disappear soon after 50 orbits and the field geometry returns to a uniform $`B_z`$. The $`\alpha `$ parameter is dominated by the Reynolds stress, and the saturation level of $`w_R/P_0`$ is larger than those in the zero net flux $`B_z`$ models (SC and SD). We also calculated a model with $`Re_{M,50}=0.01`$ and $`X_{50}=100`$ ($`T=1400`$ K). In this case, the critical wavelength for the MRI is longer than the scale height of the disk (i.e., the simulation box size), because ohmic dissipation stabilizes the shorter-wavelength fluctuations. The saturation level of the Maxwell and Reynolds stress are $`w_M/P_0=3.6\times 10^6`$ and $`w_R/P_0=2.6\times 10^3`$, respectively. The Reynolds stress still has the largest amplitude, but this is a remnant of the turbulence during outburst. Once again, we emphasize the amplitude of the Reynolds stress may be different in global calculations.

### 4.3 Decay Timescale of MHD Turbulence

When the magnetic Reynolds number is less than unity, MHD turbulence immediately dies away after the transition. We define the decay time $`t_{\mathrm{dec}}`$ of the turbulence as the time it takes for the magnetic energy to fall to 1/e of its value at 50 orbits. The decay time for models with $`Re_{M,50}=1`$ (SC and ZC) is $`t_{\mathrm{dec}}/t_{\mathrm{rot}}=0.36`$ and 0.12, where $`t_{\mathrm{rot}}`$ is the orbital time. The $`Re_{M,50}=0.1`$ runs show much faster decay; $`t_{\mathrm{dec}}/t_{\mathrm{rot}}=0.028`$ and 0.0077 for models SD and ZD, respectively. The decay time in the high-resolution runs is slightly shorter than that in the standard runs ($`t_{\mathrm{dec}}/t_{\mathrm{rot}}=0.16`$ for model SC and $`t_{\mathrm{dec}}/t_{\mathrm{rot}}=0.012`$ for model SD). Thus all the decay times we obtained are less than an orbit.

The dissipation timescale of the magnetic field is given by $`t_{\mathrm{dis}}=L^2/\eta `$, where $`L`$ is a typical length scale. The turbulent motion is driven most effectively at the most unstable wavelength for the MRI ($`\lambda v_\mathrm{A}/\mathrm{\Omega }`$), which is much smaller than the disk scale height. We find that the decay times $`t_{\mathrm{dec}}`$ obtained by the numerical simulations are equivalent to the dissipation timescale with $`L=v_\mathrm{A}/\mathrm{\Omega }`$;

$$\frac{t_{\mathrm{dis}}}{t_{\mathrm{rot}}}=\frac{v_\mathrm{A}^2/\mathrm{\Omega }^2}{\eta }\frac{\mathrm{\Omega }}{2\pi }=0.16Re_{M,50}.$$
(9)
This relation is valid only when the magnetic Reynolds number is less than unity.

## 5 DISCUSSION

### 5.1 The Critical Temperature and Transition Timescale

When the magnetic Reynolds number falls to unity, MHD turbulence is suppressed by ohmic dissipation independent of the amplitude of the Hall term (Sano & Stone 2002a,b). The critical temperature corresponding to $`Re_M=1`$ is about 2000 K for our model. As seen from equation (5), however, the critical temperature depends on the neutral density $`n_n`$ and the distance from the central star $`r`$. Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks shows the critical temperature $`T_{\mathrm{crit}}`$ as a function of $`r`$ for the cases of the neutral density $`n_n=10^{17}`$, $`10^{18}`$, and $`10^{19}`$ cm-3. We assume the Alfvén speed $`v_A=10^5`$ cm s-1 and the mass of the central star $`M=M_{}`$. At a given temperature, the electron abundance decreases as the neutral density increases, so that the critical temperature is slightly higher when the density is higher. The MRI wavelength ($`v_\mathrm{A}/\mathrm{\Omega }`$) is shorter at a smaller radius, because the angular velocity is larger. Therefore, the MRI can be suppressed with a smaller amount of diffusivity $`\eta `$ in the inner regions of the disk. The critical temperature is about 2500 and 1500 K at $`r=10^9`$ and $`10^{11}`$ cm, respectively. Since the magnetic Reynolds number, or the electron abundance, has a steep dependence on the temperature, the range of the critical temperature is quite narrow and the difference is at most a factor of 2 over the entire disk.

In our numerical analysis, we assume a fixed temperature during quiescence. Thus the electron abundance at quiescence is assumed to be constant. In a real disk, the critical temperature could be a little higher due to a runaway decay of MHD turbulence (Menou 2000). The decrease of the temperature makes both the electron abundance lower and nonideal MHD effects more important. The turbulent stress is therefore suppressed further, so that the temperature falls even more. However, a very weak dependence of the stress on the temperature can be seen in our models when $`T3000`$ K. This indicates that the runaway temperature should be less then 3000 K. Although a self-consistent study of the evolution of the temperature and electron abundance is needed to obtain a real $`T_{\mathrm{crit}}`$, Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks may give a fairly good estimate of the critical temperature in dwarf nova disks.

We also assume the transition to quiescence occurs instantaneously. In actual systems, however, this may proceed at the thermal timescale $`t_{\mathrm{th}}`$, where $`t_{\mathrm{th}}/t_{\mathrm{rot}}1/\alpha 100`$ (Frank, King, & Raine 1992). Suppose the temperature decreases gradually. As long as the temperature is higher than the critical value, we have shown that the accretion stress is nearly constant. But, once the temperature drops below this critical value, the MRI is suppressed, the turbulence decays within an orbit, and the $`\alpha `$ parameter drops by more than 2 orders of magnitude. Therefore, even though the decrease of the temperature is gradual, the transition timescale of the stress may be very much shorter.

### 5.2 Onset of the Next Outburst

Dwarf nova systems undergo recurrent outbursts. The matter from the secondary star accumulates in the disk during quiescence. When the surface density exceeds a critical value at some region in the disk, the thermal instability sets in and the temperature goes up. Here we examine the behavior of the MRI at the transition to the hot state.

Figure A Local One-Zone Model of MHD Turbulence in Dwarf Nova Disks shows the time evolution of the Maxwell and Reynolds stress for model ZC until 100 orbits. After 100 orbits, the induction equation for the ideal MHD is used in order to imitate the onset of the next outburst. As seen from the figure, the growth of the MRI starts immediately at the transition to the outburst and the magnetic field is amplified exponentially. In a few orbits, the Maxwell stress reaches the same saturation level as that before 50 orbits. The key to this evolution is the existence of a net vertical flux. During quiescence (from 50 to 100 orbits), MHD turbulence is suppressed, but the magnetic flux within the shearing box is conserved throughout the evolution. Thus the field geometry at 100 orbits is a nearly uniform vertical field. When the nonideal MHD effects are turned off, the linear growth of the MRI starts and MHD turbulence is initiated.

This evolution may be different if the system has no net magnetic flux. For example, if the magnetic field is completely dissipated during quiescence (e.g., model SD), MHD-driven activity cannot occur during the next outburst. However, if enough magnetic field remains for there to be unstable modes even in a small region of the disk, MHD turbulence develops in that region and finally may spread over a large part of the disk (Hawley & Balbus 1991). Therefore, it is quite important to understand the global structure of the magnetic field during quiescence, including the magnetosphere of the white dwarf and the magnetic field supplied by accreting material from the secondary star (Meyer & Meyer-Hofmeister 1999).

### 5.3 Beyond the Local Model

Since we focus on the local behavior of MHD turbulence in this paper, neither global structural effects nor global instabilities have been considered. In real disks, global effects such as spiral shocks (Sawada et al. 1986) may also transport angular momentum.

From our numerical results, the dependence of the $`\alpha `$ parameter associated with local turbulence on temperature in a dwarf nova disk can be summarized as follows. The contribution of the MRI to the $`\alpha `$ parameter is significant during outburst. When the temperature is higher than the critical value $`T_{\mathrm{crit}}2000`$ K, the accretion stress is dominated by the Maxwell stress, with $`\alpha 0.1`$ and 0.01 for the cases with and without net vertical flux, respectively. At quiescence, on the other hand, $`\alpha `$ is very sensitive to the temperature of the disk. If $`TT_{\mathrm{crit}}`$, the stress is more than 2 orders of magnitude smaller than that during outburst.

There are several reasons why global effects may influence these results. Firstly, the presence of a net vertical flux can affect the amplitude of the turbulence stress during outburst, that is $`\alpha _{\mathrm{hot}}`$. In agree with the previous works using the local shearing box (e.g., Hawley et al. 1995; 1996), our simulations have shown that an sufficient amount of angular momentum transport in the hot state requires the existence of a net vertical field. On the other hand, recent global simulations obtained a higher stress ($`\alpha 0.1`$) even without a net magnetic flux (Stone & Pringle 2001; Hawley, Balbus, & Stone 2001). There is no observational constraint yet in the structure and strength of the magnetic field in dwarf nova disks. Thus, theoretical modeling of the global structure of the magnetic field is quite important for a more detailed analysis of $`\alpha `$ at the hot state. Secondly, the ionization fraction of the disk gas, and therefore the importance of non-ideal MHD effects, is extremely sensitive to temperature. Thus global models in which the radial and vertical temperature profiles in the disk are calculated self-consistently with the local heating rate due to turbulent stresses and cooling rate due to radiation are required. Furthermore, at the surface of the disk, the radiation from a hot white dwarf could be important source of heating and ionization (Hameury, Lasota, & Dubus 1999; Menou 2002). Finally, even when Ohmic dissipation completely suppresses the MRI, small amplitude motions remain in the disk and can give a non-zero Reynolds stress. In the local model, these motions are damped on a viscous timescale, which may be very long. However, they may be affected by global density waves in the disk or tidally induced spiral shocks. Thus, global MHD simulations are an important next step for understanding the dynamics of dwarf nova disks.

We would like to thank John Cannizzo, Kristen Menou, and Neal Turner for helpful discussions. Computations were carried out on the VPP5000 at the National Astronomical Observatory of Japan, and the VPP700 at the Subaru Telescope, NAOJ. This work was supported by grants from the NASA OSS program and the NSF.


format: text

text:
# Solar neutrino oscillation parameters after first KamLAND results

## I Introduction

The year 2002 is likely to be remembered as the annus mirabilis of solar neutrino physics. On April 20, direct and highly significant evidence for $`\nu _e`$ flavor change into active states was announced by the Sudbury Neutrino Observatory (SNO) experiment AhNC , crowning a four-decade long Deca  series of beautiful observations Cl98 ; Ksol ; Ab02 ; Ha99 ; Ki02 ; Fu01 ; Fu02 ; AhCC ; AhDN  of the solar $`\nu _e`$ flux deficit Ba89 ; BP00 . On October 8, the role of solar neutrino physics in shaping modern science was recognized through the Nobel Prize jointly awarded to Raymond Davis, Jr., and Masatoshi Koshiba, for their pioneering contributions to the detection of cosmic neutrinos Nobe . Finally, on December 6, clear “terrestrial” evidence for the oscillation solution to the solar neutrino deficit was reported by the Kamioka Liquid scintillator AntiNeutrino Detector (KamLAND), through the observation of long-baseline reactor $`\overline{\nu }_e`$ disappearance Kam1 . The seminal idea of studying lepton physics by detecting solar Po46 ; Al49  and reactor Po46 ; Al49 ; Re53  neutrinos keeps thus bearing fruits after more than 50 years.

The KamLAND observation of $`\overline{\nu }_e`$ disappearance Kam1  confirms the current interpretation of solar neutrino data AhDN ; Fu02 ; AllS ; GetM ; Las3  in terms of $`\nu _e\nu _{\mu ,\tau }`$ oscillations induced by neutrino mass and mixing Pont ; Maki , and restricts the corresponding parameter space $`(\delta m^2,\theta _{12})`$ within the so-called large mixing angle (LMA) region. In this region, globally favored by solar neutrino data Smir , matter effects Matt ; Barg  in adiabatic regime Adia ; Matt  are expected to dominate the dynamics of flavor transitions in the Sun (see, e.g., LasA ). The KamLAND spectral data appear to exclude some significant portions of the LMA solution Kam1 , where the predicted spectrum distortions Mu02 ; Marf ; Barb ; Go01 ; Vi02 ; GoKL ; Al02 ; Ho02 ; Go02 ; Las3  would be in conflict with observations Kam1 .

In this paper we analyze the first KamLAND spectral data Kam1  and combine them with current solar neutrino data Cl98 ; Ab02 ; Ha99 ; Ki02 ; Fu02 ; AhNC ; AhCC ; AhDN , assuming two- and three-flavor oscillations of active neutrinos Las3 , in order to determine the surviving sub-regions of the LMA solution. In the analysis we include the CHOOZ reactor data CHOO  and, in the $`3\nu `$ case, also the relevant constraints on the larger mass gap $`\mathrm{\Delta }m^2\delta m^2`$, coming from the Super-Kamiokande (SK) atmospheric Shio  and KEK-to-Kamioka (K2K) accelerator K2Ke ; K2Ks  neutrino experiments, according to the approach developed in Las3 .111Our notation for the squared neutrino mass spectrum is $`(m_1^2,m_2^2,m_3^2)=(\delta m^2/2,+\delta m^2,\pm \mathrm{\Delta }m^2)`$ as in Las3 , the sign of $`\mathrm{\Delta }m^2`$ being associated to normal $`(+)`$ or inverted $`()`$ hierarchy. Unless otherwise noticed, normal hierarchy is assumed. The mixing angles $`\theta _{ij}`$ are defined as in the quark case Las3 ; PDGr .

In the $`2\nu `$ case, we find that the inclusion of the KamLAND spectrum basically splits the LMA solution into two sub-regions at “small” and “large” $`\delta m^2`$, which we call LMA-I and LMA-II, respectively (the LMA-I solution being preferred by the data). Such regions are only slightly modified in the presence of $`3\nu `$ mixing, namely, for nonzero values of the mixing angle $`\theta _{13}`$. We also present updated bounds in the $`3\nu `$ parameter space $`(\delta m^2,\theta _{12},\theta _{13})`$.

The structure of this paper is as follows. In Sec. II we describe our analysis of the KamLAND data. In Sec. III we combine such measurements with solar and CHOOZ data, assuming $`2\nu `$ mixing. In Sec. IV we extend the analysis to $`3\nu `$ mixing, including the SK+K2K terrestrial neutrino constraints. We draw our conclusions in Sec. V.

## II Kamland data input

In our KamLAND data analysis, we use the absolute spectrum of events reported in Kam1 , taken above a background-safe analysis threshold of 2.6 MeV in visible energy $`E`$.222The “visible” or “prompt” energy $`E`$ is defined as $`E=E_\nu (M_nM_p)+m_eE_\nu 0.8`$ MeV. The events below such threshold might contain a significant component of geological $`\overline{\nu }_e`$’s GeoN , whose large normalization uncertainties are poorly constrained at present by the KamLAND data themselves Kam1 . Above 2.6 MeV, a total of 54 events is found (including at most one possible background candidate), against 86.8 events expected from reactors Kam1 .

The observed energy spectrum of events is analyzed as in Las3 , with the following improvements. According to Kam1 , we adopt: (1) Relative fuel components 235U : 238U : 239Pu : 241Pu = 0.568 : 0.078 : 0.297 : 0.057; (2) Absolute normalization of 86.8 events for no oscillations; (3) Energy resolution width equal to $`7.5\%/E^{1/2}`$; (4) Thirteen bins in visible energy above 2.6 MeV, with 0.425 MeV width. The experimental and theoretical (oscillated) number of events in each bin are denoted as $`N_i^{\mathrm{expt}}`$ and $`N_i^{\mathrm{osc}}`$, respectively $`(i=1,\mathrm{},13)`$.

Concerning the systematic errors, the information in Kam1  does not allow to trace the magnitude and the bin-to-bin correlations of each component. We have then approximately grouped such uncertainties into two main components, labelled by the index $`k=1,\mathrm{\hspace{0.17em}2}`$: (1) energy scale error; and (2) overall normalization error. The first (second) uncertainty basically shifts the spectrum along the energy (event) coordinate. More precisely, the effect of the energy scale error $`(k=1)`$ is estimated by taking the (symmetrized) fractional differences $`c_i^1`$ between the values of $`N_i^{\mathrm{osc}}`$ calculated with and without shifts of $`\pm 1.91\%`$ Kam1  in the true visible energy, for each point in the oscillation parameter space. The normalization uncertainty $`(k=2)`$ is here assumed to collect all the other error components listed in Kam1 , up to a total of about 6%: $`c_i^2=0.06`$. The systematic shifts $`c_i^k`$ of the $`N_i^{\mathrm{osc}}`$’s are then implemented through the “pull approach” described in GetM , namely, through linear deviations $`N_i^{\mathrm{osc}}N_i^{\mathrm{theo}}=N_i^{\mathrm{osc}}(1+_kc_i^k\xi _k)`$, where the $`\xi _k`$’s are univariate Gaussian random variables. In this way, systematic correlations are also accounted for GetM .

Concerning the statistical errors, the presence of bins with few or zero events requires Poisson statistics, that we approximately implement through the $`\chi ^2`$-like recipe suggested in the Review of Particle Properties PDGr . A quadratic penalty in the pulls $`\xi _k`$ is introduced to account for the systematics GetM ; Lind . The final $`\chi ^2`$ function for KamLAND (KL) is thus

$$\chi _{\mathrm{KL}}^2=\underset{\{\xi _k\}}{\mathrm{min}}\left[\underset{i=1}{\overset{13}{}}2\left(N_i^{\mathrm{theo}}N_i^{\mathrm{exp}}+N_i^{\mathrm{exp}}\mathrm{ln}\frac{N_i^{\mathrm{exp}}}{N_i^{\mathrm{theo}}}\right)+\underset{k=1,2}{}\xi _k^2\right],$$
(1)
where the $`i`$-th logarithmic term is dropped if $`N_i^{\mathrm{exp}}=0`$ PDGr . By expanding the first sum in the $`\chi ^2`$ at first order in the shifts $`c_i^k`$, the minimization in Eq. (1) becomes elementary.

Some final remarks are in order. A detailed comparison of the above function $`\chi _{\mathrm{KL}}^2`$ with the one adopted by the KamLAND collaboration is not currently possible, since the latter is given in a symbolic form in Kam1 . However, one can observe that: (1) The KamLAND collaboration splits the total rate and spectrum shape information, while we prefer to use the absolute spectrum; (2) The KamLAND collaboration can take into account a (very small) background component above 2.6 MeV and four sources of systematics while, for a lack of information, we neglect such small background, and group the systematics into two main sources. For such reasons, we do not expect perfect quantitative agreement between the KamLAND official oscillation analysis and ours. In fact, we find that the 95% C.L. contours of the rate+shape analysis in Fig. 6 of Kam1  are accurately reproduced by our Eq. (1) at a slightly lower C.L. ($`90\%`$). This loss of statistical power appears tolerable to us, in the light of the limited experimental information which is currently available.333A better evaluation and decomposition of systematic effects and an improved $`\chi _{\mathrm{KL}}^2`$ definition will be possible when more detailed KamLAND information on spectrum shape errors will become public.

## III $`2\nu `$ analysis

The updated $`2\nu `$ analysis of current solar+CHOOZ neutrino data, as performed in Las3 , is presented here for the sake of completeness. The fit includes 81 solar neutrino observables GetM ; Las3  and 14 CHOOZ spectrum bins CHOO ; Las3 , for a total of 95 data points. The best-fit point and its $`\chi _{\mathrm{min}}^2`$ value are given in the first row of Table I. The $`\mathrm{\Delta }\chi ^2`$ expansion around the minimum, relevant for the estimation of the oscillation parameters $`(\delta m^2,\mathrm{sin}^2\theta _{12})`$, is shown in Fig. 1, where we have restricted the $`\delta m^2`$ range to the only three decades relevant for the LMA solution and for the following KamLAND analysis. Notice that the scale is linear in the $`\mathrm{sin}^2\theta _{12}`$ variable.

The $`2\nu `$ analysis of KamLAND is performed by using Eq. (1), which gives the absolute $`\chi _{\mathrm{min}}^2`$ value reported in the second row of Table I. For later purposes, we also quote the second best fit parameters. Expansion around the absolute minimum gives the C.L. contours shown in Fig. 2.444No oscillations (basically, oscillations with $`\delta m^210^6`$ eV2) provide a very bad fit, $`\mathrm{\Delta }\chi ^2=14.1`$. There appears to be a “tower” of solutions which tend to merge and become indistinguishable for increasing $`\delta m^2`$; the lower three ones are, however, rather well separated at 90% C.L. Notice that our allowed regions are slightly larger (i.e., less constraining) than those in the rate+shape analysis of Kam1 , as explained in Sec. II. One of the two octant-symmetric best fits points in Fig. 2 (black dots) is remarkably close to the best fit in Fig. 1 (see also Table I). The difference in location with respect to the KamLAND official best-fit point at $`(\delta m^2/\mathrm{eV}^2,\mathrm{sin}^2\theta _{12})=(6.9\times 10^5,\mathrm{\hspace{0.17em}0.5})`$ Kam1  is not statistically significant, amounting to a variation $`\mathrm{\Delta }\chi ^21`$.

The combination of the solar+CHOOZ results in Fig. 1 with the KamLAND results in Fig. 2 gives the global $`2\nu `$ results shown in Fig. 3, which represent the main result of this work. Two rather distinct solutions, that we label LMA-I and LMA-II, are seen to emerge. They are basically located at the intersection of the LMA solution in Fig. 1 with two of the well-separated KamLAND solutions in Fig. 2, and are characterized by the mass-mixing parameters and $`\chi ^2`$ values given in the last two rows of Table 1. The LMA-I solution is clearly preferred by the data, being close to the best fit points of both solar+CHOOZ and KamLAND data. The LMA-II solution is located at a $`\delta m^2`$ value about twice as large as for the LMA-I, but is separated from the latter by a modest $`\mathrm{\Delta }\chi ^2=5.4`$ difference (dominated by solar neutrino data). Indeed, if we conservatively demand a 99.73% C.L. for the allowed regions, the LMA-I and LMA-II solutions appear to merge (and extend towards $`\delta m^23\times 10^4`$ eV2) in a single broad solution. In any case, at any chosen C.L., the allowed regions of Fig. 3 are significantly smaller than those in Fig. 1. Therefore, with just 54 initial events, the KamLAND experiment is not only able to select the LMA region as the solution to the solar neutrino problem, but can also significantly restrict the corresponding oscillation parameter space. With several hundred events expected in the forthcoming years, there are thus very good prospects to refine the $`(\delta m^2,\mathrm{sin}^2\theta _{12})`$ parameter estimate Las3

The most important task for the next future appears to be the confirmation of one of the two solutions in Fig. 3 through higher statistics and, possibly, lower analysis threshold. Figure 4 shows the absolute KamLAND energy spectra of reactor events predicted at the LMA-I and LMA-II global best-fit points (last two rows of Table I), together with the no-oscillation spectrum, with the current event normalization. The KamLAND data (with statistical errors only) are also superposed above the analysis threshold (2.6 MeV). It can be seen that the main difference between the two oscillated spectra is the position of the spectrum peak, which is roughly aligned with the no-oscillation position for the LMA-II case, while it is shifted at higher energies for the LMA-I case. This feature might be a useful discrimination tool in the next future. From Fig. 4, it appears also that the largest difference between the LMA-I and LMA-II spectra in KamLAND occurs just in the first bin below the current analysis threshold. Therefore, a better understanding of the background below 2.6 MeV will be highly beneficial. In any case, an increase in statistics by a factor of a few appears necessary to disentangle the two oscillated spectra in Fig. 4.

## IV $`3\nu `$ analysis

The updated $`3\nu `$ analysis of solar+CHOOZ neutrino data, including the constraints on the large (“atmospheric”) squared mass gap $`\mathrm{\Delta }m^2`$ coming from SK atmospheric Shio  and K2K accelerator neutrino data K2Ke ; K2Ks , as performed in Las3 , is reported here for the sake of completeness.555In Las3 , the combination of SK and K2K constraints on the $`\mathrm{\Delta }m^2`$ parameter was obtained by adding the preliminary $`\chi ^2`$ functions from the two collaborations, as graphically presented in Shio ; K2Ke . We anticipate that such combination appears to be in good agreement with a more refined and joint reanalysis of SK atmospheric and K2K data Prog . Figure 5 shows the pre-KamLAND $`3\nu `$ results Las3  in the parameter space $`(\delta m^2,\mathrm{sin}^2\theta _{12},\mathrm{sin}^2\theta _{13})`$ relevant for the LMA solution, shown through sections at four different values of $`s_{13}^2=\mathrm{sin}^2\theta _{13}`$.666We remind that, in Fig. 5 and in the following $`3\nu `$ figures, the information on $`\mathrm{\Delta }m^2`$ is projected away by taking $`\mathrm{min}_{\mathrm{\Delta }m^2}(\chi ^2)`$, see Las3 . For later purposes, we only notice that the upper bound on $`\delta m^2`$ becomes stronger for higher $`\mathrm{sin}^2\theta _{13}`$, as a result of the CHOOZ constraints. For a discussion of such $`(\delta m^2,\mathrm{sin}^2\theta _{13})`$ anti-correlation, and for the pre-KamLAND solutions below the LMA, see Las3  and references therein.

Figure 6 shows the results of our $`3\nu `$ fit to the KamLAND data, in the same format as in Fig. 5. The allowed $`3\nu `$ regions appear to get slightly enlarged (especially in the mixing parameter $`\theta _{12}`$) for increasing values of $`\mathrm{sin}^2\theta _{13}`$, as expected GoKL . In fact, for $`\theta _{13}0`$, part of the KamLAND event disappearance is explained by averaged oscillations driven by $`\mathrm{\Delta }m^2`$ in the $`(\nu _e,\nu _3)`$ sector, so that the overall oscillation amplitude $`\mathrm{sin}^22\theta _{12}`$ in the $`(\nu _1,\nu _2)`$ sector is allowed to reach smaller values, and the $`\mathrm{sin}^2\theta _{12}`$ range is correspondingly enlarged GoKL . The absolute $`\chi ^2`$ minimum in Fig. 6 is reached for $`\theta _{13}0`$, but such preference is not statistically significant, as expected Marf ; indeed, we find a mere $`\mathrm{\Delta }\chi ^2=0.4`$ increase of the KamLAND $`\chi _{\mathrm{min}}^2`$ value for fixed $`\mathrm{sin}^2\theta _{13}=0.06`$. Concerning the variations of the $`(\delta m^2,\mathrm{sin}^2\theta _{12})`$ best-fit coordinates for fixed values of $`\mathrm{sin}^2\theta _{13}`$, we find that the $`\delta m^2`$ value is stable (and equal to $`7.3\times 10^5`$ eV2, as in the $`2\nu `$ case), while $`\mathrm{sin}^2\theta _{12}`$ decreases from 0.335 to 0.225 when increasing $`\mathrm{sin}^2\theta _{13}`$ from 0 to 0.06, in good qualitative agreement with the expectations discussed in GoKL .

However, the additional spread in $`\mathrm{sin}^2\theta _{12}`$ induced by nonzero $`\theta _{13}`$ currently does not play any relevant role when KamLAND is combined with world neutrino data, for at least two reasons: (1) Pre-KamLAND bounds from solar+terrestrial data currently dominate the constraints on $`\mathrm{sin}^2\theta _{12}`$, as evident from a comparison of Fig. 5 and Fig. 6; (2) The likelihood of genuine $`3\nu `$-induced effects rapidly decreases with increasing $`\theta _{13}`$, because of the strong upper bounds on such mixing parameter Las3 . Therefore, we do not expect any significant enlargement of the $`\mathrm{sin}^2\theta _{12}`$ allowed range from the pre- to the post-KamLAND $`3\nu `$ analysis, despite the presence of such effects in KamLAND alone.

It should be noted that, strictly speaking, Fig. 6 is not exactly a “KamLAND only” analysis, since we have implicitly taken some pieces of information from terrestrial neutrino data. In particular, we have implicitly assumed in Fig. 6 that: (i) the “atmospheric” squared mass splitting $`\mathrm{\Delta }m^2`$ is sufficiently high to be unresolved in KamLAND, and (2) the relevant values of $`\mathrm{sin}^2\theta _{13}`$ are limited in the few percent range. It is interesting to study the effect of an explicit combination of such (atmospheric + CHOOZ) information with KamLAND data. The results are given in Fig. 7. This figure show that purely terrestrial neutrino data from atmospheric (SK), accelerator (K2K) and reactor (KamLAND + CHOOZ) neutrino experiments, by themselves, are now able to put both upper and lower bounds on the solar parameters $`(\delta m^2,\mathrm{sin}^2\theta _{12})`$. As previously noted in the comment to Fig. 5, the CHOOZ upper bound on $`\delta m^2`$ becomes stronger when $`\mathrm{sin}^2\theta _{13}`$ increases.777Analogously, the upper bound on $`\mathrm{sin}^2\theta _{13}`$ becomes stronger for increasing $`\delta m^2`$ BiCo . The slight octant asymmetry for $`\mathrm{sin}^2\theta _{13}>0`$ is due to the fact the corresponding $`3\nu `$ CHOOZ probability is not invariant under the change $`\theta _{12}\pi /2\theta _{12}`$ for fixed hierarchy QEIP  (assumed to be normal in Fig. 6). The octant differences would be swapped by inverting the hierarchy QEIP . In practice, however, such octant asymmetries are numerically irrelevant in the global fit which we now discuss.

Figure 8 shows the final $`3\nu `$ combination of world (solar+terrestrial) neutrino constraints, including pre-KamLAND data (Fig. 5) and the first KamLAND data (Fig. 6). For $`\mathrm{sin}^2\theta _{13}=0`$ (absolute best fit), we get the same LMA-I and LMA-II solutions reported in Fig. 3 and in Table I, modulo the expected widening induced by one additional degree of freedom in the $`\mathrm{\Delta }\chi ^2`$’s associated to each C.L. contour. This widening gives, at 99% C.L., marginal allowance for a third solution with $`\delta m^22.5`$$`3.2\times 10^4`$ eV2, which we call LMA-III.888 Notice that the LMA-I and LMA-III solutions differ by $`\mathrm{\Delta }\chi ^211`$. All solutions rapidly shrink for increasing $`\mathrm{sin}^2\theta _{13}`$, the LMA-I being the most stable and the last to disappear. The hierarchy is assumed to be normal in Fig. 8; the differences with respect to the inverted hierarchy case (not shown) are completely negligible.

Figure 9 represents a compendium of our post-KamLAND $`3\nu `$ analysis of all data (solar+terrestrial), in terms of the $`\mathrm{\Delta }\chi ^2`$ function for each of the three variables $`(\delta m^2,\mathrm{sin}^2\theta _{12},\mathrm{sin}^2\theta _{13})`$, the others being projected (minimized) away. This figure should be compared with the pre-KamLAND one in Las3 . The comparison shows that the bounds on $`\mathrm{sin}^2\theta _{13}`$ remain basically unaltered in the post-KamLAND era, while the mixing parameter $`\mathrm{sin}^2\theta _{12}`$ is clearly nailed down around the value $`0.3`$. The absolute and second best minima (LMA-I and LMA-II) are also evident in the left panel. A third minimum (LMA-III) at higher $`\delta m^2`$ is only marginally allowed.

If we focus on the best-fit (LMA-I) solution, and adapt parabolic functions around the absolute $`\chi ^2`$ minima in terms of $`\delta m^2`$ (in linear scale), $`\mathrm{sin}^2\theta _{13}`$, and $`\mathrm{sin}^2\theta _{13}`$, we obtain from the post-KamLAND global analysis the following approximate $`\pm 1\sigma `$ estimates ($`\mathrm{\Delta }\chi ^21`$) for the relevant solar $`3\nu `$ oscillation parameters,

$$\text{LMA-I }(1\sigma ):\{\begin{array}{c}\delta m^2(7.3\pm 0.8)\times 10^5\mathrm{eV}^2,\hfill \\ \mathrm{sin}^2\theta _{12}0.315\pm 0.035,\hfill \\ \mathrm{sin}^2\theta _{13}0.017.\hfill \end{array}$$
(2)
The above ranges are meant to show that we are not far from a 10% determination (at $`1\sigma `$) of both $`\delta m^2`$ and $`\mathrm{sin}^2\theta _{12}`$, but should not be quoted as a “summary” of the current $`3\nu `$ situation. The technically correct reference summary for our analysis is represented by the $`\chi ^2`$ functions of Fig. 9, which include the possibility of a second global best fit (LMA-II), and maybe of a third best fit (LMA-III) in $`\delta m^2`$. In particular, we stress that the LMA-II solution is currently quite acceptable from a statistical point of view, and that the history of the solar neutrino problem teaches us that one should not exclude a priori what may appear as the second best-fit at a particular time.

Similar cautionary remarks also apply to the use of the initial KamLAND data to refine the emerging indications of adiabatic matter effects in the Sun Adia ; Matt  through the approach developed in LasA . As stressed in LasA , such analysis will become sufficiently stable and constraining when a single solution (either the LMA-I or the LMA-II) will be (hopefully) clearly selected by KamLAND. In the meantime, further SNO data are also expected to increase our confidence in the occurrence of solar matter effects, see LasA .

We conclude our $`3\nu `$ analysis by observing that, from the point of view of future (very) long baseline (LBL) experiments, the “solar” neutrino parameters $`\delta m^2`$ and $`\theta _{12}`$ enter the subleading oscillation amplitudes of both the “golden” channel $`\nu _e\nu _\mu `$ Gold  and the “silver” channel $`\nu _e\nu _\tau `$ Silv  only through the combination $`\kappa =\delta m^2\mathrm{sin}2\theta _{12}`$, up to $`O(\kappa ^2)`$ terms included Gold ; Silv . Therefore, we think it useful to show in Fig. 10 the global (solar+terrestrial) bounds on this specific combination of parameters (in units of $`10^5`$ eV2 and in linear scale), all other $`3\nu `$ variables being projected away. The lower and higher minima in Fig. 10 correspond to the LMA-I and LMA-II cases, respectively. By adapting a parabola in the $`\pm 2\sigma `$ region around each of the two minima, we get the following $`1\sigma `$ ($`\mathrm{\Delta }\chi ^2=1`$) estimates for $`\kappa =\delta m^2\mathrm{sin}2\theta _{12}`$: $`\kappa 6.7\pm 0.7`$ (LMA-I) and $`\kappa 14.7\pm 1.7`$ (LMA-II), in units of $`10^5`$ eV2. In both cases, $`\kappa `$ appears to be already constrained with a $`11\%`$ accuracy. These approximate estimates—or, better, the accurate bounds in Fig. 10—may prove useful in prospective studies of $`\kappa `$-dependent subleading oscillation effects in LBL projects.

## V Summary and prospects

The KamLAND experiment has clearly selected the LMA region as the solution to the solar neutrino problem, and has further reduced the $`(\delta m^2,\mathrm{sin}^2\theta _{12})`$ parameter space for active neutrino oscillations. In this work, we have studied such parameter reduction in the the context of $`2\nu `$ and $`3\nu `$ oscillations, by using the limited KamLAND experimental information which is publicly available. In the $`2\nu `$ case, we find that the post-KamLAND LMA solution appears to be basically split into two sub-regions, LMA-I and LMA-II. The LMA-I solution, characterized by $`\delta m^27\times 10^5`$ eV2 and $`\mathrm{sin}^2\theta _{12}0.3`$, is preferred by the global fit. The LMA-II solution represents the second best fit, at about twice the value of $`\delta m^2`$ (see Table I). This situation is not significantly changed in the $`3\nu `$ case, for which we present a global post-KamLAND analysis of solar and terrestrial data in the $`(\delta m^2,\mathrm{sin}^2\theta _{12},\mathrm{sin}^2\theta _{13})`$ parameter space. There are good prospects to separate the LMA-I and LMA-II cases with future, higher-statistics KamLAND data, by looking at the peak of the energy spectrum, and by lowering the current analysis threshold (2.6 MeV) by at least $`0.5`$ MeV.

###### Acknowledgements.

This work was supported in part by INFN and in part by the Italian Ministero dell’Istruzione, Università e Ricerca through the “Astroparticle Physics” research project. We are grateful to G. Gratta for prompt information about the release of the first KamLAND results.


text:
# The effect of a nucleating agent on lamellar growth in melt-crystallizing polyethylene oxide

## Abstract

The effects of a (non co-crystallizing) nucleating agent on secondary nucleation rate and final lamellar thickness in isothermally melt-crystallizing polyethylene oxide are considered. SAXS reveals that lamellae formed in nucleated samples are thinner than in the pure samples crystallized at the same undercoolings. These results are in quantitative agreement with growth rate data obtained by calorimetry, and are interpreted as the effect of a local decrease of the basal surface tension, determined mainly by the nucleant molecules diffused out of the regions being about to crystallize. Quantitative agreement with a simple lattice model allows for some interpretation of the mechanism.

Nucleating agents have important applications in the industrial use of polymeric materials. As an example, mechanical properties such as hardness and elastic moduli can be controlled with these additives, which are able to affect crystallization kinetics rather dramatically 1 . In fact, not only the primary nucleation rate is influenced, but also the lamellar growth process via secondary nucleation. In the context of bulk crystallization, the latter mechanism has been envisaged for the first time within the framework of a phenomenological interpretation of the free growth kinetics’ data of indigo-nucleated isotactic polypropylene 2 ; 3 . Later, similar observations in nucleated high density polyethylene and polyethylene oxide (PEO) have been interpreted in the same way, with the further support of lattice calculations for the basal lamellar interface 4 .

The basic idea is that during growth, the nucleant molecules present in regions which are about to crystallize, diffuse towards the forming basal interfaces prior to crystallization (in all the cases above, the nucleating agent did not co-crystallize). These molecules cause a local decrease of the basal surface tension $`\sigma _e`$, so they tend to remain confined in the interfacial regions because diffusion towards the amorphous bulk would require an extra energy to restore a higher value of $`\sigma _e`$. In this situation the linear growth velocity $`v`$ in low to moderate undercooled melts, is still related to the secondary nucleation free enthalpy $`\mathrm{\Delta }G`$ by 5

$$v\mathrm{exp}\{\mathrm{\Delta }G/kT\}$$
(1)
(with $`T`$ the temperature), but the expression

$$\sigma _e=\sigma _{e0}\left(1+\frac{a}{\mathrm{\Delta }T}\right),$$
(2)
has to be substituted for $`\sigma _{e0}`$ in the formula

$$\mathrm{\Delta }G=4b_0\sigma \sigma _{e0}\frac{v_cT_m}{H_f\mathrm{\Delta }T}$$
(3)
customarily used for pure systems. In the above Eqs. 2 and 3, $`b_0`$ is the stem diameter, $`\sigma `$ is the lateral surface tension, $`\sigma _{e0}`$ is the basal surface tension in the absence of the nucleant, $`T_m`$ is the equilibrium melting point, $`H_f/v_c`$ is the enthalpy of fusion per unit volume, $`\mathrm{\Delta }T=T_mT`$ and $`a`$ ($`<0`$) is a coefficient proportional to the concentration of nucleating agent 2 ; 3 ; 4 .

Fig.1 shows the overall crystallization rates obtained by accurate DSC measurements on PEO nucleated with a number density $`n2.6\times 10^{19}`$ cm-3 of indigo, i.e. 1% wt fraction, in the interval $`331T334`$ K 4 . Fittings with Eqs. 1-3 are also reported \[i.e. assuming a constant nucleation density $`N(\mathrm{\Delta }T)`$ in this small range\] for the cases where either $`a`$ is let free to change, or is set equal zero. With $`a1.8`$ K ($`T_m=338`$ K) the fitting improves significantly 4 ; a close value is also obtained with a simple lattice model 4 . All this, and the fact that $`\sigma _{e0}\sigma _e`$ decreases with $`\mathrm{\Delta }T`$, implies that the mechanism described by Eq. 2 overwhelmes the effect of an increasing $`N(\mathrm{\Delta }T)`$ in this range (the irrelevance of other possibly competing mechanisms has been briefly discussed in 4 ).

Consider now the morphological features controlled by the nucleant in crystallized PEO. At fixed $`\mathrm{\Delta }T`$, number and final size of the spherulites are influenced by indigo through the primary nucleation rate. About a morphological counterpart of the secondary nucleation route to crystallization kinetics’ control, a fundamental role is played by Eq. 2. Since the lamellar fold length is related to $`\sigma _e`$ by

$$l=2\sigma _e\frac{v_cT_m}{H_f\mathrm{\Delta }T},$$
(4)
it is expected that the lamellar thickness of nucleated samples should be smaller than in pure samples crystallized at the same $`\mathrm{\Delta }T`$, and for $`\left|a/\mathrm{\Delta }T\right|`$ not too large:

$$\frac{l_0l}{l_0}\frac{a}{\mathrm{\Delta }T},$$
(5)
where $`l_0`$ is given by Eq. 4 with $`a=0`$.

SAXS measurements carried out on pure and nucleated PEO (with the same content of indigo as in Fig. 1) in fact confirm this prediction.

To assure homogeneity, nucleated samples were prepared from a solution of PEO 600 000 (Aldrich) and the appropriate amount of indigo in acetonitrile, which was kept under stirring for one day. The solvent was then very slowly evaporated at about 50 C. Pure PEO samples were prepared with similar procedure. Before crystallization, the samples were kept above the melting temperature for half an hour or more, and then cooled as rapidly as possible to the preset crystallization temperature.

The SAXS measurement apparatus consists of a Rigaku rotating anode generator, with a wavelength of 1.54 $`\mathrm{\AA }`$. The counter was 76 cm far from the sample (held the latter at room temperature), and the dimension of the beam at the counter was $`3\times 0.5`$ mm2, so that desmearing was necessary (point collimation measurements, with a $`0.5\times 0.5`$ mm2 beam, were performed as a check, and gave the same results, although with less statistics).

The lamellar thicknesses in the pure ($`l_0`$) and nucleated ($`l`$) samples have been estimated from the evaluation of the correlation function $`K(s)`$ 6 . The relatively small electron density difference between amorphous and crystalline regions in PEO rendered the data analysis a bit cumbersome. The spurious background forward scattering characterizing these systems (which is not unusual in general 7 ) was modelled by a local (i.e. $`q=0`$) gaussian distribution of density fluctuations, whose amplitude and width (the latter always found in the range $`4.57.5\times 10^3\mathrm{\AA }^1`$) were obtained by means of best fitting procedures. Fig. 2 shows the interface distributions $`K\mathrm{"}(s)`$ for some of the PEO samples. The first maximum in $`K\mathrm{"}(s)`$ is assigned to the interlamellar thickness, since i) its position is weakly dependent on $`\mathrm{\Delta }T`$, as can be seen from Tab. 1, and ii) the linear crystallinity thus found gets closer to the DSC one ($``$ 65 % at a rate -1 C/min in both nucleated and pure PEO). The resulting Gibbs-Thompson plot is in Fig. 3. The linear fit to the pure PEO data points extrapolates to a melting temperature $`T_c^{\mathrm{}}=339.3`$ K, i.e. rather close to the above value of $`T_m`$.

As the crystallization temperature decreased, complete thermalization was ever more difficult to achieve in PEO/indigo samples before the onset of crystallization. For nucleated samples crystallized at large $`\mathrm{\Delta }T`$ the preset temperature was only apparent. For example, at $`T=322`$ K the lamellar thickness of the pure and nucleated samples coincide, although a 0.1 relative difference had to be expected.

The quantity $`a/\mathrm{\Delta }T`$ as obtained from SAXS, fitting of Fig. 1 data (Cal.) and from lattice calculations, is reported in Tab. 1. The good agreement between the values of this parameter as obtained by SAXS (grown lamellae) and best fitting (growing lamellae), indicates that in the $`T`$-range considered the mechanism described by Eq. 2 plays a main role in the influence of the nucleant on the growth process (in particular, effects of nucleation density changes are negligible). Moreover, it supports the hypothesis that the nucleant remains confined at the basal interfaces during lamellar growth, without appreciable increase of the nucleant concentration in the remaining amorphous phase.

Some insight in the effect can be gained from the lattice description 4 . In the crystalline region the chain segment fractional occupation of the lattice sites is $`\varphi _c=1`$, and the orientational entropy associated to bonds connecting adjacent segments in the chain is lowest. On the other hand, in the amorphous bulk the fractional occupation $`\varphi _a`$ (i.e. the amorphous-to-crystalline density ratio) is less than $`\varphi _c`$, and is related to an appropriate Flory-Huggins segment/vacuum interaction parameter $`\chi _{p0}`$ (=1.8 for PEO); correspondingly, the bond orientational entropy is highest. In the interfacial region the fractional occupation $`\varphi `$ is higher than $`\varphi _a`$, i.e. the chain segments are somewhat compressed, while a relatively large bond orientational entropy is still maintained. This contributes to the overall $`\sigma _e`$. Now, when a nucleant molecule is put in place of a segment in a site adjoining the folding plane, then i) a nucleant/segment solvation energy arises, which is described by a parameter $`\chi _s`$ ($`0`$), ii) a local segment-segment compression work correspondingly disappears and iii) steric hindrance causes a decrease in the bond orientational entropy associated to the neighboring chain segments. These factors contribute to to the change in $`\sigma _e`$, which can be approximated by 8 :

$$\left(\frac{d\sigma _e}{d\phi }\right)_{n=0}\frac{kT}{b_0^2}\left[\left(\chi _s\chi _{p0}\right)\frac{1}{2}\mathrm{ln}\left(\frac{z}{2}\right)\right],$$
(6)
where $`z`$ $`(=6)`$ is the lattice coordination number; $`\phi `$ is the average number of nucleant molecules diffusing out of the volume of a forming stem, towards one of the ends of the stem itself (for finite $`\mathrm{\Delta }T`$ it is given by $`\phi nb_0^2l_0/2`$); the value $`0.5`$ can be taken for $`\chi _s`$, corresponding to a ”solvation” energy of $`1.3`$ kJ/mol typical of this kind of ”solutions”. Note that in this (rough) approximation $`(\sigma _e/\phi )_{n=0}`$ is independent of $`\sigma _{e0}`$. In fact, the calculation of the full expression doesn’t show a dramatic dependence on $`\sigma _{e0}`$. As a practical example, assuming $`\sigma _{e0}=20`$ erg cm-2 9  one finds $`a1.7`$ K from the expression $`anb_0^2T_mv_cH_f^1(\sigma _e/\phi )_{n=0}`$ 4 ; on the other hand, taking $`H_f/v_c210`$ J cm-3 9 ; 10 , SAXS data imply $`\sigma _{e0}70`$ erg cm-2 (cfr. 10 ), and the calculation yields $`a1.2`$ K in this case. In Table 1 the values of $`a/\mathrm{\Delta }T`$ for both these $`\sigma _{e0}`$ values are reported.

In summary, we have shown that in the present case Eq. 2 is supported by growth kinetics’ data, SAXS and lattice calculations; the latter also providing a microscopic picture of the mechanism.


text:
# References

IMPLICATIONS OF A TIME-VARYING
FINE STRUCTURE CONSTANT

Antonio Alfonso-Faus

*E.U.I.T. Aeronáutica, Spain*

> Much work has been done after the possibility of a fine structure constant being time-varying. It has been taken as an indication of a time-varying speed of light. Here we prove that this is not the case. We prove that the speed of light may or may not vary with time, independently of the fine structure constant being constant or not. Time variations of the speed of light, if present, have to be derived by some other means and not from the fine structure constant. No implications based on the possible variations of the fine structure constant can be imposed on the speed of light.

PACS numbers 04 04.60.-m

The strength of the electromagnetic interaction between photons and electrons is measured by the fine structure constant, $`\alpha =e^2/(4\pi ϵ_0\mathrm{}c)`$. A time variation of this constant has been recently reported as a slow increase over cosmological timescales \[1-4\]. Non-standard cosmological theories have been proposed invoking a varying speed of light \[5-8\] or a varying electronic charge . Using the entropy of a black hole  theoretical reasons have been proposed to favor a varying $`c`$ over a varying $`e`$. However, we have already proved elsewhere  that the speed of light does not enter in the formulation of the fine structure constant. Hence no time variations of the speed of light can be inferred from time variations of the fine structure constant.

All constants of physics are gravitational, with a varying degree of importance, due to the fact that gravitation is present in all the Universe, and nothing can escape that. The Equivalence Principle, and therefore Local Lorentz Invariance (LLI) and/or Local Position invariance (LPI) , need not necessarily be violated by space and/or time variations of certain constants, like the fine structure constant, provided first principles are preserved, as the Action Principle for example. Preservation of Local Lorentz Invariance (LLI) as well as Local Position Invariance (LPI) , imply that it is necessary and sufficient for this preservation that the permittivity $`ϵ_0`$ and permeability $`\mu _0`$ be both equal to $`1/c`$ . Using a TH $`ϵ\mu `$ formalism  one can implement the Einstein s Equivalence Principle and prove that a necessary and sufficient condition for both LLI and LPI to be valid is given by

$$ϵ_0=\mu _0=\left(H_0/T_0\right)^{1/2}$$

for all events. Since the product $`ϵ_0\mu _0`$ is equal to $`c^2`$, then the above relation implies

$$ϵ_0=\mu _0=\frac{1}{c}$$

Hence, in these units we can rephrase Maxwell s equations  giving Coulomb s law as $`F=ce^2/r^2`$ and for the fine structure constant

$$\alpha =\frac{e^2}{4\pi \mathrm{}}$$

so that the speed of light does not explicitly enter in the constitution of the fine structure constant. The arguments using the entropy of a black hole to favor a varying $`c`$ over a varying  $`e`$ are then drastically changed, under LLI and LPI preservation, to favor a varying $`\mathrm{}`$ and/or $`e`$ instead, and not $`c`$.


text:
# The Evolutionary Minority Game with Local Coordination

## INTRODUCTION

There are a great number of situations in which a many agent system self organizes by coordinating individual actions. Such coordination is usually achieved by agents with partial information about the system, and in some cases optimizing utility functions that conflict with each other. A similar situation is found in many particles, physical systems. The word “coordination” used in a social or economic context is then replaced by “ordering”. Examples are the growth of a crystalline structure or a transition leading to some specific magnetic phase.

Interesting situations arise when the optimal configurations for different individuals do collide with each other. This can be due to the nature of the interactions between the particles as in a spin glass, or by boundary conditions which prevent a global ordering or by the constitutive rules of a system of multiple players that prevent that all agents can win. In these cases it is said that the system displays some degree of *frustration*. An example of a frustrated multi agent system is given by the evolutionary minority game (EMG) MGame  in which many players have to make a binary choice and the winning option is the one made by the minority. The similarities between some variants of the Minority Game and spin glasses have been discussed in great detail in Ref. tano .

A macroscopic signature of frustration is that the system can not accommodate into a single, optimal state in which the energy is a minimum but it relaxes instead to one of many, suboptimal, equivalent configurations that are local minima in the energy landscape. In the random relaxation dynamic that is used for the EMG each player continually modifies her choice searching for a winning option. The final result is that the population is segregated into two parties that take opposite actions. This partition is not unique and also tends to reduce the frustration as much as possible by minimizing the number of losers.

The relaxation process is usually assimilated to the search of a solution of a combinatorial optimization problem in which it is possible a strategy of “divide and conquer” libro , i.e. circumstances in which one can attempt to divide the system into parts and search for separate optima in each part. Frustration arises when such local solutions can not be reassembled into a global optimum also fulfilling the boundary conditions.

A relevant example of the study of the global outcome of a local coordination strategy (i.e. involving only a fraction of the system) is Schelling’s segregation model schelling . Agents of two kinds are placed in a square grid. The system relaxes to equilibrium allowing any two agents of different kinds to exchange places if they are surrounded by, say, a majority of agents of the opposite kind. In the present paper we discuss a relaxation dynamics for the EMG in which we impose a local coordination strategy. We borrow the picture of Schelling’s models and place the players in a lattice. It is then possible to associate a neighborhood to each player and thus implement a local coordination strategy letting each player to adjust her decision to the situation in her neighborhood. We call this model the Local Evolutionary Minority Game (LEMG).

A previous work in this direction is suiza  in which players are also located on the nodes of a grid but are endowed with (two) strategies that are selected on the basis of their successful use. This work further imposes that both strategies have to be anticorrelated, i.e. they tend to produce opposite actions with the same input. We stress the *evolutionary* nature of the present model: players bear no memory of past actions and do not have any strategy in the sense of suiza , to guide their actions. In spite of this difference the effects of local coordination produce a similar ordered pattern. We pay special attention to the effects of such a local coordination in the optimization process.

## I The rules of the LEMG

We first consider the traditional EMG. This involves $`N`$ players that make one binary decision (0 or 1). Each player has a probability $`p_i;i=1,2,\mathrm{},N`$ of choosing, say, 0. Each player receives one point if her decision places her in the minority, and loses a point otherwise. When her account of points falls below 0, she changes $`p_ip_i^{}`$ with $`p_i^{}[p_i\mathrm{\Delta }p,p_i+\mathrm{\Delta }p]`$, at random, and $`\mathrm{\Delta }p1`$. Reflective boundary conditions are imposed at $`p_i=0,1`$. All agents are assumed to update the corresponding $`p_i`$’s synchronically. It is customary to display the self-organization of the system through the probability density function $`P(p)`$ obtained in a statistical ensemble of systems that are allowed to relax to equilibrium. The value of $`P(p)dp`$ is the fraction of the population having a probability between $`p`$ and $`p+dp`$ of choosing, say, 0. When the probabilistic relaxation is used, the asymptotic function $`P(p)`$ has a U-shape with two symmetric peaks at $`p0`$ and $`p1`$ thus indicating that the $`N`$ agents have segregated into two parties making opposite decisions. The relaxation process corresponds to the minimization of an “energy” function Thermal  given by the standard deviation $`\sigma `$ defined by:

$$\sigma ^2=\underset{A}{}𝒫(A)(AN/2)^2$$
(1)
where $`𝒫(A)`$ is the probability distribution of parties of $`A`$ agents that have chosen 0. The value of $`\sigma ^2`$ depends upon the properties of $`P(p)`$. In Thermal  it is proven that:

$$\frac{\sigma ^2}{N}=N(\overline{p}1/2)^2+(\overline{p}\overline{p^2})$$
(2)
where $`\overline{p^s}=p^sP(p)𝑑p`$. At equilibrium the linear dependence of $``$ on N disappears, and $`\sigma ^2`$ turns out to be an extensive magnitude proportional to $`N`$. A minimization of $``$ is equivalent to find a distribution $`P(p)`$ with the smallest possible number of losers. In fact $`\sigma ^2`$ is related to the number of losers because

$$\sigma ^2=<(AN/2)^2>=\frac{<(w\mathrm{})^2>}{4}=\frac{<(N2\mathrm{})^2>}{4}$$
(3)
where $`w`$ ($`\mathrm{}`$) is the number of winners (losers) and $`<\mathrm{}>`$ represents an ensemble average.

If one assumes naïvely $`P(p)=\delta (p1/2)`$ corresponding to a symmetric random walk (and thus eliminating the term $`O(N)`$ in Eq.) one gets $`=1/4,`$ while $`P(p)=`$ constant yields $`=1/6`$. A better result is obtained with the usual random relaxation dynamics for the EMG. This yields Thermal  $`1/8`$. Energy and frustration remain linked to each other. For the EMG we can define frustration as $`=\mathrm{}/N`$; which fulfills $`01`$. This definition may also be used for any system involving a game with multiple players. The value $`=0`$ corresponds to a situation such as the “majority game” in which a player is a winner if her decision is the same as the majority. This leads to situations that can be assimilated to a ferromagnetic phase (all the players (spins) have chosen the same option (orientation)). In the EMG there are less winners than losers, and therefore $`1/2<_{EMG}1`$. The lowest possible frustration for the EMG is reached when the $`N`$ (odd) agents are coordinated to produce the largest possible minority, i.e. $`(N1)/2`$. Thus the lowest possible frustration for a finite minority game is $`^{}=(1+1/N)/2`$.

We now turn to the LEMG in which the $`i`$-th player makes her decision depending upon the situation in her neighborhood $`𝒩_i`$. In order to define the neighborhoods we assume three possible spatial orderings. Two of them correspond respectively to a one-dimensional (1D) or a square two- dimensional (2D) regular array with periodic boundary conditions (i.e. respectively a ring and a torus). In the third arrangement the agents are placed in the nodes of a random undirected graph with a fixed number of neighbors for each agent so that a reciprocity relationship is automatically fulfilled (if node $`i`$ is taken to be linked to node $`j`$, the reciprocal is also true). All neighborhoods are assumed to have the same (odd) number $`n`$ of agents (we consider that $`i𝒩_i`$).

The rules of the LEMG are the same as for the EMG except for the important difference that an agent wins or loses points depending whether she is, or she is not, in the minority *of her own neighborhood*. No attention is paid to the agents that do not belong to $`𝒩_i`$. The LEMG coincides with the usual EMG when $`𝒩_i`$ coincides with the complete $`N`$ -agent system. In the regular orderings the neighborhoods are respectively a segment or a square with an odd number of agents. The only agent that updates her $`p_i`$ is located at the center of the square or segment. Notice that an agent may be in the minority (a winner) in her neighborhood and in the majority (a loser) when the entire system is considered, and vice versa.

Let us consider the simple example of an infinite linear chain of agents and a neighborhood with $`n=3`$. We define $`R_i`$ to be the probability that the $`i`$-th agent belongs to the minority of $`𝒩_i`$. We can thus write:

$$R_i=(1p_{i1})p_i(1p_{i+1})+p_{i1}(1p_i)p_{i+1}$$
(4)
where $`\mathrm{}i\mathrm{}`$. The probability that all agents are winners is $`R=_iR_i`$. Obviously $`R=1`$ if and only if $`R_i=1,i`$. This is possible only if $`p_i=1`$ and $`p_{i\pm 1}=0`$. This corresponds to a pattern in which 0’s and 1’s alternate with a period of 2. Larger neighborhoods give also rise to periodic solutions with larger (even) periods. Any finite ring with an odd number of agents is frustrated because such periodicity can not fit along the chain.

Equation  can be used to construct a (deterministic) relaxation dynamics to adjust the $`p_i`$’s climbing along the gradient of $`R_i`$. We thus assume $`p_i(t+1)=p_i(t)+\mathrm{\Delta }p_i`$ and set

$$\mathrm{\Delta }p_i=\eta R_i/p_i=\eta (1p_{i1}p_{i+1})$$
(5)
with $`1\eta >0`$. A stationary ($`\mathrm{\Delta }p_i=0i`$) solution for this dynamics is $`p_i=1/2;i`$. This solution is unstable because any random perturbation of any $`p_i`$ leads to a situation in which $`\mathrm{\Delta }p_i0i`$ foot2 . This dynamics stabilizes a pattern of 0’s and 1’s that alternate with each other. In fact if $`p_{i+1}`$ and $`p_{i1}`$ are both greater (smaller) than 1/2, then $`\mathrm{\Delta }p_i<0(>0)`$ thus forcing $`p_i<1/2(>1/2)`$. This relaxation dynamics is therefore expected to lead to distributions $`P(p)`$ that vanish at $`p=1/2`$.

## II RESULTS

### Self-segregation

In Fig.1 we show the results of $`=\sigma ^2/N`$ as a function of the size parameter $`S`$ defined as the ratio $`S=n/N`$ obtained in several numerical experiments. The value for $`S=1`$ corresponds to $`_{EMG}`$ obtained for the EMG. In this section we only discuss results for $`N=121`$. We have considered the topologies of a ring, a “square” torus (with $`N=11\times 11`$), and of random graphs. In all the cases considered, the values of $`_S`$ with $`S1`$ fulfill $`_S<_{EMG}`$. This feature is stressed in Fig.1 with an horizontal line drawn at the value of $`_{EMG}`$. The value of $`_S`$ for regular 1D and 2D lattices grows with $`S`$ and for $`S.5`$ becomes even larger that $`_{EMG}`$.

The corresponding density distributions $`P(p)`$ are shown in Fig.2 for several values of $`S`$. These are compared with the distribution $`P_{EMG}(p).`$ We observe that $`P_S(p)`$ with $`S1`$ are always symmetric and U-shaped as $`P_{EMG}(p)`$ but they differ from this in the fact that they vanish around $`p=1/2`$. This agrees with the discussion given above for the linear chain. As we shortly discuss, this turns out to be a highly relevant and general feature of the LEMG. Such distributions are a better approximation to an ideal distribution

$$P(p)=\frac{N1}{2N}[\delta (p)+\delta (p1)]+\frac{1}{N}\delta (p1/2)$$
(6)
that yields the optimal value of $`_{EMG}^{}=1/4N`$ (and $`_{EMG}^{}=(1+1/N)/2`$). A noticeable dip is produced for $`S=(N2)/N`$. This can be understood in the following way. Assume that a symmetric distribution $`P(p)`$ has already developed and two agents are removed in order that the $`i`$th player can check her decision in her neighborhood of $`N2`$ agents. If the two agents that have been removed have $`p>1/2`$ ($`p<1/2`$) the $`i`$th agent has the single winning option of choosing $`p_i1`$ ($`p_i0`$). In the other cases (one player with $`p>1/2`$ and the other with $`p<1/2`$) her choices of approaching 0 or 1 have equal probability. The net result once all players have updated her respective $`p`$’s in the same fashion is to force the resultant $`P(p)`$ to drop at p=1/2 and grow at $`p=0,1`$ as discussed for the linear chain. This argument can be extended for neighborhoods of other sizes.

### The Optimization Problem

The shape of $`P_S(p)`$ changes with $`S`$. For intermediate values (for instance $`S=81/121`$) this distribution has radically changed from the shape of a U to a two wing profile with secondary maxima at both sides of $`p=`$1/2 still keeping the fact that $`P_S(1/2)=0`$. This is associated to an increase in the number of “local winners” . In fact after some time there are left almost no players that need to update their $`p_i`$’s (see Fig.3) while for $`S0`$ or $`S1`$ the relaxation process reaches a dynamical equilibrium in which few players continuously update their $`p_i`$’s.

For such intermediate values of $`S`$ it is found that a minority is clearly defined in most neighborhoods and the corresponding agents are unambiguously induced to take one winning option. They therefore continue to accumulate points and cease to change their $`p_i`$’s preventing the system to reach a more efficient self-segregation.

A situation like this has been extensively discussed in Ref.Quenching . In these circumstances the relaxation process ceases to be effective to lower the energy and the system freezes in a configuration that is far from a better *local* optimum. Although the frozen microstates do depend upon initial conditions, the asymptotic density $`P(p)`$ is independent both because this is a density function that is associated with a macrostate of the system and because possible random variations are averaged out by repeating the relaxation process for a large ensemble of systems. Moreover, succesive runs to obtain $``$ changing the random initial configuration, yield essentially the same result. This is indeed to be expected because the energy is defined as the average of a statistical fluctuation, as given in Eq..

The procedure to regain the true optimal self-segregation pattern, is to force the relaxation procedure by periodically removing the points that have been accumulated by every player, resetting their accounts to 0. This procedure changes *only* the situations in which the system is frozen but leaves unchanged situations in which this does not happen such as for instance for $`S0`$ or $`S1`$.

These “annealing” episodes melt the system thus making it possible to reach the best local configurations. We have performed an annealed relaxation (with a fixed annealing protocol) for all three topologies. The results are displayed with open symbols in all three panels of Fig.1. The fluctuations in the lower curve of panel a) for $`S0.6`$ disappear for $`N500`$ thus indicating that it is a finite size artifact of the model.

A remarkable result displayed in Fig.1 is that the composition of local optima always yields a better coordination than the one obtained within the framework of the EMG in which all agents are involved in the same relaxation process.

There are actually two situations to consider. One in which $`S0`$ or $`S1`$, and the other where $`S`$ is within these two extreme values. In the first case the LEMG always yields remarkably lower values of $``$. This is indeed a general result, because holds true no matter the topology and the size of the system in which the players are located. In these cases no annealing is required because the system never gets quenched.

Outside the neighbourhood of $`S0`$ or $`S1`$, the system gets quenched, and for regular topologies a poorer value of $``$ (*i.e.* greater than $`_{EMG}`$) is obtained. However, for such values a better (lower) value of $``$ can always be obtained if the annealing procedure is used to force the relaxation process beyond the frozen states. This actually means that a better coordination, i.e. a lower frustration, is achieved whenever the system is found in a configuration that corresponds to a composition of good local optima, and this is true in the whole range of values of $`S`$. If good local optima are guaranteed one finds values of $`_S`$ that are significantly lower that $`_{EMG}`$. A typical value of $`_S`$ 1/16 is obtained in this way that is half the value $`_{EMG}`$ 1/8.

Except for finite size effects, all the results presented in this subsection do not differ from those obtained for many other values of $`N`$ that we have investigated.

### The 2D Case

A much richer situation is found for the case in which players are located on a grid with the topology of 2D torus. An example is shown in Fig.4 in which $`N=31\times 31`$. The values of $`p`$ are associated to shades of gray. Frustration can be perceived in the fact that there is not a single global ordering of black ( $`p0`$) and white ($`p1`$) stripes for the whole array. These are instead grouped in domains with different orientations or with the same orientation but shifted with respect to each other. The relaxation process is fast in an initial stage and slows down once the domains have fully developed. The domain walls are a source of frustration. In fact, when such stage has been reached all the agents that have 0 points and continue to update their $`p_i`$’s are located in the domain walls giving rise to a slow dynamics in which walls move enlarging or shrinking domains. The whole picture resembles a crystallization process; for a value of $`S`$ grater than a critical threshold, all domains collapse into a single pattern of stripes. Frustration shows up as an indented (fuzzy) border in some of the stripes. These results will be discussed in detail elsewhere.

## III Conclusions

We have studied the organization pattern achieved by many agents playing an EMG with local coordination. We find important differences between the coordination achieved when the whole ensemble of agents participates in the same relaxation process or when local coordination is imposed.

According to the present results, the LEMG is an example in which a better coordination or, what is the same, a lower frustration, can be achieved provided that the self organization is ruled by a local process, i.e. if agents govern their actions paying no attention to global trends of the system but rather to her immediate neighborhood. In the LEMG, this statement holds true even in the case in which very few agents are removed ($`S(N2)/N`$ ) from the whole ensemble. The fact that a more efficient coordination is achieved by *ignoring* what happens to the total ensemble of players is expected to be a feature of a special kind of multi-agent systems. Other coordination problems may not behave in this fashion, thus allowing a classification of coordination games into classes linked to the type of optimization problem that is being “solved” by the ensemble of agents. Further investigation should be devoted to classify multi agent games into those that fulfil this property and those that can not be optimized by braking them into pieces.

It is found that the LEMG displays some of the features that are typical of antiferromagnetic systems including the emergence of domains, frustration, fast and slow dynamics, etc., while also keeping the essence of multiagent models, as applied to social or economic organization. It therefore sheds light on the connection between those two bodies of knowledge. E.B. has been partially supported by CONICET of Argentina, PICT-PMT0051.