text:
# Lattice Boltzmann Study of Velocity Behaviour in Binary Mixtures Under Shear

## I Introduction

Recently, the dynamical behaviour of complex fluid and binary mixtures under the action of applied flows has been investigated in many theoretical and experimental works . In phase separating binary mixtures the morphology and the growth properties of the domains of the two components are greatly affected by the presence of a flow. This flow is generally imposed from the boundaries onto the system and takes some time to reach its steady value. If this time is comparable with other typical time-scales in the system, transient effects have to be carefully considered. For example, in the case of an oscillating shear flow, the ratio between the period of the applied flow and the relaxation time of a triangular steady shear profile is relevant for the effectiveness of the flow on the bulk of the system. In spinodal decomposition, if that ratio is small, the growth remains mainly isotropic and elongated domains can be observed only in the boundary layers .

The estimation of relaxation times of applied flows is generally based on analogies with simple fluids. However, a careful checking of the laws valid for simple fluids in cases when mesoscopic structures are present in the system - e.g. interfaces between domains in spinodal decomposition, cannot be found in literature. In this paper we evaluate numerically the relaxational properties of a triangular shear flow applied on phase separating binary mixtures. At a certain time during phase separation the top and the bottom boundaries of our system start to move with velocities $`w`$ and $`w`$ respectively and after a while a triangular profile will appear. We study the evolution of the velocity field in terms of the viscosity and of the size of the system and find that the results for a simple fluid also hold in the case of binary mixtures with interfaces inside.

We simulate the binary mixture system by applying lattice Boltzmann methods which allow to mimic the behaviour of a system described by Navier-Stokes and convection-diffusion equations. The lattice Boltzmann algorithm is based on a collision step and a propagation step occurring on time $`\mathrm{\Delta }t`$. Only a finite sets of local velocities are allowed. In particular we use a scheme based on a free-energy approach which has the advantage that the equilibrium thermodynamics of the system is “a priori” known. The role of the simulation time step $`\mathrm{\Delta }t`$ will be also considered. With fast relaxing velocity profiles, the usual choice $`\mathrm{\Delta }t=1`$ is not appropriate for obtaining smoothly relaxing profiles and smaller $`\mathrm{\Delta }t`$ have to be considered.

In the next section we briefly describe the methods used in the simulations; section III contains our results and section IV some conclusions.

## II The model

We consider a two-dimensional binary fluid with components A and B of number density $`\rho _A`$ and $`\rho _B`$, respectively. Such a system can be modeled by the following free energy functional,

$$F=𝑑𝐫[\frac{1}{3}\rho \mathrm{ln}\rho +\frac{a}{2}\phi ^2+\frac{b}{4}\phi ^4+\frac{\kappa }{2}(\phi )^2]\text{}$$
(1)
where $`\rho =\rho _A+\rho _B`$ is the local total density and $`\phi =\rho _A\rho _B`$ is the local density difference and the order parameter of the system; the term in $`\rho `$ gives rise to a positive background pressure and does not affect the phase behavior. The terms in $`\phi `$ correspond to the usual Ginzburg-Landau free energy typically used in studies of phase separation. The polynomial terms are related to the bulk properties of the fluid. The gradient term is related to the interfacial properties. The parameter $`b`$ is always positive, while the sign of $`a`$ distinguishes between a disordered ($`a>0`$) and a segregated mixture ($`a<0`$) where the two pure phases with $`\phi =\pm \sqrt{a/b}`$ coexist. In this paper we will consider quenches into the coexistence region with $`a<0`$ and $`b=a`$, so the equilibrium values for the order parameter are $`\phi =\pm 1`$. The initial state in simulations will be random configurations corresponding to the high temperature disordered phase.

The thermodynamic properties of the fluid follow directly from the free energy (1). The chemical potential difference between the two fluids is given by

$$\mathrm{\Delta }\mu =\frac{\delta F}{\delta \phi }=a\phi +b\phi ^3\kappa ^2\phi \text{}$$
(2)
The pressure is a tensor $`P_{\alpha \beta }`$ since interfaces in the fluids can exert non-isotropic forces. A suitable choice is

$$P_{\alpha \beta }=p_0\delta _{\alpha \beta }+\kappa _\alpha \phi _\beta \phi \text{,}$$
(3)
where the diagonal part $`p_0`$ can be calculated using thermodynamics relations from (1):

$`p_0`$ $`=`$ $`\rho {\displaystyle \frac{\delta F}{\delta n}}+\phi {\displaystyle \frac{\delta F}{\delta \phi }}f(\rho ,\phi )`$ (4)
$`=`$ $`{\displaystyle \frac{\rho }{3}}+{\displaystyle \frac{a}{2}}\phi ^2+{\displaystyle \frac{3b}{4}}\phi ^4\kappa \phi (^2\phi ){\displaystyle \frac{\kappa }{2}}(\phi )^2`$ (5)
Our simulations are based on the lattice Boltzmann scheme developed by Orlandini et al. and Swift et al. . In this scheme the equilibrium properties of the system can be controlled by introducing a free energy which enters properly into the lattice Boltzmann model. The scheme used in this paper is based on the D2Q9 lattice: A square lattice is used in which each site is connected with nearest and next-nearest neighbors. The horizontal and vertical links have length $`\mathrm{\Delta }x`$, the diagonal links $`\sqrt{2}\mathrm{\Delta }x`$, where $`\mathrm{\Delta }x`$ is the space step. Two sets of distribution functions $`f_i(𝐫,t)`$ and $`g_i(𝐫,t)`$ are defined on each lattice site $`𝐫`$ at each time $`t`$. Each of them is associated with a velocity vector $`𝐞_i`$. Defined $`\mathrm{\Delta }t`$ as the simulation time step, the quantities $`𝐞_i\mathrm{\Delta }t`$ are constrained to be lattice vectors so that $`|𝐞_i|=\mathrm{\Delta }x/\mathrm{\Delta }tc`$ for $`i=1`$, $`2`$, $`3`$, $`4`$ and $`|𝐞_i|=\sqrt{2}c`$ for $`i=5`$, $`6`$, $`7`$, $`8`$. Two functions $`f_0(𝐫,t)`$ and $`g_0(𝐫,t)`$, corresponding to the distribution components that do not propagate ($`𝐞_0=0`$), are also taken into account. The distribution functions evolve during the time step $`\mathrm{\Delta }t`$ according to a single relaxation time Boltzmann equation:

$$f_i(𝐫+𝐞_i\mathrm{\Delta }t,t+\mathrm{\Delta }t)f_i(𝐫,t)=\frac{1}{\tau }[f_i(𝐫,t)f_i^{eq}(𝐫,t)],$$
(6)
$$g_i(𝐫+𝐞_i\mathrm{\Delta }t,t+\mathrm{\Delta }t)g_i(𝐫,t)=\frac{1}{\tau _\phi }[g_i(𝐫,t)g_i^{eq}(𝐫,t)],$$
(7)
where $`\tau `$ and $`\tau _\phi `$ are independent relaxation parameters, $`f_i^{eq}(𝐫,t)`$ and $`g_i^{eq}(𝐫,t)`$ are local equilibrium distribution functions. Following the standard lattice Boltzmann prescription, the local equilibrium distribution functions can be expressed as an expansion at the second order in the velocity $`𝐯`$:

$$\begin{array}{ccc}f_0^{eq}\hfill & =\hfill & A_0+C_0v^2\hfill \\ f_i^{eq}\hfill & =\hfill & A_I+B_Iv_\alpha e_{i\alpha }+C_Iv^2+D_Iv_\alpha v_\beta e_{i\alpha }e_{i\beta }+G_{I,\alpha \beta }e_{i\alpha }e_{i\beta }\hfill \\ & & i=1,2,3,4,\hfill \\ & & \\ f_i^{eq}\hfill & =\hfill & A_{II}+B_{II}v_\alpha e_{i\alpha }+C_{II}v^2+D_{II}v_\alpha v_\beta e_{i\alpha }e_{i\beta }+G_{II,\alpha \beta }e_{i\alpha }e_{i\beta }\text{,}\hfill \\ & & i=5,6,7,8\text{}\hfill \\ & & \text{ }\hfill \end{array}$$
(8)
and similarly for the $`g_i^{eq}`$, $`i=0`$, $`\mathrm{}`$, $`8`$. The expansion coefficients $`A_0`$, $`A_I`$, $`A_{II}`$, $`B_I`$, $`\mathrm{}`$ are determined by using the following relations

$$\begin{array}{cccc}_if_i^{eq}(𝐫\text{}t)=\rho \text{,}\hfill & & & _ig_i^{eq}(𝐫\text{}t)=\phi \text{,}\hfill \end{array}$$
(9)
$$\begin{array}{cccc}_if_i^{eq}(𝐫\text{}t)𝐞_i=\rho 𝐯\text{,}\hfill & & & _ig_i^{eq}e_{i\alpha }=\phi v_\alpha \text{}\hfill \end{array}$$
(10)
$$\begin{array}{cccc}_if_i^{eq}e_{i\alpha }e_{i\beta }=c^2P_{\alpha \beta }+\rho v_\alpha v_\beta \text{}\hfill & & & _ig_i^{eq}e_{i\alpha }e_{i\beta }=c^2\mathrm{\Gamma }\mathrm{\Delta }\mu \delta _{\alpha \beta }+\phi v_\alpha v_\beta \text{,}\hfill \end{array}$$
(11)
where $`P_{\alpha \beta }`$ is the pressure tensor, $`\mathrm{\Delta }\mu `$ is the chemical potential difference between the two fluids and $`\mathrm{\Gamma }`$ is a coefficient related to the mobility of the fluid. We stress that we are considering a mixture with two fluids having the same mechanical properties and, in particular, the same viscosity. The second constraint in Eq. (10) expresses the fact that the two fluids have the same velocity.

A suitable choice of the coefficients in the expansions (8) is shown in Ref.. Such a lattice Boltzmann scheme simulates at second order in $`\mathrm{\Delta }t`$ the continuity, the quasi-incompressible Navier-Stokes and the convection-diffusion equations with the kinematic viscosity $`\nu `$ and the macroscopic mobility $`\mathrm{\Theta }`$ given by

$$\begin{array}{cccc}\nu =\mathrm{\Delta }t\frac{c^2}{3}(\tau \frac{1}{2})\text{}\hfill & & & \mathrm{\Theta }=\mathrm{\Gamma }\mathrm{\Delta }tc^2(\tau _\phi \frac{1}{2})\text{}\hfill \end{array}$$
(12)
The shear flow can be imposed by introducing boundary walls on the top and bottom rows of lattice sites. The velocities of the upper and lower walls are along the horizontal direction and their values are $`w`$ and $`w`$, respectively, where $`w=\gamma (L1)/2`$ and $`\gamma `$ is the shear rate imposed on the system. The bounce-back rule is adopted for the distribution functions normal to the boundary walls. In order to preserve correctly mass conservation, a further constraint, related to the distribution functions at the previous time step $`(t\mathrm{\Delta }t)`$, is used. Details of the scheme are given in Ref. .

Finally, we observe that if we set $`a=b=\kappa =0`$ in the free energy functional (1), the present lattice Boltzmann methods can be used to simulate simple fluids.

## III Results

We numerically check the relaxation behaviors of the horizontal velocity profile for symmetric and asymmetric binary fluids. We consider two cases: (i) switch on the shear from the beginning of the phase separation and (ii) switch on the shear during the phase separating process and after the interfaces between the two fluids have formed. We denote the time at which the shear is switched on as $`t_{on}`$. We focus on the effects of viscosity, so we vary $`\tau `$ in the simulations. The other parameters, if not differently stated, are fixed at the values: $`a=b=1.252\times 10^4`$, $`\kappa =8\times 10^5`$, $`\mathrm{\Gamma }=40`$, $`\mathrm{\Delta }x=1`$, $`\mathrm{\Delta }t=0.2`$, $`\gamma =0.005`$. Similar results are obtained for other sets of parameters.

Figure 1 shows the relaxation of the horizontal velocity profile for a symmetric binary fluid with $`\tau =2`$, $`L=64`$, $`t_{on}=5000`$. The four lines correspond to $`t=5020`$, $`t=5040`$, $`t=5060`$, $`t=5120`$. From this figure we can see how the shear flow comes into the bulk of the system. When $`t>5120`$ the velocity profile shows an almost linear behavior with $`y`$. Before switching on the shear the domains grow isotropically. After $`t_{on}`$ anisotropic behavior come into the bulk of the system with time. Figure 2 shows the configuration at the time $`t=5120`$. Domains separated by well formed interfaces can be observed to incline to the flow direction with time. For the case $`t_{on}=0`$ we observed the same behavior for the velocity profile, while the interfaces have not reached their equilibrium shape. That means the existence of the interfaces does not influence the shear effects coming into the system. For the asymmetric case, we studied the binary fluid with the ratio $`80:20`$. After switching on the shear the velocity profile shows the same behavior.

To understand better the time evolution of the velocity profile of binary mixtures, we take the simple fluid as a reference. We consider the Newtonian viscous flow between two infinite plates with a distance of $`L`$. We use $`y`$ to denote the coordinate in the vertical direction and $`y=0`$ in the middle of the system. The velocities of the upper and the lower plates are $`w`$ and $`w`$, respectively . The shear rate imposed on the system is $`\gamma `$. The motion equation is

$$\rho \frac{u}{t}=\eta \frac{^2u}{y^2}\text{ }$$
(13)
The velocity profile can be obtained by standard methods and is given by

$$u=\gamma y\underset{n}{}(1)^{n+1}\frac{\gamma L}{n\pi }\mathrm{exp}(\frac{4n^2\pi ^2\nu }{L^2}t)\mathrm{sin}(\frac{2n\pi }{L}y)\text{.}$$
(14)
When $`t`$ is large enough, the modes with $`n2`$ can be neglected, which is confirmed by simulations. So we can define a relaxation time for the velocity profile in the following way,

$$T_R=\frac{L^2}{4\pi ^2\nu }=\frac{3L^2}{2\pi ^2\mathrm{\Delta }tc^2(2\tau 1)}\text{.}$$
(15)
It is interesting to see if or not such a definition also works for binary fluids. To numerically check the relation between $`T_R`$ and $`\nu `$ (or $`\tau `$) or $`L`$, we calculate $`T_R`$ in the following way,

$$T_R=\frac{t}{\mathrm{ln}(\mathrm{\Xi })}\text{}$$
(16)
where

$$\mathrm{\Xi }=\frac{\pi (\gamma yu)}{\gamma L\mathrm{sin}(\frac{2\pi }{L}y)}\text{}$$
(17)
In order to calculate $`T_R`$, we use a time $`t`$ at which the velocity profile has been almost linear, so that terms with $`n2`$ in Eq.(14) can be neglected.

Figure 3 shows the simulation results for $`T_R`$ as a function of $`\nu `$, where $`L=128`$, $`\mathrm{\Delta }t=0.2`$, the vertical axis is for $`1/T_R`$ and the horizontal axis is for $`\nu `$ so that we expect a linear behavior. The dotted line with points correspond to the simulation results and the solid line corresponds to the expected value from the definition (15). The simulation results confirm the existence of the exponential behavior in the relaxation process of the velocity profile and the validity of the definition of $`T_R`$.

Figure 4 shows $`T_R`$ as a function of $`L`$, where the vertical axis is for $`T_R`$ and the horizontal axis is for $`L^2`$ . The dotted lines with symbols are simulation results for the cases of $`\nu =4.17`$, $`\nu =5.83`$ and $`\nu =7.5`$. The solid lines are the expected values from the definition (15). The four points in each case correspond to $`L=64`$, $`128`$, $`256`$, and $`512`$. The simulation step in this figure is $`\mathrm{\Delta }t=0.2`$. We find the expected linear behavior between $`T_R`$ and $`L^2`$.

There is a bending tendency in the simulation results for $`1/T_R`$ in Fig.3, so that the simulation results deviate from the expected values. If we continue the simulation up to a much higher viscosity region, we will find a more pronounced bending behavior. We emphasize that this is an artificial phenomenon, which can be mainly attributed to the finite size of the simulation step $`\mathrm{\Delta }t`$. If the viscosity is larger the velocity profile relaxes more quickly and we should use a smaller time step $`\mathrm{\Delta }t`$ for the simulation. When the viscosity is very high and $`\mathrm{\Delta }t`$ is not small enough we can not observe smooth velocity profiles in the simulation. To confirm this numerical analysis we show two simulation results in Fig.5, where the definition value of $`1/T_R`$ is also shown to guide the eyes. Compared with the case of $`\mathrm{\Delta }t=0.2`$, when we use $`\mathrm{\Delta }t=0.1`$ the linear behavior region is almost doubled.

Finally, it is interesting to compare the relaxation behaviors of binary and simple fluids. Fig.6 shows the simulation results of the velocity profiles for the case with $`L=256`$, $`\tau =2`$ ($`\nu =2.5`$), $`tt_{on}=600`$. The solid line corresponds to Eq. (14) with $`n=1`$. From this figure the following remarks are evident: (i) the relaxation process of the binary mixtures follows the same behavior as that of simple fluids; (ii)for binary mixtures the relaxation behavior of the velocity profile is independent of the component ratio and the time when shear is switched on, which means that the forming of the interfaces does not evidently influence the shear effects coming into the system.

## IV Conclusions

In this paper we study the relaxation behavior of the velocity profile of binary mixtures under steady shear with lattice Boltzmann methods. Following the simple Newtonian viscous flow, we define a relaxation time of the velocity profile $`T_R`$ whose value is inversely proportional to the viscosity and proportional to the square of the size of the system. The simulation results show that the shear effects come into the system in the same way for binary mixtures and simple fluids, which confirm the validity of the definition of $`T_R`$ in this and previous studies. For binary mixtures, the shear behavior is independent of the component ratio and the time at which the shear is switched on. The presence of interfaces between the two fluids has negligible influence on the relaxation process of the velocity profile.


format: text

text:
# Secular variability, geodetic precession and moment of inertia of binary pulsars

## I Introduction

In the first gravitational two-body equation including spinsbo , the precession velocity of the orbital angular momentum vector, $`𝐋`$, was expressed as around a vector, combined by the spin angular momenta of the pulsar, $`𝐒_1`$, and its companion $`𝐒_2`$. This velocity has been ignored, since it is much smaller, 2 Post-Newtonian order (PPN), than the precession velocity of the pulsar spin axis, 1.5PPN in binary pulsar systems.

However, the small orbital precession velocity relative to the vector combined by $`𝐒_1`$ and $`𝐒_2`$ doesn’t mean that it is also small relative to the line of sight (LOS). Then in the study of observational effect of the orbital precession due to the coupling of the spin induced quadrupole moment with the orbital angular momentum sb ; lai ; wex2 , the precession of $`𝐋`$ is expressed as relative the total angular momentum vector, $`𝐉`$, which is at rest to LOS (counting out the proper motion effect); rather than the vector ($`𝐒_1`$ and $`𝐒_2`$), which are precessing relative to the LOS rapidly themselves.

Thus, there seems a contradictory between Barker and O’Connell’s equation and that of other authors mention above on which direction $`𝐋`$ should precess around. Actually they can be consistent in the following scenariosb ; hs . In which $`𝐋`$, $`𝐒_1`$ and $`𝐒_2`$ all precess around $`𝐉`$ rapidly (1.5PPN), while the relative velocities of $`𝐋`$ relative to $`𝐒_1`$ and $`𝐒_2`$ are very small (2PPN).

The coupling of the spin induced quadrupole moment and the orbit is 2 PPN bo , so the quadrupole moment of the companion star has to be very large to make the effect observable. Therefore, this model is suitable for such binary pulsars, as neutron star-main sequence star (NS-MS), and neutron star-black hole (NS-BH)lai ; wex2 .

Unlike NS-MS and NS-BH binaries, neutron star-white dwarf (NS-WD) and neutron star-neutron star (NS-NS) binaries have much smaller quadrupole moment, and the spin angular momenta of the two stars can be comparable. To explain the observations of such binaries, alternative model is necessary.

Apostolatos et alapo  and Kidderkid  studied the precession of $`𝐋`$ as relative to $`𝐉`$ in discussing the modulation of gravitational wave by the spin-orbit (S-L) coupling in coalescing binary systems (actually NS-BH). Similar to the former quadrupole-orbit coupling, this S-L coupling only considers the spin of the companion star with the orbit, and the spin of the pulsar is ignored. The difference is that the former quadrupole-orbit coupling corresponds to 2 PPN, whereas the latter S-L coupling corresponds to 1.5 PPN. Therefore, the latter is easier to explain the measurements of NS-WD and NS-NS binaries, which have smaller quadrupole moments and comparable spins.

This paper develops the equation of orbital precession of Apostolatos et al and Kidder, making it applicable for general binary pulsars. We point out that once a binary system has two spins instead of one (the masses of the binary system can be different), the precession of the spin angular momenta of the two stars, $`𝐒_1`$ and $`𝐒_2`$, lead to a variable $`𝐒`$ ($`𝐒𝐒_1+𝐒_2`$), which in turn tilts the orbital plane. Thus the additional motion of the orbital plane can explain both $`\dot{x}`$ and $`\dot{P}_b`$ measured in NS-WDs, i.e., PSR B1957$`+`$20aft  and PSR J2051$``$0827do , which have been interpreted separately by different modelswex2 ; as . Moreover, the derivatives, $`\ddot{P}_b`$, $`\stackrel{\dot{}\dot{}\dot{}}{P}_b`$, and $`\ddot{x}`$ can also be naturally interpreted.

The geodetic precession induced $`\dot{x}`$ and $`\dot{P}_b`$ include an angle, $`\lambda _{LJ}`$ (the angle between $`𝐋`$ and $`𝐉`$), which represents the intensity of S-L coupling ($`\lambda _{LJ}S/L`$). Comparing with the measured $`\dot{x}^{obs}`$, the magnitude of S/L can be obtained. Since L can be obtained from the binary parameters measured, we can thus estimate the magnitude of $`𝐒_1`$ and $`𝐒_2`$, and in turn the moment of inertia of the pulsar, through the measured pulsar period. The obtained moment of inertia is about $`10^{44}`$g cm2 and $`10^{45}`$g cm2, which is in the range of strange stars and NSs.

On the other hand, the geodetic precession model is supported not only by the secular variabilities measured, but also by theoretical predictions of the moment of inertia of pulsars.

Section II introduces the orbital precession in special cases. Section III discusses the orbital precession in general cases and its relationship with the orbital precession in special cases. Section IV applies the general orbital precession to two NS-WD binaries, explains the significant secular variabilities measured and estimates the moment of inertia of pulsars in the two binaries.

## II orbital precession in special cases

The motion of a binary system can be seen as the precession of three vectors, the spin angular momenta of the pulsar and its companion star, $`𝐒_1`$ and $`𝐒_2`$, and the orbital angular momentum $`𝐋`$. The change of the orbital period due to the gravitational radiation is 2.5 PPN; whereas the geodetic precession corresponds to 1.5 PPN. So the influence of gravitational radiation on the motion of a binary system can be ignored in the discussion of dynamics of a binary pulsar system. Therefore, the total angular momentum, $`𝐉=𝐋+𝐒_1+𝐒_2`$, can be treated as invariable both in magnitude and direction ($`\dot{𝐉}=0`$). With $`\mathrm{\Omega }_0`$, the precession rate of $`𝐋`$ around $`𝐉`$, the S-L coupling can be expressed asbo

$$𝛀_0\times 𝐋=𝛀_1\times 𝐒_1𝛀_2\times 𝐒_2,$$
(1)
where $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ represent the precessions of the pulsar and its companion star, respectively. Ignoring those terms that are over 2 PPN, $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ can be written asbo

$$\mathrm{\Omega }_1=\frac{L}{2r^3}(4+\frac{3m_2}{m_1}),\mathrm{\Omega }_2=\frac{L}{2r^3}(4+\frac{3m_1}{m_2}),$$
(2)
where $`m_1`$ and $`m_2`$ are masses of the pulsar and the companion star, respectively, and $`r`$ is the separation of $`m_1`$ and $`m_2`$. Notice $`Lr^{1/2}`$, $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ are 1.5 PPN.

In the case of one spinning body, i.e., $`𝐒_2=0`$, $`𝐒=𝐒_1`$, $`𝐋`$ and $`𝐒`$ precess about the fixed vector $`𝐉`$ at the same rate with a precession frequency approximatelykid ,

$$\omega _p=\frac{|𝐉|}{2r^3}(4+\frac{3m_2}{m_1}).$$
(3)
Eq(3) is also correct if the two bodies have equal masses ($`𝛀_1=𝛀_2`$ in Eq(1)).

Section III shows that the orbital precession in a general binary pulsar system ($`S_20`$, $`S_20`$, $`m_1m_2`$) will automatically lead to the significant variability in $`\dot{P}_b`$ and $`\dot{x}`$.

## III orbital precession and its effects in general cases

Eq(1) indicts that at any instance the variation of the angular momenta $`𝐋`$ (left) equals the variation of the angular momenta $`𝐒_1`$ and $`𝐒_2`$ (right). In magnitude, the right-hand side of Eq(1) is (recall $`𝐒𝐒_1+𝐒_2`$)

$$(\tau _1)_R=\mathrm{\Omega }_2S\mathrm{sin}\lambda _{LS}+(\mathrm{\Omega }_1\mathrm{\Omega }_2)S_1^{}\mathrm{sin}\lambda _{LS_1}.$$
(4)
$`𝐋`$ can react to the torque, $`(\tau _1)_\mathrm{R}`$, by precessing around $`𝐉`$, with a very small opening angle of the precession cone $`\lambda _{LJ}S/L`$. Since the triangle ($`𝐋`$, $`𝐉`$ and $`𝐒`$) constraint much be satisfied at any instant. Thus the left-hand side of Eq(1) can be written as

$$(\tau _1)_L=|\mathrm{\Omega }_0𝐧_J\times L𝐧_0|=\mathrm{\Omega }_0L\mathrm{sin}\lambda _{LJ},$$
(5)
where $`𝐧_J`$ and $`𝐧_0`$ are unit vectors of $`𝐉`$ and $`𝐋`$, respectively. By Eq(4) and Eq(5), we have the precession rate of $`𝐋`$ around $`𝐉`$,

$$\mathrm{\Omega }_0=\mathrm{\Omega }_2\mathrm{sin}\lambda _{LS}+(\mathrm{\Omega }_1\mathrm{\Omega }_2)\frac{S_1^{}}{S}\mathrm{sin}\lambda _{LS_1},$$
(6)
where $`S_1^{}=S_1\mathrm{cos}\eta _{SS1}`$, which denotes the component of $`𝐒_\mathrm{𝟏}`$ in the plane determined by $`𝐒`$ and $`𝐉`$. Note that $`L\mathrm{sin}\lambda _{LJ}S`$ is used in Eq(6), since $`S/L1`$. The difference between Eq(3) and Eq(6) is obviously due to different number of spins included. Notice that the right-hand side of Eq(6) can as well be written by replacing subscribes 1 with 2 and 2 with 1.

Thus the orbital precession velocity of Eq(6), $`\mathrm{\Omega }_0`$, is derived via two assumptions: the conservation of the total angular momentum; and the triangle constraint on the precession cone of the orbit, $`\lambda _{LJ}`$.

Notice that $`\mathrm{\Omega }_0`$ of Eq(6) is 1.5 PPN. Which can be absorbed by the advance of the precession of the periastron, $`\dot{\omega }`$. The measured $`\dot{\omega }`$ is givensb ; wex2  by:

$$\dot{\omega }^{obs}=\dot{\omega }^{GR}+\mathrm{\Omega }_0\mathrm{cos}\lambda _{LJ}\dot{\omega }^{GR}+\mathrm{\Omega }_0.$$
(7)
Since $`𝐒_\mathrm{𝟏}`$ and $`𝐒_\mathrm{𝟐}`$ precess with different velocities, $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ respectively ($`m_1m_2`$), $`𝐒`$ varies both in magnitude and direction ($`𝐒_\mathrm{𝟏}`$, $`𝐒_\mathrm{𝟐}`$ and $`𝐒`$ form a triangle). From the triangle of $`𝐒`$, $`𝐋`$ and $`𝐉`$, in reaction to the variation of $`𝐒`$, $`𝐋`$ must vary in direction ($`|𝐋|=`$const), which means the variation of $`\lambda _{LJ}`$ ($`𝐉`$ is invariable). Therefore, by Eq(7), a variable $`\dot{\omega }^{obs}`$ is expected.

The change of $`\lambda _{LJ}`$ also means the orbital plane tilts back and forth, as shown in Fig 1, in turn both $`\lambda _{LS}`$ and $`\lambda _{JS}`$ vary with time. Therefore, by Eq(6), the derivative of the rate of orbital precession can be given by,

$$\dot{\mathrm{\Omega }}_0=\mathrm{\Omega }_2\mathrm{\Omega }_{12}X_3X_4\mathrm{\Omega }_{12}X_1(\mathrm{\Omega }_{01}X_2+\mathrm{\Omega }_{12}X_3),$$
(8)
where $`\mathrm{\Omega }_{12}=\mathrm{\Omega }_1\mathrm{\Omega }_2`$, $`\mathrm{\Omega }_{01}=\mathrm{\Omega }_1\mathrm{\Omega }_0`$, $`X_1=\frac{S_1^{}}{S}\mathrm{sin}\lambda _{LS_1}`$, $`X_2=\mathrm{tan}\eta _{ss1}`$, $`X_3=\frac{S_{V1}S_{V2}}{S^2}\frac{\mathrm{sin}\eta _{s1s2}}{\alpha \mathrm{sin}\lambda _{JS}}`$, and $`X_4=\frac{\mathrm{cos}^2\lambda _{LS}}{\mathrm{sin}\lambda _{LS}}`$; with $`\alpha =\mathrm{sin}\lambda _{JS}+\frac{\mathrm{cos}^2\lambda _{LS}}{\mathrm{sin}\lambda _{LS}}`$, $`S_{V1}=S_1\mathrm{sin}\lambda _{JS1}`$ and $`S_{V2}=S_2\mathrm{sin}\lambda _{JS2}`$ represent components of $`𝐒_\mathrm{𝟏}`$ and $`𝐒_\mathrm{𝟐}`$ that are vertical to $`𝐉`$.

Note that $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ are unchanged, and $`\lambda _{LS_\alpha }`$ are unchanged, since they decay much slower than that of the orbitapo . The second derivative of $`\mathrm{\Omega }_0`$ are given by,

$$\ddot{\mathrm{\Omega }}_0=\mathrm{\Omega }_2\mathrm{\Omega }_{12}(\dot{X}_3X_4+X_3\dot{X}_4)\mathrm{\Omega }_{12}\dot{X}_1(\mathrm{\Omega }_{01}X_2+\mathrm{\Omega }_{12}X_3)$$

$$\mathrm{\Omega }_{12}X_1(\mathrm{\Omega }_{01}\dot{X}_2+\mathrm{\Omega }_{12}\dot{X}_3),$$
(9)
where $`\dot{X}_\alpha `$ is the first derivative of $`X_\alpha `$. $`\stackrel{\dot{}\dot{}\dot{}}{\mathrm{\Omega }}_0`$ can be easily obtained from Eq(9).

$`\dot{\mathrm{\Omega }}_0`$, the derivative of $`\mathrm{\Omega }_0`$, can be absorbed by $`\dot{P}_b`$. The variation in the precession velocity of the orbit results in a variation of orbital frequency ($`\nu _b=2\pi /P_b`$), $`\nu _b^{}\nu _b=\dot{\mathrm{\Omega }}_0\mathrm{\Delta }t`$, then we have $`\dot{\nu }_b=\dot{\mathrm{\Omega }}_0`$ and $`\ddot{\nu }_b=\ddot{\mathrm{\Omega }}_0`$, therefore,

$$\dot{P}_b=\frac{\dot{\mathrm{\Omega }}_0P_b^2}{2\pi },\ddot{P}_b=\frac{\ddot{\mathrm{\Omega }}_0P_b^2}{2\pi }\frac{\dot{\mathrm{\Omega }}_0P_b\dot{P}_b}{\pi }.$$
(10)
From Eq(8) and Eq(10), we can see that the contribution of the orbital precession to $`\dot{P}_b`$ is 1 PPN, which is much larger than the contribution of GR to $`\dot{P}_b`$. Since it depends on parameters, $`X_1`$, $`X_2`$, $`X_3`$, and $`\lambda _{LS}`$ of Eq(8), then it is also possible that the geodetic precession induced $`\dot{P}_b`$ is comparable to, or even smaller than the prediction of GR in special cases (special combination of $`\lambda _{LS}`$, $`\lambda _{LS1}`$, $`\eta _{SS1}`$ and $`S_1/S`$).

The effect of S-L coupling on secular evolution of the orbital inclination, $`i`$, can be given as,

$$\mathrm{cos}i=\mathrm{cos}\lambda _{LJ}\mathrm{cos}I\mathrm{sin}\lambda _{LJ}\mathrm{sin}I\mathrm{cos}\eta _0,$$
(11)
where $`I`$ is the angle between the total angular momentum, $`𝐉`$, and LOS, and $`\eta _0=\mathrm{\Omega }_0t+\eta _i`$ ($`\eta _i`$ is the initial phase) is the phase of precession of $`𝐋`$. Thus $`i`$ is also a function of time. By Eq(11) we have,

$$di/dt=\frac{\mathrm{\Omega }_0S}{L}\mathrm{sin}\eta _0+O(r^{7/2}).$$
(12)
Thus all parameters that are related with the orbital inclination will change with the geodetic precession. Hence the projected semi-major axis, $`x=a_1\mathrm{sin}i/c`$, will vary as the precession of the orbital plane. By Eq(12), the derivatives of the projected semi-major axis are

$$\dot{x}=x\mathrm{\Omega }_0\mathrm{sin}\lambda _{LJ}\mathrm{sin}\eta _0\mathrm{cot}i,$$
(13)
$$\ddot{x}=\dot{x}(\dot{\mathrm{\Omega }}_0/\mathrm{\Omega }_0+\mathrm{\Omega }_0\mathrm{cot}\eta _0).$$
(14)
As discussed between Eq(7) and Eq(8), the change of $`S`$ leads to the change of $`\mathrm{\Omega }_0`$, and in turn $`\dot{P}_b`$ and $`\ddot{P}_b`$. If there is only one spin, then $`S`$ is a constant, thus there will be only a static addition ($`\mathrm{\Omega }_0=const`$) to the apsidal motion, as shown in Eq(7). In this case, the orbital precession is also static, which will only cause a variation in $`\dot{x}`$, as shown in Eq(13), whereas, $`\dot{P}_b`$ and $`\ddot{P}_b`$ are not influenced at all ($`\dot{\mathrm{\Omega }}_0=0`$ and $`\ddot{\mathrm{\Omega }}_0=0`$). This is why the quadrupole-orbit model (also the one spin S-L coupling model) cannot explain $`\dot{P}_b`$ and $`\ddot{P}_b`$ in PSR J2051$``$0827 and PSR B1957$`+`$20.

## IV effects of geodetic precession induced $`\mathrm{\Omega }_0`$

The S-L coupling induced $`\mathrm{\Omega }_0`$ leads to the precession of the orbital plane, then the orbital inclination, $`i`$, varies, as in Eq(12), and in turn the projected semi-major axis, $`\dot{x}`$, varies, as in Eq(13). Comparing the predicted $`\dot{x}`$ with the observational $`\dot{x}^{obs}`$, we can constrain the angle $`\lambda _{LJ}S/L`$, and therefore the moment of the inertia of the pulsar.

### IV.1 PSR J2051$``$0827

The measured orbital period derivatives of PSR J2051$``$0827do  are list in Table I, the derivatives of the semi-major axis are $`\dot{x}^{obs}=23(3)\times 10^{14}`$ and $`\ddot{x}^{obs}=7(15)\times 10^{22}`$ s-1. The measured ones are much larger than the corresponding predictions of GR, $`\dot{P}_b^{GR}=3.8\times 10^{14}`$ and $`\dot{x}^{GR}=1.3\times 10^{19}`$, as well as the proper motion induced effects, $`\dot{x}^{pm}=4(2)\times 10^{17}`$ko .

With $`m_1=1.40M_{}`$, $`m_2=0.03M_{}`$ and $`P_b=0.09911`$ddo , we have the semi-major axis, $`a=(GM/\nu _b^2)^{1/3}=7.08\times 10^{10}`$ cm, and orbital angular momentum, $`L=\nu _b\mu a^2(1e^2)^{1/2}=2.17\times 10^{50}`$g cm2s-1, with $`\mu `$ the reduced mass. From Eq(2), we have $`\mathrm{\Omega }_2=3.25\times 10^9`$s-1, and from Eq(6), $`|\mathrm{\Omega }_0|\mathrm{\Omega }_23\times 10^9`$s-1. With $`x=x^{obs}=0.045052`$, $`|\dot{x}^{obs}|=2.3\times 10^{13}`$, Eq(13) becomes

$$0.045\times (3\times 10^9)\mathrm{sin}\lambda _{LJ}|\mathrm{sin}\eta _0\mathrm{cot}i|=2.3\times 10^{13}.$$
(15)
Thus $`\mathrm{sin}\lambda _{LJ}|\mathrm{sin}\eta _0\mathrm{cot}i|2\times 10^3`$. If $`|\mathrm{sin}\eta _0\mathrm{cot}i|1\times 10^1`$, then $`\lambda _{LJ}2\times 10^2`$, which means $`S/L2\times 10^2`$. With $`L`$ obtained above, we have $`S_1S4\times 10^{48}`$g cm2s-1. Having the measured pulsar period, $`P=4.50864`$ms, the moment of inertia of the pulsar is $`I_13\times 10^{45}`$g cm2.

### IV.2 PSR B1957$`+`$20

The negative orbital period derivative changes to positive during the observation of PSR B1957$`+`$20aft . The measured derivatives of the orbital period are shown in Table 2. The measured upper-limit to $`\dot{x}`$ is $`\dot{x}^{obs}<3\times 10^{14}`$. The measured ones are also much larger than the corresponding GR predictions, $`\dot{P}_b^{GR}=3.4\times 10^{15}`$ and $`\dot{x}^{GR}=6.1\times 10^{21}`$, as well as the proper motion induced effect, $`\dot{x}^{pm}=2.4\times 10^{16}`$ko .

Similarly, with $`m_1=1.40M_{}`$, $`m_2=0.025M_{}`$ and $`P_b=33001.91484`$saft , we have $`a=1.74\times 10^{11}`$ cm, $`L=2.82\times 10^{50}`$g cm2s-1, $`\mathrm{\Omega }_2=3.43\times 10^{10}`$s-1, and in turn $`|\mathrm{\Omega }_0|\mathrm{\Omega }_23\times 10^{10}`$s-1. From Eq(13), with $`x=x^{obs}=0.0892253`$ and $`|\dot{x}^{obs}|<3\times 10^{14}`$, we have

$$0.089\times (3\times 10^{10})\mathrm{sin}\lambda _{LJ}|\mathrm{sin}\eta _0\mathrm{cot}i|<3\times 10^{14},$$
(16)
which leads to , $`\mathrm{sin}\lambda _{LJ}|\mathrm{sin}\eta _0\mathrm{cot}i|<1\times 10^3`$. If $`|\mathrm{sin}\eta _0\mathrm{cot}i|1\times 10^1`$, then $`\lambda _{LJ}<1\times 10^2`$. Which means $`S/L<1\times 10^2`$, and by the obtained $`L`$, we have $`S_1S<3\times 10^{48}`$g cm2s-1. Finally, with the measured pulsar period, $`P=1.607`$ms, the moment of inertia of the pulsar is $`I_1<8\times 10^{44}`$g cm2 .

The simple estimation above indicates that PSR B1957$`+`$20 has similar moment of inertia as that of PSR J2051$``$0827.

In these two estimations, the assumptions $`|\mathrm{\Omega }_0|\mathrm{\Omega }_2`$ and $`|\mathrm{sin}\eta _0\mathrm{cot}i|1\times 10^1`$ are used. Consider the deviation of these two assumptions from the corresponding true values, the most likely moment of inertia from the observational constraints is about $`10^{44}`$g cm2$`10^{45}`$g cm2.

## V effects of geodetic precession induced $`\dot{\mathrm{\Omega }}_0`$, $`\ddot{\mathrm{\Omega }}_0`$, $`\stackrel{\dot{}\dot{}\dot{}}{\mathrm{\Omega }}_0`$

From Eq(8) and Eq(9), we can see that $`\dot{\mathrm{\Omega }}_0`$, $`\ddot{\mathrm{\Omega }}_0`$ and $`\stackrel{\dot{}\dot{}\dot{}}{\mathrm{\Omega }}_0`$ can vary in wider and wider range, which can be written approximately as (for the convenience of estimation),

$$|\dot{\mathrm{\Omega }}_0|c_1\mathrm{\Omega }_2^2,|\ddot{\mathrm{\Omega }}_0|c_2\mathrm{\Omega }_2^3,|\stackrel{\dot{}\dot{}\dot{}}{\mathrm{\Omega }}_0|c_3\mathrm{\Omega }_2^4,$$
(17)
where $`c_1`$ represents the contribution of $`X_\alpha `$ to $`\dot{\mathrm{\Omega }}_0`$, $`c_2`$ corresponds to the contribution of $`X_\alpha `$ and $`\dot{X}_\alpha `$ to $`\ddot{\mathrm{\Omega }}_0`$, and $`c_3`$ corresponds to the contribution of $`X_\alpha `$, $`\dot{X}_\alpha `$ and $`\ddot{X}_\alpha `$ to $`\stackrel{\dot{}\dot{}\dot{}}{\mathrm{\Omega }}_0`$. Obviously $`\stackrel{\dot{}\dot{}\dot{}}{\mathrm{\Omega }}_0`$ can vary in a wider range than $`\ddot{\mathrm{\Omega }}_0`$, and $`\ddot{\mathrm{\Omega }}_0`$ can vary in a wider range than $`\dot{\mathrm{\Omega }}_0`$.

For PSR J2051$``$0827, $`\mathrm{\Omega }_2^29.8\times 10^{18}`$s-2, $`\mathrm{\Omega }_2^33.2\times 10^{26}`$s-3 and also (equation not shown) $`\mathrm{\Omega }_2^41.0\times 10^{34}`$s-4. Assuming $`c_1=c_2=c_3=0.2`$, and by Eq(17) and Eq(10), the derivatives of $`P_b`$ can be obtained, $`|\dot{P}_b|2\times 10^{11}`$, $`|\ddot{P}_b|8\times 10^{20}`$s-1 and $`|\stackrel{\dot{}\dot{}\dot{}}{P}_b|2\times 10^{28}`$s-2. Which can be well consistent with the corresponding observations, as shown in Table I.

For PSR B1957$`+`$20, the measured derivatives of $`P_b^{obs}`$ is larger than that of PSR J2051$``$0827. By assuming $`c_1=1`$, $`c_2=1\times 10^3`$ and $`c_3=1\times 10^4`$, and through the same treatment, we have $`|\dot{P}_b|2\times 10^{11}`$, $`|\ddot{P}_b|7\times 10^{18}`$s-1 and $`|\stackrel{\dot{}\dot{}\dot{}}{P}_b|2\times 10^{26}`$s-2, which can be consistent with the corresponding observations, as shown in Table II.

It is likely that in PSR B1957$`+`$20 the magnitude of $`S_1`$ and $`S_2`$ are close, so that $`S`$ can be small, which leads to a large $`X_\alpha `$, and in turn large derivatives of $`P_b`$ as shown in Eq(6) and Eq(8)-Eq(10).

## VI discussion and conclusion

Geodetic precession induced pulsar spin in binary pulsar PSR B1913+16 has been studied through structure parameters of pulsar profilekr ; wt . Actually the variation of pulsar spin axis (former) and the measured secular variabilities discussed in this paper are induced by the same physics underlying in binary pulsar systems, S-L coupling. The former corresponds to the reaction of the pulsar spin to the torque in the S-L coupling; while the latter corresponds to reaction of the orbit to the same torque.

By adding one more spin into the S-L coupling of Apostolatos et al and Kidderapo ; kid , we establish the relationship among the secular variability, geodetic precession and moment of inertia of binary pulsars. The new model provides: (a) a unified model that explain both $`\dot{x}`$ and $`\dot{P}_b`$, which has been interpreted separately by different models, (b) a new method which can extract the moment of inertia through pulsar timing measurement, (c) a new test of the geodetic precession which has strong effects to pulsar timing.


text:
# The Gap-Tooth Method in Particle Simulations

## Abstract

We explore the gap-tooth method for multiscale modeling of systems represented by microscopic physics-based simulators, when coarse-grained evolution equations are not available in closed form. A biased random walk particle simulation, motivated by the viscous Burgers equation, serves as an example. We construct macro-to-micro (lifting) and micro-to-macro (restriction) operators, and drive the coarse time-evolution by particle simulations in appropriately coupled microdomains (“teeth”) separated by large spatial gaps. A macroscopically interpolative mechanism for communication between the teeth at the particle level is introduced. The results demonstrate the feasibility of a “closure-on-demand” approach to solving hydrodynamics problems.

modeling; multiscale; closure-on-demand; gap-tooth; lifting; restriction  preprint: Physics Letters A

Traditional approaches to solving physical problems that manifest separation of scales involve first (a) deriving a set of reduced equations to describe the system, and subsequently (b) solving the equations and analyzing their solutions. Recently an “equation-free” approach has been proposed Kevrekidis02  that sidesteps the necessity of first deriving explicit reduced equations. The approach relies instead on microscopic simulations, enabling them through a computational superstructure to perform numerical tasks as if the reduced equations were available in closed form. Both macroscopically coarse-grained equations and atomistic/stochastic simulations can be regarded as “black boxes” from the point of view of appropriately formulated numerical algorithms. They constitute alternative realizations of the same macroscopic input-output mapping. For example, a crystal’s elastic response can either be evaluated using elastic constants, or estimated by a high-accuracy electronic structure program based on density functional theory, which, for a given strain, computes the stress on-the-fly. The advantage of a simulator-based approach is that it can be used generally, beyond the region of validity of any given closure - e.g. providing the correct nonlinear elastic responses in the above example. Equation-free methods hold the promise of combining direct physics-based simulation with the strength and scope of traditional numerical analysis on coarse variables (bifurcation, parametric study, optimization) for certain problems - problems for which coarse equations conceptually exist, but are not available in closed form. An example is the so-called interatomic potential finite-element method (IPFEM)Li02 , a subset of the more general quasi-continuum methodPhillips99 , used to identify elastic instabilities leading to defect nucleation in nanoindentation, for which no accurate closed-form constitutive relation is currently available due to the complex triaxial stress state at the critical site of instability.

Microscopic simulations cannot be used directly to attack problems with large spatial and temporal scales (“macrodomains” in space and time); the amount of computation is prohibitive. If, however, the actual behavior can be meaningfully coarse-grained to a representation that is smooth over the macrodomain, the microscopic systems need only be directly simulated in small patches of the macrodomain. This is done by interpolating hydrodynamic variables between the patches in space - the gap-tooth method (see Kevrekidis00 ) - and extrapolating from one or more patches in time - projective integrationGear01a ; Gear02a . In this paper, we use this “closure-on-demand” approach to solve for the coarse-grained behavior of a particular microscopic system. The illustrative example is the biased random walk of an ensemble of particles, motivated by the viscous Burgers equation,

$$u_t+uu_x=\nu u_{xx},$$
(1)
a 1D version of the hydrodynamics equations used under various conditions to model boundary layer behavior, shock formation, turbulence, and transport. Here, $`\nu >0`$ is the viscosity; periodic boundary condition is used for simplicity, and only non-negative solutions $`u(x,t)>0`$ are considered. A particular microscopic dynamics is constructed, motivated by Eq.(1), interpreting $`u`$ as the density field of the random walkers; $`u𝑑x=1`$ corresponds to $`Z`$ walkers, where $`Z`$ is a large normalization constant. In the micro-simulation, random walkers move on $`[1,1)`$ at discrete timesteps $`t_n=nh`$. At each step, an approximation to the local density, $`\rho _i`$, is computed from the positions of neighbors. Then every walker is moved by $`\mathrm{\Delta }x_iN(h\rho _i/2,2\nu h)`$, a biased Gaussian distribution. $`x_i`$’s are then wrapped around to $`[1,1)`$, and the process repeats. Since $`\rho _i`$ is a local estimate of $`u`$, this process achieves a coarse-grained flux analogous to $`ju^2/2\nu u_x`$ in Eq.(1) by assigning each walker a drift velocity of $`\rho _i/2`$.

The gap-tooth scheme, first discussed in Kevrekidis00 , covers space with teeth and intervening gaps as shown in Fig. 1 for one dimension. The microscopic evolution is simulated in the interior of each tooth. Clearly appropriate boundary conditions have to be provided at the edges of each tooth. Tooth boundaries coincident with external boundaries have the boundary conditions specified externally, while internal boundary conditions must be generated by the gap-tooth scheme itself. Because this example uses periodic boundary conditions, there are no external boundaries: the teeth can be viewed as equally spaced on a circle.

The microscopic simulation operates on the position of each particle. We are interested in a meaningful coarse description, a finite-dimensional approximation to the density of particles, $`u(x)`$, possibly averaged over several realizations of the computational experiment Kevrekidis02 . The lifting operator that maps a given $`u(x)`$ to consistent particle positions is straightforward in this case. From the density function over a tooth we can compute its integral, so we know the number of particles that should be present in that tooth. The indefinite integral of the density function over the tooth provides the cumulative distribution function for that tooth which permits the particles to be placed as a discrete representation of that functioncwgdist . If the density approximation is constant in each tooth (as has been found to be adequate in the examples here) this simply means that the particles are uniformly spaced in each tooth according to the density in that tooth.

In our particular stochastic simulation, the evolution rules require a local density estimate. This should be done by choosing a particle density influence function $`\sigma (d)`$ that specifies the contribution of each particle to the local density at a distance $`d`$ from that particle. If this function is constant for $`d<w`$ and zero elsewhere, the local density function at $`x`$ is found by just counting the number of particles within distance $`d`$ of $`x`$. Ref. Chertock01 ; Chertock02 , seeking some level of differentiability, use a Gaussian spreading function for each particle. Since we do not require differentiability, we will count particles within distance $`d`$. By making $`d`$ twice the tooth size, we can then simply count the particles in each tooth. We have also used higher-order approximations, but it is not clear that they yield sufficient improvements in accuracy to justify the additional computational effort. The technique used for higher-order approximations is based on the fact that the sample cumulative distribution function in each tooth is known from the particle positions. We can then fit a polynomial to it within each tooth. The density function over a tooth is simply given by the derivative of this polynomial.

The mapping of a phase point or points (particle positions and velocities history) to coarse fields is called a restriction operator. In addition to the density field ($`0`$th-moment), smooth velocity ($`1`$st-moment) and temperature ($`2`$nd-moment) fields can be extracted from molecular dynamics based on maximum likelihood inference Li98 . If the interior of a tooth were to be simulated by solving a PDE, we would need to prescribe appropriate boundary conditions at each tooth at each timestep. The same is still true when the tooth is realized using particle-based microscopic simulations. Creating an appropriate match between the coarse fields at the boundaries and the particles in the teeth is an area of intense researchLi99 ; E01 . Sometimes one knows so little about the nature of the coarse equation that even the correct order for imposing well-posed boundary conditions at the teeth is unknown. This issue is addressed in Li03 .

Here we use an alternative approach suggested in Gear02b . In a 1D particle based random walk simulation we can distinguish two “fluxes” - left-going and right-going. The particle simulation in the interior of each tooth generates outgoing fluxes, that is, the left\[right\]-going fluxes at the left\[right\] boundaries, directly. Boundary conditions are needed to provide matching incoming (right\[left\]-going) fluxes at the same boundaries. In $`d`$-dimensions, there will be $`2^d`$ boundaries to deal with and the corresponding incoming fluxes to provide.

Consider the estimation of the right-going incoming flux, $`\mathrm{I}_{r,1}`$, as shown in Fig. 2.

Assuming macroscopic flux smoothness suggests that we can interpolate its values from neighboring outgoing fluxes, in this case $`\mathrm{O}_{r,0}`$ and $`\mathrm{O}_{r,1}`$. If we use linear interpolation, we can write

$$\mathrm{I}_{r,i}=\alpha \mathrm{O}_{r,i1}+(1\alpha )\mathrm{O}_{r,i}.$$
(2)
The interpolation coefficients depend (in this case through $`\alpha `$) only on the gap-tooth geometry.

However, the “fluxes” under discussion here are not real-valued quantities, but discrete events as particles cross a boundary, so Eq.(2) needs a different interpretation. Consider instead the role played by each outgoing flux in the interpolation for ingoing fluxes. An interpretation of Eq.(2) for $`i=1`$ and $`i=2`$ would be the portion $`(1\alpha )`$ of $`\mathrm{O}_{r,1}`$ contributes to the flux $`\mathrm{I}_{r,1}`$ while $`\alpha `$ of it contributes to $`\mathrm{I}_{r,2}`$. A similar procedure applies to the left-going fluxes. Thus, rather than thinking in terms of flux interpolation we can think in terms of flux redistribution. Interpreting the linear interpolation stochastically, (on a regularly spaced gap-tooth scheme) we direct $`\alpha `$ of the outgoing particles as input to the neighboring tooth, and redirect $`(1\alpha )`$ of them back to the other boundary of the same tooth as shown in Fig. 3.

Flux redistribution has to recognize the position of a particle after it leaves a tooth. If it had moved to a distance $`\delta `$ beyond the boundary of the tooth, it must be inserted a distance $`\delta `$ inside the receiving tooth. If $`\delta `$ were larger that the tooth width it would have left a tooth boundary again and a further redistribution would be required following the same rule. (In multiple dimensions, the boundaries in each dimension are treated independently so that a particle will be distributed for each boundary that it crosses until it is inside a tooth.)

The above method implements effective linear interpolation. As discussed in Gear02b , linear interpolation is not adequate for second-order problems: at least quadratic interpolation must be used. A possible quadratic interpolation formula is

$$\mathrm{I}_{r,i}=\frac{\alpha (1+\alpha )}{2}\mathrm{O}_{r,i1}+(1\alpha ^2)\mathrm{O}_{r,i}\frac{\alpha (1\alpha )}{2}\mathrm{O}_{r,i+1}.$$

As before we consider the impact of each outgoing flux on incoming fluxes. The fractions of output $`\mathrm{O}_{r,1}`$ should be sent to the inputs as follows: $`(1\alpha ^2)`$ to $`\mathrm{I}_{r,1}`$; $`\alpha (1+\alpha )/2`$ to $`\mathrm{I}_{r,2}`$: and $`\alpha (1\alpha )/2`$ to $`\mathrm{I}_{r,0}`$.

Note that the last value is negative. Any linear higher-order interpolation formula contains negative coefficients. Our solution is to direct anti-particles to the appropriate teeth. There they must annihilate with regular particles - we simply annihilate with the nearest regular particle. With this approach, the $`\mathrm{O}_{r,1}`$ is redistributed as follows: $`(1\alpha (1+\alpha )/2)`$ to $`\mathrm{I}_{r,1}`$; $`\alpha ^2`$ to $`\mathrm{I}_{r,2}`$; and $`\alpha (1\alpha )/2`$ are cloned to get two regular particles sent to $`\mathrm{I}_{r,1}`$ and $`\mathrm{O}_{r,2}`$ and one anti-particle sent to $`\mathrm{I}_{r,0}`$ It is interesting to observe that this scheme conserves the total number of particles in the computation.

The microscopic evolution rules were simulated at conditions corresponding to $`\nu =0.05`$ and timestep $`h=0.002`$ in Eq.(1) for 1,000 time steps using the gap-tooth scheme with 40 equally spaced teeth in the interval $`[1,1)`$. The tooth-to-spacing ratio was $`\alpha =0.2`$. The coarse initial condition was $`1\mathrm{sin}(\pi x)/2`$. Fig. 4 shows the results at $`t=2`$ when 100,000 particles are used (upper figure) and 500,000 particles (lower figure). The dots indicate the particle density approximations within each tooth (constant in this case); they are connected by a spline interpolant. The smooth curve (to guide the eye) is the analytical solution Cole51 ; Hopf50  of Eq.(1).

The problem was also run with $`\alpha =1`$, i.e., no gaps, to compare with a conventional, full space, direct microscopic simulation. Fig. 5 shows the result using 500,000 particles and a constant density approximation in each tooth. This is the same average particle density per spatial unit as with $`\alpha =0.2`$ and 100,000 points. The results are comparable. The $`\alpha =1`$ case with 40 teeth and 500,000 particles has also been run with the polynomial density interpolants described above of degrees 1 and 2. No significant difference in the solutions were observed, and neither were differences observed when the higher order polynomial fits were used with $`\alpha =0.2`$.

We have demonstrated that the gap-tooth scheme can be successful in solving some problems using microscopic models based on the stochastic simulation of particle motion; we also introduced a novel approach for dealing with the inter-tooth boundary conditions.

In earlier work we have proposed combining this with projective integration Gear01a  and some preliminary experiments have been performed in this direction. However, projective integration requires smoothness of the time derivative estimates. The stochastic nature of the microscopic model leads to significant noise. In this simulation we used a least-squares linear estimate from a large number of time steps to get a reasonably accurate time derivative estimate. However, by then, the total time step at the microscopic level was large relative to the size of a projective step in time (which is limited by the smoothness of the solution). If the stochastic noise is reduced by using a much larger number of particles, or a large number of “copies” of the simulation, the time derivative estimates would allow the application of projective integration. In some sense there is a tradeoff between saving computation in the spatial domain with fewer teeth and fewer particles, and saving computation in the time domain by getting more accurate estimates of the time derivatives. We believe that for problems where there is a significant gap between the timescales of the microscopic dynamics and those of the expected macro-scale behavior, projective intregration would be a useful addition. We will report on such “patch dynamics” experiments in a future paper.

This work was partially supported by the Air Force Office of Scientific Research (Dynamics and Control) and an NSF ITR grant (C.W.G.,I.G.K.). J.L. acknowledges support by Honda R&D Co., Ltd. (RF# 744252) and the OSU Transportation Research Endowment Program.


text:
# 1 Introduction

## 1 Introduction

Random matrix theory has matured into a field with applications in many branches of physics and mathematics. A large number of physicists and mathematicians have been fascinated by this subject that has turned out to be surprisingly rich and far reaching. Paraphrasing Dyson, random matrix theory is a new kind of statistical mechanics where the realization of the system is not relevant. Instead of having an ensemble of states we have an ensemble of Hamiltonians. Ergodicity is now the equivalence of spectral averaging and the averaging over this ensemble.

Random matrix theory has been particularly successful in three areas: first, in describing level correlations on the scale of the average level spacing; second, in providing the generating function for combinatorial factors of planar diagrams; and, third, as an exactly solvable model with intricate connections to the theory of integrable systems. One of the reasons for the success of random matrix theory is universality: eigenvalue correlations on the scale of the average level spacing do not depend on the probability distribution. This property is at the very foundation of random matrix theory. It suggests that random matrix theory correlations of eigenvalues should be the rule rather than the exception. However, the most important reason for studying random matrix theory is that its predictions do occur in nature in systems as varied as nuclear energy levels, zeros of the Riemann $`\zeta `$ function and sound waves in quartz crystals. Another of the roles of random matrix theory is that the large $`N`$ limit of its partition function is a generating function for planar diagrams which have played an important role in quantum field theory. For example, they are the leading contributions to QCD with a large number of colors, and they are dual to triangulations of a random surface and thus describe two-dimensional quantum gravity. In addition to this, random matrix theory has attracted a great deal of interest because of the mathematical challenges it poses. The problems are highly nontrivial, but, with sufficient effort, many of the questions that arise in this field can be answered in full.

It is hard not to be fascinated by random matrix theory. Everyone who works in this field has experienced the amazement of obtaining truly universal behavior by diagonalizing a large random matrix. Nowadays, as we can see from the contributions to this special edition, the excitement about the subject is still as much alive as when it was first created. In the introduction of this special edition we give a short summary of the history of random matrix theory. The historical perspective is certainly colored by our personal experience, and, for a somewhat different perspective, we refer to the introduction in the review by the Heidelberg group . We conclude with a short overview of the papers that will appear in this special issue.

## 2 History of Random Matrix Theory

Random matrix theories have fascinated both mathematicians and physicists since they were first introduced in mathematical statistics by Wishart in 1928 . After a slow start, the subject gained prominence when Wigner  introduced the concept of statistical distribution of nuclear energy levels in 1950. However, it took until 1955 before Wigner  introduced ensembles of random matrices. In this paper he also introduced the large $`N`$ expansion and realized that the leading order contribution to the expectation values of moments of the random Hamiltonian is given by planar diagrams. In 1956, Wigner  derived the Wigner surmise from the level spacing distribution of an ensemble of $`2\times 2`$ matrices after level repulsion had been predicted by Landau and Smorodinksy  and observed by Gurevich and Pevsner . The idea of invariant random matrix ensembles was introduced in physics by Porter and Rosenzweig  after it had appeared earlier in the mathematical literature. A mathematically rigorous analysis of spacing distributions was first given by Gaudin  and Mehta . To analyze the eigenvalue density Mehta  invented the orthogonal polynomial method.

The mathematical foundations of random matrix theory were established in a series of beautiful papers by Dyson . He introduced the classification of random matrix ensembles according to their invariance properties under time reversal . As we all know, only three different possibilities exist: a system is not time reversal invariant, or a system is time reversal invariant with the square of the time reversal invariance operator either equal to 1 or $`1`$. The matrix elements of the corresponding random matrix ensembles are complex, real and self-dual quaternion, respectively, which from a mathematical viewpoint exhaust the distinct real commutative normed division algebras, or in effect number systems. The corresponding invariant Gaussian ensembles of Hermitian random matrices are known as the Gaussian unitary ensemble (GUE), the Gaussian orthogonal ensemble (GOE) and the Gaussian symplectic ensemble (GSE), in this order.

Dyson  also formulated the underlying philosophy of random matrix theory. In his words, “What is here required is a new kind of statistical mechanics, in which we renounce exact knowledge not of the state of the system but of the system itself. We picture a complex nucleus as a “black box” in which a large number of particles are interacting according to unknown laws. The problem then is to define in a mathematically precise way an ensemble of systems in which all possible laws of interaction are equally possible”. This was made more precise by Balian , who obtained the Gaussian random matrix ensembles from minimizing the information entropy.

A second important result of Dyson’s papers  was the relation between random matrix theory and the theory of exactly integrable systems: the partition function of a random matrix ensemble is equivalent to the partition function of a log-potential Coulomb gas in one dimension at three special temperatures, each with solvability properties not shared for general temperatures. Moreover the evolution of the eigenvalues of parameter dependent extensions of the Gaussian ensembles was shown to be controlled by a Fokker-Planck operator which also specifies the Brownian evolution of the Coulomb gas. These results were further explored by Sutherland  when he realized that the Calogero-Sutherland quantum many body system , for which the Hamiltonian can be constructed from $`N`$ independent commuting operators and so is integrable, is mathematically equivalent to Dyson’s Brownian motion model. The relation between random matrix theory and integrable systems is discussed extensively in the monograph by Forrester . A review of one dimensional integrable systems that touches on many ideas that also appear in random matrix theory is given in the book by Korepin, Bogoliubov and Izergin . A third idea that appeared in Dyson’s paper  is the application of Shannon’s information entropy to random matrix spectra.

The early developments in random matrix theory are well summarized in the first edition of the monograph by Mehta . This has been a very influential book containing many mathematical details which have been proved to be extremely useful over the years. A second significant book is by Porter . It contains reprints of the important papers on random matrix theory that were written before 1965.

About the same time as the early development of random matrix theory in nuclear physics, the field of disordered systems was born by the work by Anderson  on the localization of wave functions in one dimensional disordered systems. He considered a one dimensional lattice with a random potential at each lattice point. He found that the eigenfunctions of this system are exponentially localized. His work had strong impact on both experimental and theoretical solid state physics. Another early application of random matrix theory is the theory of small metallic particles by Gorkov and Eliasberg , which nowadays would be part of mesoscopic physics.

Random matrix theory, which was first formulated in mathematical statistics, continued to develop in mathematics independently of the developments in physics. Important results with regards to the integration measure of invariant random matrix ensembles were obtained by L.K. Hua . His results of more than a decade of work are summarized in his book that appeared in 1959 but which remained largely unknown. Only a small number of mathematicians worked on integrals that appear in random matrix theory. One very important result was obtained by Harish-Chandra , who evaluated a unitary matrix integral that is now known as the Harish-Chandra-Itzykson-Zuber integral . Zinn-Justin and Zuber  review this topic in the present special issue. Also the work of Selberg  is well-known, not in the least because Madan Lal Mehta devoted a chapter of the second edition of his book  to this subject. Another noteworthy contribution is the introduction of zonal polynomials by James . The 1982 book of Muirhead  ties together matrix integrals and zonal polynomials as they are relevant in mathematical statistics. Girko has written a number of mathematical books (see e.g. ) relating to analytic properties of the eigenvalue distribution of large random matrices. Voiculescu  used random matrices as a primary example of the concept of free non-commutative random variables in operator algebras. However, the mathematical literature remained largely unnoticed by physicists until recently.

What is more surprising is that the theory of disordered systems and the application of random matrix theory in nuclear physics proceeded more or less independently until the seminal work by Efetov on the supersymmetric method  and its application  to the theory of small metallic particles and to localization theory  . This is even more remarkable since both the papers by Anderson and Dyson were written at Princeton.

The main developments in random matrix theory in the decade after the appearance of the first edition of Mehta’s book were applications to nuclear physics. In particular, the statistical theory of $`S`$-matrix fluctuations received a great deal of attention. The first work in this direction dates back to Wigner , who simultaneously studied the distribution of the widths and the spacings of nuclear resonances, and to Porter and Thomas , who introduced the Porter-Thomas distribution for nuclear decay widths. Correlations of cross-sections at two different energies were considered in  and are now known as Ericson fluctuations. The formulation of the theory of $`S`$-matrix fluctuations was completed in the work of Agassi et al.  . In this paper the authors introduced resummation techniques which later received much more attention in the field of impurity scattering, introduced earlier in the book by Abrikosov, Gorkov and Dzyaloshinski . The problem of the distribution of poles of $`S`$-matrices was also the motivation of Ginibre  for introducing what is now known as the Ginibre ensemble with eigenvalues uniformly distributed inside a disk in the complex plane. His paper initiated the subfield of nonhermitian random matrix theory which is reviewed by Fyodorov and Sommers  in this issue.

In 1973 Montgomery  made a conjecture for the asymptotic limit of the two-point correlation function of the zeros of the Riemann $`\zeta `$ function on the critical line. Together with Dyson he realized that his conjectured result is the two-point function of the GUE. The connection was extended to higher correlation functions of the Riemann zeros by Hejhal  and Rudnick and Sarnak , although the full correspondence of the correlation functions with random matrix theory has still not been proved. A heuristic derivation of these results using the Hardy-Littlewood conjecture for the correlation between primes, was given by Bogomolny and Keating in 1995 . Mathematically rigorous results relating the two-point functions for the zeros of families of finite field zeta functions and eigenvalues of random matrices from the classical groups are the topic of the monograph by Katz and Sarnak .

The conjectured correspondence of the statistics of these zeros of the Riemann $`\zeta `$ function with the $`n`$-point correlation function of random matrix eigenvalues has recently meant that random matrix theory has become very useful for conjecturing quantities in number theory that were previously unattainable by any method. These include mean values of the Riemann zeta function and other $`L`$-functions , the order of vanishing at special values of $`L`$-functions , as well as discrete moments of the derivative of the Riemann zeta function  and the horizontal distribution of the zeros of the derivative . For more details there is a review in this issue by Keating and Snaith .

In the period 1975-1985, random matrix theory developed rapidly and became unified with the theory of disordered systems. The first step in this direction was made by Edwards and Anderson  who, in their influential paper on spin glasses, introduced the replica trick. This provided a natural framework for a field theoretical formulation of the Anderson model which was introduced a few years later by Wegner in 1979 . In this formulation, symmetries and the spontaneous breaking of symmetries led to a new paradigm in the theory of Anderson localization . It was soon realized that the replica formulation only works well for perturbative calculations. This problem was solved by the introduction of the supersymmetric method . In this method the determinants in the generating function of the resolvent are quenched by taking a ratio of two determinants instead of the $`n0`$ limit of the $`n`$’th power of the determinant. Relying on earlier work by Wegner  using the replica trick, Efetov showed that the partition function of a disordered system is given by a supersymmetric nonlinear $`\sigma `$-model. He identified a domain of energy differences where the kinetic term of the non-linear $`\sigma `$ model can be neglected. In this domain the two-point correlation functions coincide with the results derived by Dyson. The energy scale below which the partition function is dominated by zero momentum modes is known as the Thouless energy .

The supersymmetric method has been very fruitful. Efetov  obtained new results for one-dimensional disordered wires. Exact results were obtained for the theory of $`S`$-matrix fluctuations . Relations between the orthogonal and symplectic symmetry classes were derived  from the supersymmetric partition function. In the subsequent years many more new results were derived by means of the supersymmetric method. Among others we mention results for parametric correlations  where eigenvalue correlations for different values of an external parameter are considered. An elaborate discussion of applications of the supersymmetric method to disordered systems is given in the book by Efetov .

Exact results for $`S`$-matrix fluctuations were obtained in a completely independent way by a Mexican group . The exact distribution function of $`S`$-matrices was found starting from the three assumptions of analyticity, ergodicity and maximizing the information entropy . Another effort in nuclear physics was the introduction of random matrix ensembles that reflected the few-body nature of the interaction. In particular, French and co-workers have pursued this direction of research (see  for a review). In this issue Benet and Weidenmueller  review recent progress in this field.

A major development was the experimental discovery of universal conductance fluctuations by Webb and Washburn in 1986  after having been predicted theoretically by Altshuler  and Stone and Lee . This discovery started the new field of chaotic quantum dots. The transport properties of these quantum dots could be described by the supersymmetric nonlinear $`\sigma `$-model that had been used for the theory of $`S`$-matrix fluctuations in compound nuclei. In fact, a compound nucleus is a chaotic quantum dot (see  and  in this issue for reviews).

A few years before the discovery of universal conductance fluctuations, random matrix theory was applied to quantum field theory. Through the work of ’t Hooft  we know that in the limit of a large number of colors, the QCD partition function is dominated by planar diagrams. This is also the case for the large $`N`$ limit of random matrix theory. In  this was exploited to calculate the combinatorial factors that enter in the large $`N_c`$ limit of QCD by means of random matrix theory. A second innovative idea that appeared in this paper is the formulation of the calculation of the resolvent in random matrix theories as a Riemann-Hilbert problem. This approach has received more attention in the recent mathematical literature .

Random matrix theory has had impact on several areas of quantum field theory: lattice QCD, two dimensional gravity, the Euclidean Dirac spectrum and the Seiberg-Witten  solution of two dimensional supersymmetric gauge theories. An important result is the Eguchi-Kawai  reduction. These authors show that in the limit of a large number of colors, certain gluonic correlation functions of pure Yang-Mills theory can be reduced to an integral over 4 unitary matrices. In two spatial dimensions this reduction results in an integral over a single unitary matrix which can be evaluated in the large-$`N`$ limit.

A unitary matrix integral also occurs in the low-energy limit of QCD. Because of the spontaneous breaking of chiral symmetry, its low-energy degrees of freedom are the Goldstone modes which are parameterized by a unitary matrix valued field . Below the Thouless energy for this system the kinetic term of the effective Lagrangian can be neglected and the low-energy limit of the QCD partition function is given by the unitary matrix integral . In this domain the eigenvalues of the Dirac operator are correlated according to a random matrix theory with the additional involutive (chiral) symmetry of the QCD Dirac operator . The same symmetry is also found in two-sublattice disordered systems where hopping only occurs in between the sublattices . The eigenvalue spectrum around zero of these chiral ensembles was first derived in . An important difference between two-sublattice systems and QCD is the topology of the random matrix (i.e. the number of exact zeros) and the fermion determinant. In two-sublattice systems one is only interested in quenched results at zero topology whereas in QCD the fermion determinant and its zero modes are essential. Also in the case of the chiral ensembles we have three different symmetry classes depending on the reality content of the matrix elements. Most of the work on chiral random matrix theory and its applications to the Dirac spectrum in QCD was done in the second half of the nineties (see  for a review).

In the theory of disordered superconductors four more random matrix ensembles can be introduced , thus distinguishing a total of 10 random matrix ensembles. It was noticed by Dyson  that each of the three Wigner-Dyson random matrix theories corresponds to a symmetric space. Zirnbauer  showed that this observation can be generalized to all 10 symmetry classes of random matrix ensembles with a one-to-one correspondence to each of the large families in the Cartan classification of symmetric spaces.

There have been other attempts to derive QCD from a matrix model. Perhaps best known is the induced QCD partition function of Kazakov and Migdal  where the lattice gauge field is coupled to an adjoint scalar field. The gauge field can be integrated out by means of the Harish-Chandra-Itzykson-Zuber integral resulting in a partition function for the eigenvalues of the adjoint scalar field. This partition function can be evaluated by saddle point methods in the large $`N`$ limit. More recently, it has been shown that that the so-called prepotential of $`N=2`$ supersymmetric theory can be derived from the large $`N`$ limit of a random matrix theory .

The partition function of 2d gravity is a sum over random surfaces which can be described by means of a triangulation . The sum over triangulated surfaces can be written in terms of a random matrix theory partition function. It has been conjectured  that the double scaling limit of this theory describes the continuum limit of the 2d gravity partition function. This field brought two new ideas into random matrix theory: universality , i.e. that observables are independent of the probability potential, and the connection with integrable systems. In the context of quantum gravity it is natural to consider an arbitrary polynomial probability potential. Integrable hierarchies were obtained from differential equations in the coefficients of the probability potential . A good review of this topic was given by Di Francesco, Ginsparg and Zinn-Justin .

Earlier integrable hierarchies entered in a completely different way. In 1980 it was found by the Kyoto school  that the probability of a gap free interval in the infinite GUE is a $`\tau `$-function for a completely integrable system specifying the isomonodromy deformation of a coupled system of linear differential equations. This had the consequence that the spacing distribution could be expressed in terms of a Painlevé V transcendent. Later it was found that the distribution of the largest eigenvalue in the GUE is given by the solution of a Painlevé II equation . This development found application in the solution of a long standing mathematical problem: specifying the limiting distribution of the longest increasing subsequence length of a random permutation . In fact the sought distribution is the same as that for the largest eigenvalue in the GUE. The increasing subsequence problem can equivalently be formulated as the polynuclear growth model in $`1+1`$ dimension , and similar relationships with random matrix fluctuations are also known for certain tiling problems .

The question of why random matrix theory works has been addressed from many different points of view. It was realized early on that the detailed properties of eigenvalue correlations do not depend on the specifics of the probability distribution. One important reason for random matrix theory to work is already mentioned in a work of Dyson , asserting that if a system is sufficiently complex, the state of the system is no longer important. However, it took until the early eighties before it was realized that the key reason is that the corresponding classical system is chaotic. Although there have been a few earlier studies relating random matrix theory correlations to classical chaos , it was formulated explicitly in a ground breaking paper by Bohigas, Gianonni and Schmit  who, based on a numerical study of the Sinai billiard , conjectured that level correlations on the scale of the average level spacing are given by random matrix theory if the corresponding classical system is chaotic. This conjecture has been confirmed for numerous systems. The reverse was also shown numerically to be true: if the system is not completely chaotic, the spectral correlations are not given by the Wigner-Dyson ensembles . Although a complete proof of this conjecture is still lacking, a considerable amount of analytical understanding has been obtained on the basis of a semiclassical analysis . These inter-relations mean that random matrix theory plays an essential role in the study of quantum chaos, a fact which is given prominence in the books by Haake  and Stoeckmann .

In this short historical overview we have seen that random matrix theory has been applied to wide ranging fields. Its scope has by far not yet been exhausted as illustrated by recent publications that are as varied as applications to financial correlations  and wireless communication .

## 3 The Special Issue

The present special edition represents a broad cross section of the current activity in random matrix theory. We have subdivided the forty-four contributions into eight different groups: Applications to Number Theory, Applications to Statistical Mechanics, Integrable Systems and Random Matrix Theory, Integration Formulas, Mesoscopic Physics and Disordered Systems, Non-Hermitian Random Matrix Theories, Quantum Chaos, and Special Random Matrix Ensembles.

The contributions to the section “Applications to Number Theory” address the connection between random matrix theory and the Riemann zeta function and other $`L`$-functions. In  and  (this second in conjunction with companion papers , published elsewhere) this connection is examined through the comparison of statistics of the $`L`$-function zeros with random matrix eigenvalues, while ,  and  make use of a comparison between the characteristic polynomial of a random matrix and the $`L`$-functions themselves. The present status of the application of random matrix theory to $`L`$-functions is reviewed by Keating and Snaith .

In section two we discuss applications of random matrix theory to statistical mechanics. Three of the papers  can be regarded as outgrowths of the solution of the largest increasing subsequence problem for a random permutation. Also there are papers on financial correlations , random projections relating to phase retrieval , and on energy landscape statistics .

Already since the work of Dyson random matrix theories have been interpreted as one-dimensional integrable systems. In particular, the Calogero-Sutherland model could be solved by means of techniques with their origin in random matrix theory. In this section we have one paper on the BC-type Calogero model . More recently it was found that some random matrix theory partition functions are $`\tau `$ functions of integrable hierarchies. This is explored in . Other relations with integrable structures are discussed in .

Different types of integration formulae relevant to random matrix theory are discussed in section four. A very important class of integrals are the Harish-Chanda-Itzykson-Zuber integrals which are reviewed by Zinn-Justin and Zuber . Other unitary matrix integrals are discussed in . New results for ratios of determinants are obtained in .

The fifth section is on mesoscopic physics and disordered systems. Although this field is strongly rooted in experimental condensed matter physics we did not receive any such contributions. However, this section contains contains two theoretical reviews, one on transport through a chaotic quantum dot by Polianski and Brouwer  and one on wave function statistics by Merlin, Evers and Mildenberger . Three different disordered systems are discussed, a two-sublattice system , a one-dimensional disordered system  and a microwave cavity .

The sixth section is on Non-Hermitian Random Matrix Theories. We can distinguish two different cases: weak non-hermiticity, which is reviewed by Fyodorov and Sommers , and strong non-hermiticity. In the first case the imaginary part of the eigenvalues is of the order of the spacing of the real parts. In the second case the real and imaginary parts of the eigenvalues are of the same order of magnitude which allows for much richer eigenvalue spectra than for Hermitian matrices (see , ). In addition to the review by Fyodorov and Sommers, eigenvalue correlations of different types of non-Hermitian matrices are studied in ,  and . The level spacings of a class of non-Hermitian matrices with real eigenvalues is studied in , while  treats the statistical properties of the zeros of a random analytic function which has close analogies with complex eigenvalues.

During the past two decades, quantum chaos has been one of the main anchors of random matrix theory. This special issue also has a section on this topic. One important issue has been to show (mostly numerically) that properties of simple chaotic systems can be described by random matrix theory. In this section we have a discussion of the eigenvalue spectra of the cat map , of scattering matrices on quantum graphs  and of action correlations and random matrix theory . Wave functions of chaotic quantum systems are discussed in . We also included in this section a paper on the calculation of parametric correlations  and a paper on random waves .

The applications of random matrix theories are by far not exhausted by the above mentioned sections. In the last section we have included random matrix theories that are not included in the previous section. We have included a review on embedded random matrix ensembles by Benet and Weidenmueller , and papers on distance matrices , the Penner model , the Wishart ensemble  and on critical random matrix models .

Acknowledgments

This work was partially supported by the Australian Research Council (P.J.F.), the Royal Society Dorothy Hodgkin Research Fellowship scheme (N.S.), and the US DOE grant DE-FG-88ER40388 (J.J.M.V.).


text:
# On the Nature of X-Ray Flashes

## 1. Introduction

The recently identified X-ray flashes (Kippen et al, 2002) have non thermal spectra with $`E_\mathrm{p}<50`$ keV and weak gamma-ray fluxes ($`\mathrm{}<0.2`$ ph.cm-2.s-1 in the 50 – 300 keV energy range). They were discovered by Beppo-SAX and are currently observed by HETE 2. They can be normal GRBs at very high redshifts (which however does not seem compatible with their duration distribution similar to that of long bursts) or viewed off-axis (Yamasaki et al, 2002) but their soft spectra can also have an intrinsic origin. We have tested this possibility in the context of the internal shock model where it is supposed that a central engine generates a relativistic wind with a non uniform distribution of the Lorentz factor.

## 2. A toy model for internal shocks

The evolution of this relativistic wind should be followed with a fully hydrodynamical calculation. This can be done (Daigne and Mochkovitch, 2000) but requires large amounts of computing time which prevents to consider a large number of cases and fully explore the parameter space. These detailed calculations have however shown that a simplified approach where the wind is represented by many shells which interact by direct collisions only (all pressure waves being neglected) can also produce satisfactory results (Kobayashi, Piran and Sari, 1997; Daigne and Mochkovitch, 1998). Going a step further we have developed for this study a toy model where internal shocks are limited to the collision of two shells of equal mass $`m`$. Shell 2 (Lorentz factor $`\mathrm{\Gamma }_2`$) is generated a time $`\tau `$ after shell 1 (Lorentz factor $`\mathrm{\Gamma }_1<\mathrm{\Gamma }_2`$). The average power injected into the wind in this two shell approximation is given by

$$\dot{E}=\frac{mc^2}{\tau }(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2)=\dot{M}\overline{\mathrm{\Gamma }}c^2$$
(1)
where $`\dot{M}=2m/\tau `$ and $`\overline{\mathrm{\Gamma }}=\frac{1}{2}(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2)`$ are the average mass loss rate and Lorentz factor. Shell 2 will catch up with shell 1 at the shock radius $`r_\mathrm{s}=2c\tau \frac{\mathrm{\Gamma }_1^2\mathrm{\Gamma }_2^2}{\mathrm{\Gamma }_2^2\mathrm{\Gamma }_1^2}`$. The two shells merge and the energy dissipated in the collision $`E_{\mathrm{diss}}`$ is radiated with a characteristic broken power law spectrum. If the synchrotron process is responsible for the emission, the peak energy (maximum of $`\nu F_\nu `$) is

$$E_\mathrm{p}E_{\mathrm{syn}}\mathrm{\Gamma }_\mathrm{s}B\mathrm{\Gamma }_e^2\mathrm{\Gamma }_\mathrm{s}\rho _\mathrm{s}^xϵ_\mathrm{s}^y$$
(2)
In the last term $`\rho _\mathrm{s}`$ and $`ϵ_\mathrm{s}`$ are respectively the typical density and dissipated energy (per unit mass) in the comoving frame of the shocked material. The standard equipartition assumptions correspond to $`x=1/2`$ and $`y=5/2`$ but we consider below the possibility that $`x`$ and $`y`$ can have different values. This may be for example the case if the equipartition parameters $`\alpha _\mathrm{B}`$ and $`\alpha _e`$ which fix the fraction of dissipated energy injected into the magnetic field and the relativistic electrons are not constant but depend on $`\rho _\mathrm{s}`$ or/and $`ϵ_\mathrm{s}`$.

The physical properties of the shocked material $`\mathrm{\Gamma }_\mathrm{s}`$, $`\rho _\mathrm{s}`$ and $`ϵ_\mathrm{s}`$ can be directly related to the wind parameters $`\tau `$, $`\dot{E}`$, $`\overline{\mathrm{\Gamma }}`$ and $`\kappa =\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$ so that eq.(2) can be expressed as

$$E_\mathrm{p}\frac{\dot{E}^x\phi _{xy}(\kappa )}{\tau ^{2x}\overline{\mathrm{\Gamma }}^{6x1}}$$
(3)
where $`\phi _{xy}(\kappa )`$ is an increasing function of $`\kappa `$ for all reasonable values of $`x`$ and $`y`$. In spite of the simplicity of the two shell approximation eq.(3) predicts a duration-hardness relation and a hardness-intensity correlation (HIC) as observed in real bursts. Another interesting (and surprising) consequence of eq.(3) is that $`E_\mathrm{p}`$ is a decreasing function of $`\overline{\mathrm{\Gamma }}`$ as long as $`x>1/6`$. . Soft bursts are obtained with “clean fireballs” (large $`\overline{\mathrm{\Gamma }}`$) while winds with larger baryon loads (“dirty fireballs”) lead to a harder spectrum because internal shocks occur closer to the source at a higher density.

From the values of $`E_{\mathrm{diss}}`$, $`\tau `$, $`E_{\mathrm{break}}`$ and assuming a Band spectrum with low and high energy indices $`\alpha =1`$ and $`\beta =2.5`$ it is possible to estimate (for a given redshift) the average flux or count rate in any spectral band. The simplicity of the two shell approximation then allows to construct a large number of synthetic bursts to check if XRFs can be formed for some specific choice of the wind parameters.

## 3. A statistical approach

The redshift $`z`$ and the four wind parameters $`\tau `$, $`\dot{E}`$, $`\overline{\mathrm{\Gamma }}`$ and $`\kappa `$ have been obtained in a statistical way for each synthetic burst. If long GRBs (and XRFs) are related to the explosive death of some special class of massive stars their birth rate is directly proportional to the star formation rate $`\psi `$ and their distribution in redshift can be deduced from $`\psi (z)`$ (Porciani and Madau, 2001). The distribution of the observed duration $`t_{90}`$ for long bursts is log-normal with a maximum at $`t_{90}20`$ s. We therefore also take a log-normal distribution for $`\tau `$ (with a maximum at $`\tau _{\mathrm{max}}=10`$ s) assuming an average burst redshift $`z1`$. The last three wind parameters $`\dot{E}`$, $`\overline{\mathrm{\Gamma }}`$ and $`\kappa `$ are very poorly constrained and we adopt for them uniform distributions: between 50 and 53 for $`\mathrm{Log}\dot{E}`$, 100 and 500 for $`\overline{\mathrm{\Gamma }}`$, 0 and 1 for $`\mathrm{Log}\kappa `$.

We performed a first numerical experiment with $`10^6`$ synthetic bursts to obtain their $`E_\mathrm{p}`$ distribution. The results are shown in Fig.1 for two choices of $`x`$ and $`y`$: $`x=1/2`$ and $`y=5/2`$ i.e. synchrotron emission with standard equipartition assumptions and $`x=y=1/4`$. This last case was already considered by Daigne and Mochkovitch (2003) who have shown that it leads to very good fits of the temporal and spectral properties of individual GRB pulses. The dashed line in Fig.1 represents the distribution of $`E_\mathrm{p}`$ for all $`10^6`$ bursts while the full line corresponds to the sub-group of bursts which would have been detected by BATSE (a threshold of 0.2 ph.cm-2.s-1 in the 50 – 300 keV energy band was assumed). It can be seen that for $`x=1/2`$ and $`y=5/2`$ the distribution of $`E_\mathrm{p}`$ is much too wide while the agreement with the observations is excellent for $`x=y=1/4`$. Figure 1 also shows that a large fraction of events with $`E_\mathrm{p}<100`$ keV is not detected by BATSE but may contribute to the population of XRFs seen by Beppo-SAX and HETE 2.

To test this possibility we have generated (assuming $`x=y=1/4`$) 1000 events for which we compute both $`E_\mathrm{p}`$ and the photon flux in the 50 – 300 keV band. The results are compared to the BATSE + Beppo-SAX data (Kippen et al, 2002) in Fig.2. Small dots represent BATSE bursts while large dots are events which are not detected by BATSE but have an X-ray flux (2 – 10 keV) larger than 1 ph.cm-2.s-1. The similarity between the two diagrams clearly shows that XRFs can be produced in the context of the internal shock model.

We then obtained the distributions of $`z`$, $`\tau `$, $`\dot{E}`$, $`\overline{\mathrm{\Gamma }}`$ and $`\kappa `$ for synthetic XRFs. Their redshift distribution is nearly identical to that of GRBs and our simulation therefore confirms that XRFs are not GRBs at very high redshifts. The duration and injected power distributions are also not very different between XRFs and GRBs and it finally appears that a reduction of the contrast $`\kappa =\mathrm{\Gamma }_2/\mathrm{\Gamma }_1`$ and an increase of the average Lorentz factor $`\overline{\mathrm{\Gamma }}`$ are the most efficient ways to produce XRFs. The distribution of $`\kappa `$ has a sharp maximum at about 1.5 in XRFs while it steadily increases (and is negligible for $`\kappa <2`$) in GRBs. There are also five times more XRFs with $`\overline{\mathrm{\Gamma }}=500`$ than with $`\overline{\mathrm{\Gamma }}=100`$.

## 4. Conclusion

Our simulations have shown that XRFs can be produced by internal shocks if the Lorentz factor in the relativistic wind has variations of small amplitude only and a high average value. XRFs are not the result of a low injected power $`\dot{E}`$ but are weak and soft mainly due to a small contrast of the Lorentz factor which leads to a reduced efficiency of the dissipation in shocks.

## References

Daigne, F., Mochkovitch, R. 1998, MNRAS, 296, 275

Daigne, F., Mochkovitch, R. 2000, A&A, 358, 1157

Daigne, F., Mochkovitch, R. 2003, submitted to MNRAS

Kippen, R.M. et al 2002, astro-ph/0203114

Kobayashi, S., Piran, T., Sari, R. 1997, ApJ, 490, 92

Porciani, C., Madau, P. 2001, ApJ, 548, 522

Yamasaki, R., Ioka, K., Nakamura, T. 2002, ApJ, 571, L31