format: text

text:
# 1 Introduction

## 1 Introduction

In , it was realized that the gauge theory living on the world volume of a D3-brane probing a toric singularity is sometimes non-uniquely determined. Thus, by considering D3-branes probing non-compact, toric, singular Calabi-Yau manifolds we are lead to more than one gauge theory with the same toric geometry as its moduli space. This phenomenon is the essence of Toric Duality. It is a full equivalence between distinct $`𝒩=1`$, $`d=4`$ gauge theories in the IR limit. Microscopic theories with different matter content and interactions become indistinguishable when we consider the long distance physics they describe. This short note is a summary of recent talks given by the authors which describes the main features of this phenomenon as well as describing a formalism to conveniently compute various features of gauge theories of branes on a class of singular manifolds 111Similar ideas have been used recently in the construction of phenomenological models by wrapping D6-branes on compact, intersecting 3-cycles of Calabi-Yau manifolds ..

The organization of this note is as follows. In Section 2, we present some examples of toric dual theories. Section 3 gives a brief introduction to $`(p,q)`$ webs, which are useful in defining 5d fixed points as well as studying dynamics of 5d gauge theories but also in describing toric varieties and their associated 4d gauge theories. In Section 4, we explain how local mirror symmetry enables the computation of the gauge theory on the world volume of a D3-brane probing a toric singularity using the geometric information encoded in a $`(p,q)`$ web. We also exemplify how the $`(p,q)`$ web machinery can be used to derive toric duals. Based on the mirror Type IIA picture, we show in Section 5 how Picard-Lefschetz monodromy transformations point toward generalizations of Seiberg duality. Section 6 describes the construction of invariants for the singularities under study that generate Diophantine equations encoding the whole set of dual theories.

## 2 Toric duality

The determination of the gauge theory on a D3-brane probing a toric singularity was systematized in , by the development of the Inverse Algorithm. This procedure is based on the realization of toric varieties as partial resolutions of Abelian orbifolds. Using this technique, the cases of cones over the Zeroth Hirzebruch and toric del Pezzo surfaces have been extensively studied. These are toric singularities with a shrinking compact 4-cycle.

In various cases, the resulting gauge theory is non-unique. Let us consider the example of the Zeroth Hirzebruch surface ($`F_0`$) for which the corresponding theories are displayed in Figure 1.

Accompanying the different matter contents summarized in their quivers, these models also display very distinct interactions. These are given by the following superpotentials, which correspond to a sum over a subset of all closed polygons in the quiver

$`W_I`$ $`=`$ $`ϵ_{ij}ϵ_{mn}X_{12}^iX_{23}^mX_{31}^{jn}ϵ_{ij}ϵ_{mn}X_{41}^iX_{23}^mX_{31}^{jn}`$ (2.1)
$`W_{II}`$ $`=`$ $`ϵ_{ij}ϵ_{mn}X_{12}^iX_{23}^mX_{34}^jX_{41}^n`$
Finding the right subset is a difficult task in general and presents a technical challenge. An overall re-scaling of the gauge group ranks by a common factor $`N`$, as well as tracing over gauge indices is understood in the previous quivers and superpotentials.

## 3 $`(p,q)`$ webs

Five dimensional gauge theories can be engineered by 5-brane webs in type IIB string theory . In these constructions, the 5d gauge theories live on the $`4+1`$ common dimensions of the branes. The non-trivial intersections of the branes take place on a 2-dimensional transverse plane. The 5d theories are fully determined by the structure of the webs on this $`(x,y)`$ plane.

Every brane has an associated $`(p,q)`$ charge that dictates its tension

$$T_{p,q}=\left|p+\tau q\right|T_{D_5}$$
(3.2)
and its slope on the $`(x,y)`$ plane

$$\mathrm{\Delta }x:\mathrm{\Delta }y=p:q$$
(3.3)
where $`T_{D5}`$ is the D5-brane tension and $`\tau `$ is the complex scalar of type IIB (which we have set to be equal to $`i`$ in 3.3). Condition 3.3 determines that 8 supercharges are preserved, leading to $`𝒩=1`$ in five dimensions. Furthermore, $`(p,q)`$ charge has to be conserved at each brane intersection

$$\underset{i}{}p_i=\underset{i}{}q_i=0$$
(3.4)
The reader is referred to  for a detailed discussion of $`(p,q)`$ webs, their use in engineering five dimensional theories and explicit examples. Gauge couplings, masses of gauge bosons and quarks, BPS spectrum and monopole tension can be computed straightforwardly from the geometry of the $`(p,q)`$ web .

Alternatively, $`(p,q)`$ webs can be viewed as toric skeletons defining toric varieties  (see also  for applications of this idea along the lines that will be discussed in this note). In this interpretation, 5-branes correspond to the loci of points at which some 1-cycles of the $`T^2`$ fibrations of the toric varieties shrink to zero radius.

## 4 4d theories via local mirror symmetry

Each factor in the product gauge group of the theory on the D-brane world-volume is given by a fractional brane. These are bound states of D3, D5 and D7-branes, sharing four non-compact dimensions. D3-branes are located at points (i.e. 0-cycles) on the Calabi-Yau, while D5 and D7-branes wrap compact 2 and 4 cycles respectively.

The corresponding quiver can be obtained by looking at the mirror Type IIA geometry, in which D3-branes transverse to the original non-compact Calabi-Yau map to D6-branes wrapping a $`T^3`$ . From the homology class of the $`T^3`$

$$\left[T^3\right]=\underset{i=1}{\overset{n}{}}n_iS_in_i\text{ZZ}$$
(4.5)
we can compute the gauge group and matter content of the $`𝒩=1`$, $`d=4`$ gauge theory produced by a wrapped D6-brane

$$G=\underset{i=1}{\overset{n}{}}U\left(n_i\right)I_{ij}=^\mathrm{\#}(S_i.S_j)$$
(4.6)
Each 3-cycle $`S_i`$ wraps a 1-cycle $`C_i`$ of a smooth elliptic fiber that degenerates at a point $`z_i`$. The intersection numbers between the $`S_i`$’s are thus equal to the ones of the $`C_i`$’s which can be easily computed from their $`(p_i,q_i)`$ charges

$${}_{}{}^{\mathrm{\#}}(S_i.S_j)={}_{}{}^{\mathrm{\#}}(C_i.C_j)=det\left(\begin{array}{cc}p_i& q_i\\ p_j& q_j\end{array}\right)$$
(4.7)
Let us study how these concepts come together in the explicit example of $`dP_0`$. The corresponding $`(p,q)`$ web is presented in Figure 2.a, from where we read the following $`(p,q)`$ charges

$$\begin{array}{ccccc}(p_1,q_1)=(1,2)& & (p_2,q_2)=(2,1)& & (p_3,q_3)=(1,1)\end{array}$$
(4.8)
Using 4.7, we compute the following intersection numbers

$$\begin{array}{ccccc}{}_{}{}^{\mathrm{\#}}(C_1.C_2)=3& & {}_{}{}^{\mathrm{\#}}(C_2.C_3)=3& & {}_{}{}^{\mathrm{\#}}(C_3.C_1)=3\end{array}$$
(4.9)
which can be conveniently summarized in the quiver diagram presented in Figure 2.b.

It is interesting to note that in this case, the $`SU\left(3\right)`$ isometry of $`\mathrm{IP}^2`$ appears as a flavor symmetry in the gauge theory, the $`X_{ij}`$ chiral fields transforming in the fundamental representation for each pair of indices $`(i,j)`$. Invariance under this $`SU\left(3\right)`$ fixes the superpotential uniquely

$$W=ϵ_{\alpha \beta \gamma }X_{12}^{\left(\alpha \right)}X_{23}^{\left(\beta \right)}X_{31}^{\left(\gamma \right)}$$
(4.10)
which is the singlet in $`X_{12}X_{23}X_{31}=333=18810`$. This superpotential is also invariant under the $`\text{ZZ}_3`$ cyclic permutations of the nodes $`\left(123\right)`$. Let us demonstrate how one can compute the moduli space of vacua for this model and reproduce the right manifold we started with. The set of gauge invariant operators is given by 27 invariants, $`a^{\alpha \beta \gamma }=X_{12}^{\left(\alpha \right)}X_{23}^{\left(\beta \right)}X_{31}^{\left(\gamma \right)}`$. Using the F-term equations we find that any antisymmetric combination of the indices vanishes. Therefore $`a^{\alpha \beta \gamma }`$ is in the completely symmetric 10 dimensional representation of $`SU\left(3\right)`$. Furthermore, due to this symmetry we find the set of equations

$$\left(a^{\alpha \beta \gamma }\right)^3=a^{\alpha \alpha \alpha }a^{\beta \beta \beta }a^{\gamma \gamma \gamma },$$
(4.11)
which is a set of 7 equations for 10 variables. A quick inspection verifies that this is the set of equations for the orbifold $`C^3/Z_3`$, the manifold we started with.

In fact, $`(p,q)`$ webs are powerful computational tools in deriving toric dual theories . Del Pezzo surfaces are constructed by blowing-up up to eight generic points on $`\mathrm{IP}^2`$. A blow-up corresponds to the replacement of a point by a 2-sphere. Since the toric ($`(p,q)`$ web) representation of a 2-sphere is a segment, the blow-up of a vertex of a given $`(p,q)`$ web corresponds to its replacement by a segment. At the same time the external leg that was originally attached to the blown-up vertex is replaced by a pair of legs, whose $`(p,q)`$ charges are dictated by $`(p,q)`$ charge conservation at the new vertices. Using the $`SL(3,\text{ }\mathrm{C})`$ symmetry of $`\mathrm{IP}^2`$, the positions of up to three generic points can be mapped to vertices of the web. In this way, we see that, starting from the $`(p,q)`$ web for $`dP_0`$ given in Figure 2, we can construct all the toric duals for del Pezzo surfaces up to $`dP_3`$, by blowing-up vertices of the webs in every possible way. We apply this technique to $`dP_1`$ in Figure 3 to obtain the two toric phases of $`dP_2`$. It has been shown in  that for the Zeroth Hirzebruch and del Pezzo surfaces, Toric dual theories are indeed Seiberg duals.

The discussion in this section, in addition to the one in Section 3, allows the translation of four dimensional quantities and processes into five dimensional ones. This approach was pursued in  to link the relocation of blown-up points associated to Toric Duality to the crossing of curves of marginal stability in related five dimensional theories.

## 5 Picard-Lefschetz transformations

We have seen in Section 4 how to exploit local mirror symmetry to compute four dimensional quiver theories from the intersections of 3-cycles in the type IIA mirror picture. This geometric realization of the theories suggests how to generalize Seiberg duality by using Picard-Lefschetz (PL) monodromy transformations.

PL monodromy corresponds to the reordering of vanishing cycles. When moving a vanishing cycle $`S_j`$ around another one $`S_i`$, $`S_j`$ gets a contribution proportional to $`S_i`$, weighted by the mutual intersection number

$$S_jS_j+(S_j.S_i)S_i$$
(5.12)
while $`S_i`$ remains invariant. The 3-cycles $`S_i`$ can be represented by $`[p_i,q_i]`$ 7-branes with wrapping numbers $`n_i`$ and we can recast (5.12) in terms of $`[p,q]`$ charges (Figure 4). PL monodromies corresponds in this language to the motion of a 7-brane across the branch cut of another one

The set of Seiberg dual theories associated to a given singularity can be obtained by starting from a configuration of $`[p,q]`$ 7-branes for the geometry, acting on them with PL monodromy transformations according to the rules in (5.12) and computing the resulting quiver as explained in Section 4. In fact, the group of PL transformations is larger than the one of Seiberg dualities, and there are gauge theories attainable with monodromies that cannot be reached by performing any chain of Seiberg dualities on nodes of the quiver . Since further action with PL transformations results in Seiberg duals, these theories have been named Fractional Seiberg duals.

Figure 5 exhibits three theories related by PL transformations and Seiberg dualities for the Zeroth Hirzebruch surface. In this case, Model 2 is a new fractional Seiberg dual theory that could not have been computed by means of traditional Seiberg duality.

## 6 Diophantine equations from invariant traces

In the previous section, we have described how dual theories can be generated by performing Picard-Lefschetz transformations on the set of degenerate fibers. All along the dualization process there is a set of quantities that remain invariant, thus providing us with a powerful tool in characterizing dual theories. The existence of these invariants was studied in detail in . They are the number of the degenerate fibers, the greatest common divisor of the intersection numbers and the trace of the total monodromy matrix,

$$K=K_{z_1}K_{z_2}\mathrm{}K_{z_n}$$
(6.13)
For the configurations of degenerate fibers under study,

$$TrK=2$$
(6.14)
This equality corresponds, for each geometry, to a Diophantine equation in the intersection numbers that completely encodes all the gauge theories that can be generated by performing Picard-Lefschetz transformations. As an example, let us consider the case of $`dP_0`$. Using the $`(p,q)`$ charges in Figure 2 we arrive at

$$I_{12}^2+I_{23}^2+I_{32}^2I_{12}I_{23}I_{32}=0$$
(6.15)
Since for this specific example $`N_i=3I_{jk}`$, with $`j,ki`$, (6.15) can be turned into an equation for the allowed ranks of the gauge groups

$$N_1^2+N_2^2+N_3^23N_1N_2N_3=0$$
(6.16)
Acknowledgement We would like to thank Bo Feng, Yang-Hui He, Amer Iqbal and Angel Uranga for collaborations in the material presented in this note. A.H. Would also like to thank the organizers of the ”35th International Symposium Ahrenshoop on the Theory of Elementary Particles” for their hospitality. Research supported in part by the CTP and the LNS of MIT and the U.S. Department of Energy under cooperative agreement $`\mathrm{\#}`$DE-FC02-94ER40818. A.H. is also supported by the Reed Fund Award and a DOE OJI ward.


text:
# Simulations of amphiphilic fluids using mesoscale lattice-Boltzmann and lattice-gas methods

## 1 Introduction

The lattice-gas automata (LGA)  and the lattice-Boltzmann (LBE) methods  are relatively new approaches in computational fluid dynamics. Both methods have a mesoscopic character, as opposed to the conventional continuum approach based on numerical solution of the Navier-Stokes equations and the microscopic approach based on molecular dynamics. The key idea behind LGA is to model fluid flows by simplified kinetic equations which describe time-evolution of particles having a discrete set of velocities while moving on a regular lattice. For the purpose of computer efficiency an “exclusion principle” is imposed such that at a given time step only one particle with a given velocity can occupy any lattice site. The computationally demanding tracking of individual molecules is thus avoided and at the same time macroscopic or hydrodynamic effects naturally emerge from mesoscale lattice-gas dynamics, provided that the LGA collision operator possesses the correct and necessary conservation laws and the underlying lattice has sufficient symmetry .

Historically, the lattice-Boltzmann model was developed from lattice-gas automata by replacing the particle occupation numbers (Boolean variables) by single-particle distribution functions (real variables). Furthermore, the LGA collision operator, which is a stochastic many-body collision operator, is massively simplified in LBE methods using the assumption of molecular chaos. Usually a simple BGK approximation  is also employed, which is linear and deterministic, although more complex linearised collision operators may also be used . The spontaneous fluctuations of LGA are thus eliminated in LBE. An important feature of LGA and LBE is the locality of the updating rules which make these methods ideal for parallel processing.

Both methods have found applications in fluid dynamical problems in which more conventional continuum approaches face difficulties, such as multiphase and multicomponent flow and flow in porous media. LBE has also been applied to study colloidal suspensions , reaction-diffusion systems and other complex flow situations . An important area of application of these methods is for modelling amphiphilic fluids which consist of two immiscible phases (such as oil and water), together with an amphiphile (surfactant) species, such as detergent . The study of these fluids is of relevance to a wide variety of industrial, chemical and biological applications and is of great fundamental interest. There is, however, no general consensus on what continuum hydrodynamic equations are most appropriate for describing such complex fluids. Moreover, due to the presence of complicated interfaces which undergo topological change, continuum methods based on numerical solution of Navier-Stokes-type equations face great difficulties in dealing with such fluid mixtures. Atomistic approaches based on molecular dynamics which deal with “real” surfactant molecules are able to deal with such systems but are still computationally too demanding to access large time dynamics involved in many problems of interest, such as self-assembly of micelles and the formation of lamellar and other ordered mesophases .

In the last few years we have pioneered the development of mesoscale lattice-gas  and, more recently, lattice-Boltzmann methods  for the simulation of the non-equilibrium hydrodynamic properties of amphiphilic fluids. The common features of both models is that they are ternary vector models which treat amphiphilic molecules as a separate species with an orientational degree of freedom, and implement coupling between fluid components from the “bottom-up” by introducing interactions between particles. Our LGA studies of amphiphilic fluids, both in two  and three dimensions , demonstrate the ability of the LGA scheme to describe the non-equilibrium dynamics and hydrodynamics of these systems while our recent LBE studies in two-dimensions  show that certain complex fluids phenomenologies, such as microemulsion states can be produced by both schemes. At the same time, these studies have highlighted striking differences in the phenomenology accessible to each method .

The aim of the present article is to provide a comparison between our lattice-gas and lattice-Boltzmann methods for amphiphilic fluids, in order to elucidate the individual strengths and weaknesses of these schemes, in terms of their capability to reproduce in a consistent way the phenomenology of amphiphilic systems and from the point of view of their numerical stability, computational efficiency and speed. We also describe in detail the unified parallelization strategy which we have applied to both algorithms and present the results of extensive benchmark studies of the parallel performance efficiency and speed of the resulting codes, performed on two completely different massively parallel platforms.

The rest of this paper is organised as follows. In Section 2 we provide a short overview of our lattice-gas and lattice-Boltzmann algorithms for amphiphilic fluids and discuss the relation between the two schemes and their similarities and differences. Particular emphasis will be given to the differences in the treatment of intermolecular collisions and their consequences for algorithmic complexity of the methods. In order to be clear in our discussions of both implementational and theoretical issues, we shall refer to the lattice-Boltzmann and lattice gas methods by the initials LGA and LBE, respectively, and to our particular implementations of them by the initials ME3D for our lattice gas code and LB3D for our lattice-Boltzmann code. In Section 4 we describe in detail our parallelization strategy and present results of our benchmark studies of the parallel performance of the LGA and LBE codes on CRAY T3E and Origin 2000 parallel supercomputers. In Section 5 we discuss the high performance visualisation techniques ubiquitous in the understanding of results from both models. We describe our implementation of computational steering for ME3D and LB3D on parallel platforms and discuss our experience with the method. We close this paper in Section 6 with some conclusions and outlook for the future development of the schemes.

## 2 Amphiphilic lattice-gas and lattice-Boltzmann dynamics

In this section, we first summarise our lattice gas and lattice-Boltzmann models, and then make a specific algorithmic comparison between them.

### 2.1 Amphiphilic lattice-gas fluids

Our lattice-gas model  is based on a microscopic, bottom-up approach, where dipolar amphipile particles are included alongside the immiscible oil and water species. The model was developed as an extension of the Rothman and Keller (RK)  model of binary immiscible fluids, including its generalization by Chan and Liang , to take into account the essential characteristics and interactions of amphiphiles. In the RK model  the lattice-gas particles are tagged with two “colours”, to distinguish oil and water. The immiscible fluid effect is incorporated by modifying the LGA collision operator such that the collisions favor outcomes that send particles towards neighbouring sites dominated by other particles of the same colour (under the constraint that mass of of each species and the total momentum of both species together is locally conserved). This is achieved by sampling the collision outcome from the Gibbsian equilibrium corresponding to a local, sitewise Hamiltonian function which depends on the net colour charge associated with a given site

$$q_i(𝐱,t)n_i^R(𝐱,t)n_i^B(𝐱,t),$$
(1)
and the colour flux vector of an outgoing state

$$𝐉(𝐱,t)\underset{i=1}{\overset{b}{}}𝐜_iq_i^{}(𝐱,t),$$
(2)
where the prime denotes a post-collision quantity.

The amphiphilic model is obtained by introducing surfactant particles in the RK lattice gas as a third species. In nature, amphiphilic molecules usually possesses two different fragments, each having an affinity for one of the two immiscible components . In the case of a prototype mixture of immiscible fluids like oil and water, a typical surfactant molecule would be an amphiphile that has a hydrophilic head preferring to be in contact with water molecules, and a hydrophobic tail preferring to be surrounded by oil. This crucial property is taken into account in the lattice gas model by assigning to each surfactant particle at position $`𝐱`$ moving in the direction $`i`$ a “colour dipole vector” $`\sigma _i(𝐱,t)`$ which models the orientational degree of freedom of these molecules. In analogy with electrostatics the colour Hamiltonian of the RK model is then extended to take into account dipolar interactions among surfactant molecules themselves and between surfactant particle and colour charges .

Lattice-gas particles have velocities $`𝐜_i`$, where $`1ib`$, and $`b`$ is the number of velocities per site. A particle emerging from a collision at site $`𝐱`$ and time $`t`$ with velocity $`𝐜_i`$ streams along its velocity vector to site $`𝐱+𝐜_i\mathrm{\Delta }t`$ where it may undergo the next collision. We let $`n_i^\alpha (𝐱,t)\{0,1\}`$ denote the presence ($`1`$) or absence ($`0`$) of a particle of species $`\alpha \{R,B,A\}`$ ($`R`$, $`B`$, $`A`$ denoting red (oil), blue (water) and green (amphiphile) species respectively) with velocity $`𝐜_i`$, at lattice site $`𝐱`$ and time step $`t`$. The collection of all $`n_i^\alpha (𝐱,t)`$ for $`1ib`$ will be called the population state of the site; it is denoted by $`𝐧(𝐱,t)`$. The dipole vector of amphiphilic particles $`\sigma _i`$ is assumed to have a fixed magnitude but its orientation can vary continuously. The collection of the $`b`$ vectors $`𝝈_i(𝐱,t)`$ at a given site $`𝐱`$ and time step $`t`$ is called the orientation state.

The time evolution of the system is an alternation between an streaming or propagation step and a collision step. The propagation step is described mathematically by the replacements

$$n_i^\alpha (𝐱+𝐜_i,t+1)n_i^\alpha (𝐱,t),$$
(3)
$$𝝈_i(𝐱+𝐜_i,t+1)𝝈_i(𝐱,t),$$
(4)
for all $`𝐱`$, $`1ib`$ and $`\alpha \{R,B,A\}`$. That is, particles with velocity $`𝐜_i`$ simply move from point $`𝐱`$ to point $`𝐱+𝐜_i`$ in one time step. In the collision step, the newly arrived particles interact, resulting in new momenta and surfactant orientations. The collisional change in the state at a lattice site $`𝐱`$ is required to conserve the mass of each species present

$$\rho ^\alpha (𝐱,t)\underset{i}{\overset{b}{}}n_i^\alpha (𝐱,t),$$
(5)
as well as the $`D`$-dimensional momentum vector

$$𝐩(𝐱,t)\underset{\alpha }{}\underset{i}{\overset{b}{}}𝐜_in_i^\alpha (𝐱,t),$$
(6)
(where we have assumed for simplicity that the particles all carry unit mass). Thus, the set of population states at each site is partitioned into equivalence classes of population states having the same values of these conserved quantities.

The Hamiltonian for the amphiphilic model is much more complex than the original RK Hamiltonian and is derived and described in detail in . It is given by

$$H(s^{})=𝐉(\alpha 𝐄+\mu 𝐏)+𝝈^{}(ϵ𝐄+\zeta 𝐏)+𝐐:(ϵ𝐑+\zeta 𝐌)+\frac{\delta }{2}𝐯(𝐱,t)^2,$$
(7)
where the various quantities in this expression are defined in Table 1

The mass of the particles is taken as unity, and $`\alpha `$, $`\mu `$, $`ϵ`$, $`\zeta `$ and $`\delta `$ are coupling constants.

In the above equations the term parameterised by $`\alpha `$ models the interaction of colour charges with surrounding colour charges as in the original Rothman-Keller model ; that parameterised by $`\mu `$ describes the interaction of colour charges with surrounding colour dipoles; that parameterised by $`ϵ`$ accounts for the interaction of colour dipoles with surrounding colour charges (alignment of surfactant molecules across oil-water interfaces); and finally that parameterised by $`\zeta `$ describes the interaction of colour dipoles with surrounding colour dipoles (corresponding to interfacial bending energy or “stiffness”). The outgoing states of the collisions are sampled from

$$P(s^{})=\frac{1}{Z}\mathrm{exp}\left[\beta H(s^{})\right],$$
(8)
where $`\beta `$ is an inverse temperature, $`H(s^{})`$ is the energy associated with collision outcome $`s^{}`$, and $`Z`$ is the equivalence-class partition function,

$$Z(\rho ,𝐩,\beta )\underset{𝐧^{}E(\rho ,𝐩)}{}𝑑𝚺^{}\mathrm{exp}\left[\beta H(s^{})\right],$$
(9)
and we have defined the measure on the set of orientational states

$$𝑑𝚺\underset{i=1}{\overset{b}{}}𝑑𝝈_i.$$
(10)
In practice, we sample $`s^{}=(𝐧^{},𝚺^{})`$ by first sampling the postcollision population state $`𝐧^{}`$ from the reduced probability density

$$P\left(𝐧^{}\right)=𝑑𝚺^{}𝐌(s^{}).$$
(11)
We then sample the postcollision orientation state by sampling the $`b`$ orientations $`𝝈_i^{}`$ from each of

$$\pi _i\left(𝝈_i^{}\right)=\underset{ji}{\overset{b}{}}𝑑𝝈_j^{}𝐌(𝐧^{},𝚺^{}).$$
(12)
equation for $`1ib`$; these are independent distributions, so the $`b`$ samples may each be taken without regard for the other outcomes.

In the current implementation of the above algorithm in 3D we use a projected face-centered hypercubic (PFCHC) lattice . The motivation for using this lattice is that it is known to yield isotropic Navier-Stokes behavior for a single-phase fluid .

### 2.2 Amphiphilic lattice-Boltzmann fluids

In close analogy with the above amphiphilic lattice-gas model, our amphiphilic lattice-Boltzmann scheme is obtained by introducing surfactant molecules as a third species in the immiscible two component lattice-Boltzmann scheme of Shan and Chen . In this scheme , coupling between fluid components is achieved by introducing pair-wise interactions between particles, within a mean field approximation. The resulting force acting on particles of component $`\sigma `$ is modelled as

$$𝐅^{\sigma \overline{\sigma }}(𝐱,𝐭)=\psi ^\sigma (𝐱,𝐭)\underset{\overline{\sigma }}{}𝐠_{\sigma \overline{\sigma }}\underset{𝐱^{}}{}𝐆_{\sigma \overline{\sigma }}(𝐱,𝐱^{})\psi ^{\overline{\sigma }}(𝐱^{},𝐭)(𝐱^{}𝐱).$$
(13)
Here $`G_{\sigma \overline{\sigma }}(𝐱,𝐱^{})`$ is an interaction kernel and $`\psi ^\sigma =\psi ^\sigma (n^\sigma (𝐱,t))`$ is a function of density, whose form is chosen empirically to model various types of fluids. Finally, $`g_{\sigma \overline{\sigma }}`$ ($`>0`$ for immiscible fluids) is a coupling constant, whose magnitude controls interfacial tension of the binary immiscible fluid .

In the original Shan-Chen model the above force term is incorporated in an ad hoc fashion by adding an increment

$$\delta 𝐮^\sigma =\frac{𝐅^\sigma }{\rho ^\sigma }\mathrm{\Delta }t$$
(14)
to the velocity $`\stackrel{~}{𝐮}`$ which enters the equilibrium distribution function in the BGK collision operator. A more rigorous inclusion of the force term in LBE models, however, can be achieved by discretization of the continuum Boltzmann equation in the presence of a body force . An advantage of this approach is that, unlike the original Shan-Chen model, the force term does not enter in the equilibrium distribution function, hence making a Chapman-Enskog expansion of the LBE equations straightforward. Furthermore, it can be shown that with the choice $`\psi ^\sigma =n_\sigma `$ the phenomenological Shan-Chen model maps exactly onto a discretised version of the Boltzmann-Vlasov equations of binary interacting fluids .

Following the treatment of surfactant molecules in the amphiphilic LGA scheme the surfactant molecules are assumed to carry a dipole vector (or “director”). Instead of specifying directors for each surfactant particle, however, this degree of freedom is further coarse-grained by introducing a dipole vector field $`𝐝(𝐱,t)`$ which represents the average director of all amphiphilic particles present at site $`𝐱`$ at time $`t`$. Once again, the orientation of $`𝐝`$ is allowed to vary continuously in time but its magnitude is a fixed input of the calculations. Incorporation of dipole-carrying amphiphilic fluid particles results in two fundamentally new types of interaction forces (between amphiphilic and non-amphipilic particles and among amphiphiles themselves) which depend not only on the relative distance between particles but also on the dipolar orientations. These forces are derived from (13 by treating each amphiphilic molecule as a pair of water and oil molecules displaced by a distance $`𝐝(𝐱,t)`$ from each other and performing a Taylor expansion in $`𝐝`$ in the resulting expression for the total force . Assuming that the dipole head and tail have equal and opposite colours $`e=\pm 1`$ and only nearest-neighbour interactions considered, the additional forces are given by

$$𝐅^{\sigma ,s}(𝐱,t)=2\psi ^\sigma (𝐱,t)g_{\sigma s}\underset{i0}{}𝐝(𝐱+𝐜_i\mathrm{\Delta }t)(𝐈\frac{𝐜_i𝐜_i}{c_i^2}D)\psi ^s(𝐱+𝐜_i\mathrm{\Delta }t,t)$$
(15)
$$𝐅^{s,c}(𝐱,t)=2\psi ^s(𝐱,t)𝐝(𝐱,t)\underset{\sigma }{}g_{\sigma s}\underset{i0}{}(𝐈\frac{𝐜_i𝐜_i}{c_i^2}D)\psi ^\sigma (𝐱+𝐜_i,t)$$
(16)
and

$`𝐅^{s,s}(𝐱,t)`$ $`=`$ $`{\displaystyle \frac{4D}{c^2}}g_{ss}\psi ^s(𝐱){\displaystyle \underset{i}{}}𝐝(𝐱+𝐜_i\mathrm{\Delta }t,t)𝐝(𝐱,t):[𝐈{\displaystyle \frac{𝐜_i𝐜_i}{c_i^2}}D]𝐜_i`$
$`+`$ $`𝐝(𝐱+𝐜_i\mathrm{\Delta }t,t)𝐝(𝐱,t)𝐜_i`$
$`+`$ $`𝐝(𝐱,t)𝐝(𝐱+𝐜_i\mathrm{\Delta }t,t)𝐜_i\}\psi ^s(𝐱+𝐜_i+\mathrm{\Delta }t,t)`$
In the above equations $`𝐈`$ is the second-rank unit tensor, $`𝐅^{\sigma ,s}`$ is the force acting on non-amphiphilic particles $`\sigma `$ (water and oil) due to amphiphile dipoles, $`𝐅^{s,c}`$ is the force acting on amphiphilic particles due to all non-amphiphilic particles and $`𝐅^{s,s}`$ is the force among amphiphilic particles themselves. The coupling constants $`g_{\sigma s}`$ and $`g_{ss}`$ determine, respectively, the strength of interaction between water/oil particles and surfactant particles, and among surfactant particles themselves.

The resulting amphiphilic lattice-Boltzmann model is characterised by the following set of coupled equations :

$`f_i^\sigma (𝐱+𝐜_i\mathrm{\Delta }t,t+\mathrm{\Delta }t)f_i^\sigma (𝐱,t)`$ $`=`$ $`\mathrm{\Delta }t{\displaystyle \frac{f_i^\sigma f_i^{\sigma (eq)}}{\lambda _\sigma }}`$ (18)
$`+`$ $`{\displaystyle \underset{\overline{\sigma }}{}}{\displaystyle \underset{j}{}}\mathrm{\Lambda }_{ij}^{\sigma \overline{\sigma }}f_j^{\overline{\sigma }}`$
$`+`$ $`{\displaystyle \underset{j}{}}\mathrm{\Lambda }_{ij}^{\sigma s}f_j^s`$
$`f_i^s(𝐱+𝐜_i\mathrm{\Delta }t,t+\mathrm{\Delta }t)f_i^s(𝐱,t)`$ $`=`$ $`\mathrm{\Delta }t{\displaystyle \frac{f_i^sf_i^{s(eq)}}{\lambda _s}}`$
$`+`$ $`{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{j}{}}\mathrm{\Lambda }_{ij}^{s\sigma }f_j^\sigma `$
$`+`$ $`{\displaystyle \underset{j}{}}\mathrm{\Lambda }_{ij}^{ss}f_j^s`$
$$𝐝(𝐱,t+\mathrm{\Delta }t)\overline{𝐝}(𝐱,t)=\mathrm{\Delta }t\frac{\overline{𝐝}(𝐱,t)𝐝^{(eq)}(𝐱,t)}{\lambda _d}.$$
(20)
The first two equations describe time evolution of discrete velocity distribution functions $`f_i^\sigma `$ and $`f_i^s`$ belonging to component $`\sigma `$ ($`\sigma `$ oil, water) and surfactant (s), respectively and the third equation describes time evolution of surfactant dipoles $`𝐝(x,t)`$. In the above equations $`f_i^{\sigma (eq)}`$, $`f_i^{s(eq)}`$ and $`𝐝^{(eq)}`$ are suitably chosen local equilibrium distribution functions , $`\lambda _\sigma `$, $`\lambda _s`$ and $`\lambda _d`$ are relaxation times and $`\overline{d}`$ is the average dipole at site $`𝐱`$ prior to collisions. The first term on the right-hand sides of Eqns. (18) and (2.2) is the standard BGK collision operator . The terms $`\mathrm{\Lambda }_{ij}^{\sigma \overline{\sigma }}`$, $`\mathrm{\Lambda }_{ij}^{s\sigma }`$, $`\mathrm{\Lambda }_{ij}^{s\sigma }`$ and $`\mathrm{\Lambda }_{ij}^{ss}`$ are matrix elements of collision operators which result from mean-field interactions among different fluid components . For Example $`\mathrm{\Lambda }_{ij}^{\sigma \overline{\sigma }}`$ is given by:

$$\mathrm{\Lambda }_{ij}^{\sigma \overline{\sigma }}=\omega _i\left[\frac{1}{c_s^2}(\delta _{\sigma \overline{\sigma }}𝐜_i\alpha _{\sigma \overline{\sigma }}𝐜_j)+\alpha _{\sigma ,\overline{\sigma }}\frac{𝐜_i.𝐜_𝐣}{c_s^4}𝐜_i\right].𝐚_\sigma \mathrm{\Delta }t$$
(21)
with

$$\alpha _{\sigma \overline{\sigma }}=\frac{n^\sigma }{n^{\overline{\sigma }}}\times \frac{\frac{\rho ^{\overline{\sigma }}}{\tau _{\overline{\sigma }}}}{_{\overline{\sigma }}\frac{\rho ^{\overline{\sigma }}}{\tau _{\overline{\sigma }}}}.$$
(22)
Hydrodynamic quantities such as number densities $`n^{\sigma ,s}`$ and velocities $`𝐮^{\sigma ,s}`$ of each fluid component are obtained from velocity moments of the corresponding distribution functions. The kinematic viscosities of each fluid component are controlled by the corresponding relaxation time $`\lambda _\sigma `$ and $`\lambda _s`$ .

### 2.3 Algorithmic comparison

Each time step of our LGA and LBE algorithms consists of the following substeps: streaming of particles or distribution functions to adjacent sites, computation of interaction fields, fluxes and forces, and local computation of the post-collision state at each site.

In the streaming step, the particle occupation numbers or distribution functions are updated according to (3). LGA particles carry their directors to the neighbouring site during this step and the average director of LBE particles is updated using

$$n^s(𝐱,t)\overline{𝐝}(𝐱,t)=\underset{i}{}f_i^s(𝐱𝐜_i\mathrm{\Delta }t)𝐝(𝐱𝐜_i\mathrm{\Delta }t,t).$$
(23)
Once the streaming substep is completed, the new fields $`S`$, $`𝐄`$, $`𝐏`$, and $`𝐑`$ and $`𝐏`$ (in LGA) and $`𝐅^{\sigma \overline{\sigma }}`$, $`𝐅^{\sigma s}`$, $`𝐅^{sc}`$, $`𝐅^{ss}`$ (in LBE) are computed at each site.

The collision substep of the algorithms differs greatly. The LBE collision step requires only the calculation of the equilibrium distribution functions $`f_i^{\sigma ,(eq)}`$,$`f_i^{s,(eq)}`$ $`d_i^{\sigma ,(eq)}`$ for each species and each velocity vector. Thus in going from a two-dimensional lattice (with typically $`9`$ velocity vectors) to a three-dimensional lattice (with typically $`19`$ to $`23`$ velocity directions) the computational time spent in the collision substep (per site) increases by less than a factor of three. The LGA collisions, however, are truly particulate and involve many-body collisions between all particles entering a site. To perform the LGA collision substep, a list of all possible states that has been sorted according to their equivalence class is recomputed and stored. In 3D the storage of the full list would require a total of $`2^{52}`$ bits of memory and is not feasible on any existing computer. A method for shortening this list was described in  which drastically reduces the required memory to $`12`$ MB. The lookup table accepts the current population state $`𝐧`$ as input and returns a pointer to the initial position and the number of elements of the equivalent class $`E(\rho ,𝐩)`$ in the table of stored states . Next each site loops over this set of allowed postcollision states and computes the required fluxes and fields shown in table 1. Note that unlike the field quantities, which depend only on the precollision state and are calculated once for each site, the fluxes must be calculated for each possible outgoing state. These are then used to compute the Gibbsian probability densities from Eqn. (11) from which a final state $`𝐧^{}`$ is sampled. Once $`𝐧^{}`$ is known, the postcollision orientation states of surfactant particles are sampled from Eqn. (12).

### 2.4 Numerical stability

In conventional fluid dynamics (CFD) much attention is devoted to analysis of the numerical stability of Navier Stokes solvers. CFD algorithms typically become numerically unstable for critical values of the time step and/or grid size. Such considerations become paramount when pursuing the highest possible Reynolds number simulation, as, for example, when studying turbulence. When the FHP lattice gas was originally introduced one feature of particular interest was the unconditional numerical stability of the model. The introduction of our modifications to enable simulation of amphiphilic fluids retain this property. The additional floating point calculations required for computation of the local Hamiltonian are sufficiently straightforward that they do not introduce additional stability considerations.

On moving to the lattice-Boltzmann model, this appealing feature of the lattice gas is lost. The BGK approximation for the collision operator introduces hydrodynamic modes which may cause instabilities. In addition to these hydrodynamic modes, the force term in our discrete Boltzmann equation may become so large that it produces negative values of the distribution function. Such negative values are unphysical and cause the code to become unstable. The force term depends on the composition of the system, and so these instabilities are in principle unpredictable in nature. In practice, in extensive $`2`$ dimensional studies we have found that choices of the LB3D parameters which allow the code to run at all will allow the code to run for in excess of $`2\times 10^5`$ time steps. Such long term stability is much better than that reported for a lattice-Boltzmann model based on free energies, which is currently believed to be unconditionally unstable (the code will become unstable at sufficiently long times for any parameter choice) .

## 3 Non-equilibrium dynamics and self-assembly of amphiphilic fluids

The non-equilibrium behaviour of the LGA and LBE models described above has been previously reported in the literature . In this section we review only those aspects of the models behaviour which yield informative comparisons.

Under normal temperature and pressure conditions, water and oil are immiscible. The addition of amphiphile can bring them together with the formation of a wealth of complex structures whose construction is driven by the properties of amphiphiles. The energy of the amphiphile is lowest when it can find or create surfaces between oil and water at which it can adsorb. Lyotropic phases such as hexagonal, lamellar and cubic phases occur when the concentration of surfactant, relative to that of oil and water, is high enough to arrange the interfaces into structures with long-range order . If the concentration is not high enough to bring about the formation of lyotropic phases but more than enough to overcome the tendency of oil and water to phase separate, then an equilibrium fluid phase can be formed in which the water and oil are solubilised, with surfactant residing predominantly at the interfaces between these two components. Such a fluid is called a microemulsion . Self-assembly also occurs in two component water-surfactant systems, where depending on the concentration of surfactant, disordered micellar aggregates and ordered phases can form .

The above-described phenomenology of amphiphilic fluids is extremely rich and the occurrence of various ordered and disordered phases depends not only on the relative concentrations of water, oil and surfactant, and on temperature, but also on specific molecular properties of surfactant, such as size of the headgroups and ionic strength, inter alia . We do not expect, therefore, that minimal models of amphiphiles, like ours, could provide a complete description of all this phenomenology within the parameter space. However, we do expect that phenomena which depend critically on kinetic fluctuations, nucleation and other particle discreteness effects (such as formation of micelles) are better described by our LGA model than by our LBE model. On the other hand the lack of fluctuations in the LBE model means that self-assembly of structures with long-range order, such as lamellar and cubic phases , are more likely to be found within this model.

In order to demonstrate the above ideas we select a set of parameters for each model such that for a given concentration of oil, water and surfactant both models have as their final equilibrium state the oil-in-water droplet microemulsion phase. This phase is obtained when the concentration of oil in the system is lower than that of water and consists of finely divided spherical regions of oil, with stabilizing monolayers of surfactant surrounding them, embedded within a continuously connected water background. Starting from this common structure we then investigate the phase behavior of each model when the concentration of oil is set equal to zero. Both our LGA and LBE simulations reported below are performed on a $`64^3`$ PFCHC lattice with periodic boundary conditions.

### 3.1 Lattice-Boltzmann simulations

We used the following set of canonical parameters throughout our LBE simulations (the time step $`\mathrm{\Delta }t`$ is set to $`1`$ throughout) $`g_{\sigma \overline{\sigma }}=0.05`$, $`g_{\sigma s}=0.01`$, $`g_{ss}=0.01`$, $`\tau _\sigma =\tau _s=1`$, $`\tau _d=2`$, $`m^\sigma =1`$, $`m^s=1`$. In Figure 1 we show the initial and final states of a simulation in which the average concentrations of water, oil and surfactant are set at $`0.5`$, $`0.25`$ an $`0.25`$, respectively. It can be seen that the introduction of surfactant overcomes the tendency of oil and water to phase separate, resulting in a final equilibrium state which consists of droplets of oil in water. It can also be seen that in the final state the surfactant particles reside predominantly at the oil-water interfaces. Next we performed a simulation in which oil is totally removed from the system. As figure 2 shows, in the absence of oil-water interactions, the tendency of amphiphile to cause interfaces to self-assemble produces a lamellar structure with long range order to form from an initially homogeneous mixture.

### 3.2 Lattice Gas simulations

The following set of canonical parameters as defined in equation 7 are used consistently throughout our LGA simulations: $`\alpha =1.0,ϵ=2.0,\mu =0.75,\zeta =0.5`$. Lattice gas simulations at an oil-water-surfactant composition $`0.25`$, $`0.50`$ and $`0.25`$, identical to that used in the lattice-Boltzmann simulations described above, do not lead to the same final state. Whereas the lattice-Boltzmann morphology shown in Figure 1 is that of an ordered array of spherical droplets, the lattice gas simulations lead to a phase separated final state. This difference in behaviour is due to the particulate fluctuations inherent in the lattice gas but wholly absent in lattice-Boltzmann models. An array of spherical droplets, once formed, has no hydrodynamic mechanism of coarsening, and so the domain structure must coarsen by droplet coalescence. This coalescence occurs due to Brownian motion of the droplets. Lattice-Boltzmann simulations are inherently fluctuationless, and without the ad hoc insertion of noise sources are well known to lack such Brownian motion effects .

In Figure 3 we show the behaviour of a lattice gas simulation which reproduces a droplet microemulsion phase. The results shown are from two simulations, one in which the average concentrations of water, oil and surfactant are set at $`0.95`$, $`0.05`$ an $`0.0`$, respectively, and one in which the average concentrations of water, oil and surfactant are set at $`0.50`$, $`0.05`$ an $`0.25`$, respectively. It can be seen that the introduction of surfactant overcomes the tendency of oil and water to phase separate. However, unlike the lattice-Boltzmann case the final state does not consist of spherical droplets of oil in water. Instead, the oil forms elongated structures, which we refer to as swollen wormlike micelles. The transition from spherical to non-spherical structures in the presence of surfactant is indicative of a balance between the surface tension and stiffness of the interface with surfactant adsorbed to it.

Next we performed a simulation in which oil is totally removed from the system, and the average concentrations of water, oil and surfactant are set at $`0.75`$, $`0.05`$ an $`0.25`$ respectively. Figure 4 shows that in the early stages of the simulation the surfactant forms a network of wormlike micelles. Later in the simulation the formation of a bilayer (albeit one containing several defects, or holes) is apparent. Thus the final equilibrium state for this system in the lattice-gas is much more dynamic, with surfactant particles forming and reforming a range of aggregated structures.

## 4 Parallelization method and performance

The CPU time and memory requirements of ME3D and LB3D scale as $`M\times N^D`$ where $`M`$ is the number of discrete velocity vectors, $`N`$ is the linear size of the system and $`D`$ is spatial dimension. In two dimensions it is possible to study adequate system sizes and time scales using standard workstations. In three dimensions, however, the serial algorithm quickly becomes prohibitive in terms of computer memory and CPU time even for moderately sized systems. Fortunately, an important feature of the methods is their intrinsically parallel structure: using a unified strategy we have implemented parallel version of both algorithms which have allowed us to perform large-scale $`3D`$ studies on massively parallel platforms. In this section we describe our parallelization method and the results of benchmark studies of the two parallel codes.

### 4.1 Parallelization strategy

Parallelization was performed utilizing the single program multiple data (SPMD) model, wherein the same program is executed on all processors, and by means of a domain decomposition strategy. This strategy is very suitable for grid-based and semi-local algorithms like lattice-gas, lattice-Boltzmann, finite-difference and finite-element methods. The underlying 3D lattice is partitioned into spatial sub-domains and each sub-domain is assigned to one processor. Such a decomposition may be done in one-dimension (slices), two-dimensions (rods) or three-dimensions (boxes). We have used box decomposition which offers maximum flexibility and also results in a minimum surface:volume ratio of sub-domains, hence minimizing communication overheads. Each processor is responsible for the particles within its sub-domain and performs exactly the same operations on these particles. Two rounds of communication between neighbouring sub-domains are required: at the propagation step, where particles on a border node can move to a lattice point in the sub-domain of a neighbouring processor, and in evaluating the forces, fields and fluxes. By using a ghost layer of lattice points around each sub-domain, the propagation and collision steps can be isolated from the communication step. Before the propagation step is carried out the values at the border grid points are sent to the ghost layers of the neighbouring processor and after the propagation step an additional round of communication is performed to update the ghost layers. This additional round of communication is required because of the presence of non-local interactions in the model whose computation requires (the updated) occupation numbers/distribution functions at neighbouring sites. Figure 5 shows schematically one step of the parallel algorithm. In the current implementation of the codes we have only considered nearest-neighbour interactions. However, an efficient generalization to long-range interactions of the parallel algorithm, which makes use of parallel fast-Fourier transform (FFT) algorithms, is possible and has currently been implemented in $`2D`$ .

The parallel codes are written in standard FORTRAN90, and make use of a number of features of that language that are object-oriented in spirit. Essentially all of the subprograms used in the codes are encapsulated, in that they accept arrays of any size specified at execute time. Many of them are also overloaded in that they can accept arguments of a variety of different ranks and types, and do the right thing in each case. The codes also make use of user-defined types to represent the data. The codes utilise the standard message passing interface MPI for synchronization and communication between processors . Most of the MPI calls, however, have been embedded in wrapper functions to allow relative ease of porting the code to any other message passing protocol if desired. To port to another parallel language, the MPI calls can simply be replaced by their equivalents in the alternative protocol.

### 4.2 Parallel performance

The parallel performance of the codes was benchmarked on two parallel platforms with totally different architectures: the CRAY T3E 1200E and the SGI O2000. The Cray T3E system is a distributed memory Massively Parallel Processor (MPP) machines (Dec Alpha EV6 processors running at 600 MHz, with $`256`$ Mbytes of memory each), whereas the SGI O2000 machines are Shared Memory Parallel (SMP) machines with cache-coherent Non-Uniform Memory Architecture (ccNUMA) (192 MIPS R10000 processors running at 195 MHz, with 768 MB of high-speed cache memory). All benchmark simulations were performed on $`64^3`$ systems using the PFCHC lattice with coordination number $`24`$ for $`500`$ time steps. The ME3D lattice vectors are the $`24`$ nearest-neighbours of the site (0,0,0,0,0), projected to $`3D`$ plus two rest particles. In LB3D one rest distribution is used.

Fluid densities at each lattice site and surfactant dipole moments were written out periodically at every $`50`$ time steps of the simulations. The presence of surfactant interactions and dynamics significantly increases the memory requirements and execution time of both algorithms. Since in ME3D the occupation numbers are boolean operators the absence of surfactant at a given site can be easily detected and evaluation of the surfactant-dependent portions of forces and fluxes at such sites omitted. This feature is built into the ME3D code such that the overall computation time of the algorithm depends on the average concentration of surfactant in the fluid, and is greatly reduced when little or no surfactant is present in the system. Since the LB3D “occupation numbers” are real variables, implementing a similar feature in the LB3D algorithm is less straightforward and was not attempted. Instead, we implemented two parallel versions of the LB3D algorithm, one for binary immiscible systems and one for ternary amphiphilic systems. Parallel performance of both binary and amphiphilic version of the two algorithms was benchmarked on the distributed memory CRAY T3E 1200E massively parallel platform. To make a fair comparison possible between the performance of the amphiphilic LB3D and ME3D codes we used the same concentrations of surfactant in benchmarking these codes.

In Tables 2 and 3 we compare the performance (number of time steps per minute) of the binary and the ternary codes, respectively. It can be seen that on the Cray T3E the LB3D codes are $`5`$ times faster than the corresponding ME3D codes, regardless of the number of processors used. The only situation in which the ME3D code is able to seriously compete with the LB3D code is for the binary system on the O2000 architecture, where LB3D is only $`10\%`$ faster. The LB3D code is approximately twice as fast for the ternary system on this architecture. This is consistent with our observation above that the ME3D code must compute a large number of components of tensor fluxes and fields in the ternary case. It is likely that the better cache handling properties of the O2000 enhance the performance of the computation of these quantities, and that in the binary case where these calculations are unnecessary the ME3D code is able to maximise its competitiveness with LB3D.

We note that since ME3D is an inherently stochastic algorithm the calculation of physical quantities such as surface tension, diffusion coefficients, etc. with this scheme involves additional ensemble-averaging over a set of initial conditions in order to eliminate statistical noise, which is absent in the lattice-Boltzmann method. This increases the effective computational time of ME3D by, roughly, a factor $`4`$ or $`5`$.

Many analyses of the parallel performance of codes rest on plots of the so-called speedup of the code. We have measured a sequence of speeds $`S_{N_i}`$, measured in time steps per minute for a sequence of parallelisations across$`N_i`$ processors. The speedup is usually defined as:

$$SU_i=\frac{S_{N_i}}{S_{N_0}},$$
(24)
where $`N_0`$ is the smallest number of processors across which the system size chosen will run. Plots of speedup are therefore critically dependent on the choice of $`N_0`$. For $`N_0`$ too small the codes’ performance will be completely dominated by memory bandwidth effects, where the bottleneck is moving data in and out of on-processor memory, irrespective of processor speed and inter-processor bandwidth and latency. As one increases the number of processors the amount of data being moved through the processors’ cache is reduced until some critical unit of data will fit entirely into the on-processor cache. At this point the bottleneck becomes either computational speed, or communication time. In order to separate the parallel performance of our codes from dependence on the unrepresentative behaviour of simulations on small numbers of processors we plot the parallel efficiency:

$$E_i=\frac{S_{N_{i+1}}}{S_{N_i}}\frac{N_i}{N_{i+1}}.$$
(25)
This quantity is equal to one if doubling the number of processors doubles the performance (linear scaling), greater than one if the performance is more than doubled (superlinear scaling), and less than one if the performance is less than doubled (sublinear scaling).

The parallel efficiency of the ME3D and LB3D codes for the binary system is shown in Figure 6. The superlinear scaling of the ME3D code shows that it is significantly influenced by cache effects on the O2K, as previously supposed from its greater competitiveness with LB3D on this platform and system. The LB3D code has very high efficiency and sublinear scaling on the O2K, indicating a computation dominated code uninfluenced by cache effects. The behaviour on the T3E is more complicated. ME3D shows a high efficiency with sublinear scaling, while LB3D has a critical number of processors above which the efficiency improves.

For the ternary system LB3D shows high efficiency and sublinear scaling on both platforms, while ME3D shows a peak of superlinear performance at $`32`$ processors on the O2K, after which the efficiency scales sublinearly, and at $`16`$ processors on the T3E, beyond which the efficiency stays superlinear (Figure 7). This shows that the design of the O2K memory system, with a large 8MB secondary cache, yields significant advantages above 32 CPUs. The T3E only has a 96 kB secondary cache, but also has several other features, such as stream buffers and “E-registers” to improve memory bandwidth and reduce latency. The T3E yields a gradually decreasing advantage with decreasing data per CPU, the codes’ performance still showing possible memory bandwidth effects with 64 CPUs. It is possible that the T3E’s memory system would require specific optimizations to take advantage of its more complicated design.

## 5 High Performance Visualisation and Computational Steering

Technological advances are making parallel computers ever more powerful, and the range of problems that they can tackle is continually expanding. However, a ‘mainframe mentality’ persists in the manner of their use which dates from the nineteen-sixties. They are typically operated as batch-processors sealed off from their users in some way - usually system time must be reserved days in advance. The purpose of this method of operation is to try to maximise the flow of work through a heavily-loaded shared resource. However, this approach overlooks the fact that many jobs ultimately produce negative results.

A common feature of both our LGA and LBE models is the presence of a set of parameters which can be adjusted to produce the desired amphiphilic fluid phenomenology, and/or numerical stability properties. A prerequisite of the use of these methods is the performance of large parameter space searches, which typically involves massively parallel “taskfarms”, in which each processor independently runs a small system with a different set of simulation parameters. Such searches again require the use of massively parallel centralised computing resources.

Were there the capability to visualise the state of running programs in essentially real time many of the aforementioned failures might be detected at an early stage, freeing up the resources for another purpose. In some cases failure might be averted simply by modifying model parameters and feeding them back into the application in progress. Allowing scientists the freedom to exploit their intuition interactively can greatly reduce the computation time required to get results. Computational steering is an efficient method for conducting such searches which allows the user to home in to the desired phenomenology by “steering” his or her way in the parameter space of the models.

Much research effort has been applied to developing computational steering environments for parallel programming (e.g. ). Typically these involve a visualisation component, a communication component, and some sort of model for parallel computation. There are usually tens of thousands of lines of source code that need to be ported to one’s own system. This can be intimidating to a potential user, giving the impression that much work will have to be done. This need not be the case and, indeed, the approach we describe here is bare-bones. We do not provide any parallel programming model, but we support the industry-standard MPI library (and also OpenMP). We do not specifiy any particular visualisation system, but have found that our code fits easily into existing graphics systems, such as AVS. All we provide is the ‘glue’ for connecting the parallel program to visualisation code using three simple subroutines. The advantage of this minimalistic approach is ease of implementation. The effort required to make an existing parallel program steerable is very small, assuming that the operation of the program is well-understood.

Our model for computational steering consists of two elements: a parallel application and a steering agent. The steering agent provides a user interface for visualisation and command entry, and may run on a separate machine from the application. When prompted by the user, the agent sends a signal to the application causing it to pause, at a suitable point in its execution, and then send back a snapshot of its data. The data is then visualised and explored by the user and in due course a restart signal is sent back to the application, possibly accompanied by some modified parameters (figure 8).

The commercial visualisation package AVS (both AVS5 and AVS/Express) was used for all the figures in this paper and for much of our work on steering. As one proceeds to large ($`128^3`$ and $`256^3`$) datasets, generating and rendering complex isosurfaces becomes a computationally demanding task. In particular, one requires that the visualisation can be generated rapidly enough that one can interactively rotate, zoom, and otherwise interrogate it. What this means in practice is that for a given visualisation system and type of data (which determines the complexity of the isosurface) there is a maximum data size above which one cannot usefully perform visualisation.

We found that, for generating and viewing isosurfaces from typical amphiphilic fluid datasets larger than around $`64\times 64\times 64`$, the open-source package VTK (www.kitware.com) was much more responsive than AVS. VTK’s provision of interfaces to several languages made scripting and rapid prototyping of the visualization networks only marginally more onerous than using AVS’s very user friendly drag-and-drop graphical programming interface.

Various parallel solutions to bottlenecks in the visualisation pipeline have been proposed and implemented. In the context of AVS, the rendering bottleneck for large datasets was addressed by the AVS/Express MPE, which implements parallel rendering on multiple graphics pipes. However in our testing of AVS/Express MPE on an SGI Onyx-2 system with 2 pipes we observed no rendering speedup whatsoever. The isosurface (and, in principle, other graphics objects) generation bottleneck was addressed for AVS (and, again in principle, for other visualisation systems) by the VIPAR project . This project implements the parallelisation of graphical object generation routines (“modules” in the context of AVS). This system is currently undergoing active development by the Manchester Visualisation Centre. Parallel visualisation for VTK is also currently under development .

### 5.1 A case study: locating the spinodal point of lattice gas binary immiscible fluids

As a concrete example of the utilisation of our nascent steering techniques we describe the location of a particular point in the lattice-gas models’ phase space using steering techniques. These techniques have also been used to good effect in our lattice-Boltzmann model; however, for reasons of space we confine ourselves to the description of just one example.

An initially homogenous mixture of oil and water below a critical temperature referred to as the spinodal point will spontaneously phase separate. Our lattice-gas model is capable of simulating this behaviour and computational steering was used to locate the spinodal point. Previously such points in the phase diagram of the model have been found using large task-farm parameter searches. Such searches are computationally expensive and time consuming. Additionally, one is limited to a system which will fit on a single processor and so finite size effects may distort the results obtained. In the course of a single simulation lasting under an hour the system’s temperature was repeatedly raised above and lowered below the spinodal point and the behaviour observed by direct visualisation until the desired accuracy was obtained. The sequence of events observed in a cycle of simulation is shown in figure 9. The use of steering meant that the simulation could be stopped when the desired accuracy for the result had been obtained, and represented an enormous saving in both wallclock and CPU time.

For examples such as the above, search algorithms based on some metric measure are the alternative to the more “brute force” taskfarm. For one-dimensional parameter spaces such search algorithms are relatively straightforward. However, the extension of such algorithms to multidimensional parameter spaces is highly non-trivial. By contrast, steering enables the use of intuition to rapidly guide such parameter space searches.

## 6 Conclusions

The conclusions we derive from the work described above fall into three categories. Firstly, what have we learned about amphiphilic systems by studying them with two different models. Secondly, what have we learned about the relative merits in terms of numerical stability and performance of our two models, and thirdly what lessons may we draw from the performance of this work to enable the ease and efficiency of future high performance computation to be improved.

We have demonstrated previously that our two models can both capture a very wide range of amphiphilic behaviour. In this paper we have drawn attention to those areas where our models produce different results. It is in these areas that studying the same system with different mesoscale models yields useful information. It seems likely from our simulation results that the inability of the lattice gas model to create static structures is due to the small ($`1/26`$) fraction of rest particles. In our lattice-Boltzmann model, where no such restriction applies, static lamellar and cubic equilibrium phases have been demonstrated to exist.

Conversely, the absence of noise or stochastic driving forces in our lattice-Boltzmann model means that important modes of domain coarsening are absent. Particularly, Brownian motion effects which arise naturally from the particulate nature of the lattice gas are absent in the model. This leads to the arrest of domain growth for spherical domains which do not have hydrodynamic coarsening mechanisms. The same domains are observed to grow by Brownian collision and coalescence coarsening mechanisms in our lattice-gas model .

The unconditional numerical stability of our lattice-gas model is a feature which greatly aids computational work with this model. The stability problems of the LBE model outlined above mean that, unlike the lattice-gas model, it currently cannot reproduce the full range of amphiphilic behaviour within one set of model parameters, but requires different sets depending on the mesophase required. The hydrodynamic instabilities ubiquitous to all LBE models which use the BGK approximation may be treated with linear stability analysis, and removed by introduction of a linearised collision operator. As our forcing term is separate from the BGK operator such an approach may yield useful results here, and may enable the model to be stable to higher Reynolds numbers, a matter of some significance in the study of binary immiscible fluids. Treatment of the instabilities produced purely due to the multiphase nature of the model would require a less ad hoc approach from the outset, however. Recently, lattice-Boltzmann models based on an appropriate Lyapounov functional with unconditional nonlinear stability have been proposed . Such models could unltimately lead to unconditionally stable multiphase models similar to those described here, although this work is at a rather early stage.

Our results regarding the parallel performance of the models clearly demonstrate the large computational advantage obtained by using the lattice-Boltzmann approximation. The almost order-of-magnitude speedup, even before taking into account the obviation of the need to ensemble-average the results, accounts for the massive takeup of the lattice-Boltzmann method in the mesoscale community. Analysis of the relative influence of cache effects (large for small numbers of processors), computational load (dominant) and communication overhead (acceptably small) show that the computation of the (local) fluxes and (nonlocal) fields make the ternary lattice-gas model sensitive to cache effects. In the binary case a modification of the lattice gas model is possible which makes all computations local. Ghost particles are introduced which effectively transmit information about colour gradients during the propagation step. Our model was modified to include such particles, but no significant speedup was seen, indicating that the bottleneck is indeed in the communication to memory, and not the processor-to-processor communication. In addition, no such modification is possible for the ternary code.

The highly scalable nature of our codes and physical complexity of amphiphilic fluids means that we are able to produce large quantities of complex data. Such data can only be rapidly understood with the aid of high performance visualisation techniques. However, the use of such techniques on gigabytes of data produced in batch simulations and analysed ex post facto may take far longer than the simulations themselves. We have found that closely coupling simulation with visualisation through techniques such as computational steering may be done simply and quickly, and can lead to a large reduction in the wallclock time between conceiving a simulation and gaining an understanding of the results of that simulation.

Our research so far in this area will be extended over the next three years in a UK government (EPSRC) funded e-science testbed called RealityGrid (http://www.realitygrid.org). The mesoscale codes, visualisation and computational steering techniques described in this paper will form the basis of a much more comprehensive and focussed effort to link high performance visualisation, computing and experimental facilities through newly available high bandwidth links.

## Acknowledgements

We are indebted to numerous people and organisations for their support of this work, including Silicon Graphics Incorporated (particularly Bart van Bloemann Waanders, Daron Green and Rob Jenkins), Advanced Visual Systems (particularly Ian Currington), Oxford Supercomputing Centre, the EPSRC E7 High Performance Computing Consortium (particularly Sujata Krishna), the Edinburgh Parallel Computing Centre (particularly Mario Antonioletti), the National Computational Science Alliance in Illinois, USA, Boston University Center for Computational Science, CSAR, and Manchester Visualisation Centre (particularly Jo Leng). Resources for the CSAR T3E and Origin $`2000`$ were allocated under EPSRC grant number GR/M56234, and for the CSAR O$`2000`$ by special arrangement with CSAR. Access to a $`16`$ node SGI Onyx system was provided by the Centre for Computational Science, Queen Mary, University of London.


text:
# Avoided intersections of nodal lines

## I Introduction

The morphology of the nodal sets of wave functions depends crucially on whether the underlying classical dynamics is integrable or chaotic. This was first proposed in miller79  and was followed by the study of various features of the nodal lines, such as the distribution of its curvatures Simmel . Recently, the counting statistics of nodal domains for integrable and chaotic systems were investigated in BGS , and it was shown that in the chaotic case, the statistics follow the predictions derived by assuming that the wave functions are random superpositions of plane waves Berry77 ; bogoschmidt . Local effects due to boundary conditions, and the corresponding modifications of the random wave ensembles were also studied Berry02 ; heller02 .

The interest in the properties of the nodal set is not confined to the physics literature only. Most of the mathematics literature on this subject is concerned with solutions of the Helmholtz equation in the interior of compact domains in $`^2`$ with Dirichlet or Neumann boundary conditions. Courant Courant ; CH  and later Pleijel pleijel  pioneered these studies, and computed an upper bound on the number of nodal domains. Krahn krahn24  provided a lower bound for the area of a nodal domain. Other authors gave estimates on e.g., the length of the nodal sets Don-Fef , and addressed the general properties of the nodal network HON .

Uhlenbeck’s theorem uhlenbeck76  states that the nodal lines of “generic” wave functions do not intersect (this statement will be further discussed in the next section). An important class of exceptions of this rule, are the eigen-functions of separable systems, where the nodal lines form a grid, and the number of intersections equals the number of nodal domains. The purpose of the present work is to provide a quantitative measure of the degree by which nodal lines avoid each other. We shall associate the avoidance range to each avoided crossing, and will compute their distribution for the random waves ensemble. The avoidance range vanishes for an intersection, and therefore, the avoidance range distribution for separable functions is proportional to a Dirac $`\delta `$ at zero. We also show that eigen-functions of classically chaotic systems display an avoidance range distribution which is similar to the one obtained for a random waves ensemble, while for an intermediate system, the avoidance distribution differs substantially from either of the extreme distributions mentioned above. The nodal structure for a few representative wave functions can be seen in figure (1).

The wave function of a rectangle (a) displays a perpendicular grid of nodal lines typical of separable systems. An equilateral triangle (b) is classically integrable, but not separable, and the wave function shows a few crossings and avoided crossings of the nodal lines. The domain in (c) is a Sinai billiard where the classical dynamics is chaotic. With the exception of the boundary crossings, the nodal set displays a characteristic set of avoided crossings.

The paper is organized in the following way: We shall begin the next section by discussing the nodal structure of general solutions of the Helmholtz equation. We shall then identify avoided crossings of the nodal lines, and define the corresponding avoidance range. The avoidance range distribution will be written down explicitly. In section (III) we shall define the random wave ensemble and compute explicitly the expected distribution of the avoidance ranges. The resulting expression will be compared with the avoidance distributions obtained numerically for high lying eigen-functions of a chaotic domain. For Dirichlet problems, the boundary belongs to the nodal set, where nodal intersections occur at a density which is approximately $`2`$ intersections per wavelength. This property cannot be reproduced by the uniform random wave ensembles, which are expected to account for bulk properties. This problem, in the context of the avoided crossing distribution, is discussed at the end of section III.

## II Avoided crossings of nodal lines and the avoidance range

We consider real solutions of the Helmholtz equation in a domain $`\mathrm{\Omega }^2`$ which can be expressed as linear superpositions of regular solutions of the equation in the entire plane. They can be expressed in terms of either plane or cylindrical waves, and the two representations are equivalent. In the plane waves representation the wave function is written as

$$\mathrm{\Psi }(𝐫)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}a_n\mathrm{e}^{i𝐤_n𝐫}$$
(1)
with wave vectors ($`𝐤_n=𝐤_n,|𝐤_n|^2=k^2`$) directed at an angle $`\theta _n`$ and coefficients which ensure that the series (1) converges absolutely (we require $`a_n^{}=a_n`$ to render $`\mathrm{\Psi }(𝐫)`$ real). Using the expansion of a plane wave in cylindrical coordinates we get

$$\mathrm{\Psi }(𝐫)=\underset{l=\mathrm{}}{\overset{\mathrm{}}{}}\alpha _lJ_l(kr)e^{il\theta },$$
(2)
with

$$\alpha _l=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}a_n\mathrm{e}^{il\theta _n},\mathrm{and}\alpha _l^{}=(1)^l\alpha _l.$$
(3)
The infinite sums in both (1) and (2) can be truncated, when one considers domains of finite area. Semiclassical arguments show that the necessary number of terms is $`Lk/\pi `$, where $``$ is the perimeter of the domain. For reasons which will become clear in the sequel, we prefer to use the cylindrical wave representation in the present work.

The expansion (2) refers to a particular choice of the origin. Using Graf’s addition theorem Graf  the origin can be shifted to $`𝐫`$, and the wave function retains its form,

$$\mathrm{\Psi }(𝐫+𝝆)=\underset{m}{}\alpha _m(𝐫)J_m(k\rho )e^{im\varphi },$$
(4)
where

$$\alpha _m(𝐫)=\underset{l}{}\alpha _l(\mathrm{𝟎})J_{lm}(kr)e^{il\theta },$$
(5)
with $`\alpha _l(\mathrm{𝟎})=\alpha _l`$, and the angle $`\varphi `$ is measured from the direction defined by $`𝐫`$. The translated coefficients $`\alpha _m(𝐫)=\beta _m(𝐫)+i\gamma _m(𝐫)`$ are related to the wave-function and its derivatives computed at the point $`𝐫`$:

$`\beta _0(𝐫)`$ $`=`$ $`\mathrm{\Psi }(𝐫)`$ (6)
$`\beta _1(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{k}}\mathrm{\Psi }_r(𝐫);\gamma _1(𝐫)={\displaystyle \frac{1}{kr}}\mathrm{\Psi }_\theta (𝐫)`$
$`\beta _2(𝐫)`$ $`=`$ $`\mathrm{\Psi }(𝐫)+{\displaystyle \frac{2}{k^2}}\mathrm{\Psi }_{rr}(𝐫);\gamma _2(𝐫)={\displaystyle \frac{2}{k^2r^2}}(\mathrm{\Psi }_\theta (𝐫)r\mathrm{\Psi }_{r\theta }(𝐫))`$
$`\mathrm{\Psi }_{\theta \theta }(𝐫)`$ $`=`$ $`kr\beta _1(𝐫){\displaystyle \frac{k^2r^2}{2}}\left(\beta _2(𝐫)+\beta _0(𝐫)\right)`$
In the close vicinity of $`𝐫`$, where $`k\rho <1`$, and to second order in $`k\rho `$

$`\mathrm{\Psi }(𝐫+𝝆)`$ $``$ $`\beta _0(𝐫)\left(1\left({\displaystyle \frac{k\rho }{2}}\right)^2\right)`$
$`+`$ $`\left|\alpha _1(𝐫)\right|\left({\displaystyle \frac{k\rho }{2}}\right)\mathrm{cos}(\varphi +\varphi _1)+{\displaystyle \frac{1}{2}}\left|\alpha _2(𝐫)\right|\left({\displaystyle \frac{k\rho }{2}}\right)^2\mathrm{cos}2(\varphi +\varphi _2),`$
where $`\varphi _l`$ are the phases of $`\alpha _l(𝐫)`$. If $`\beta _0(𝐫)=0`$, $`𝐫`$ is a nodal point. It cannot be an isolated zero since the second term vanishes on the line segment through $`𝐫`$ which is oriented at the direction $`\frac{\pi }{2}\varphi _1`$. Hence, the nodal set consists of lines. Two nodal lines intersect at $`𝐫`$ if both $`\beta _0(𝐫)=0`$ and $`\alpha _1(𝐫)=0`$, while $`\alpha _2(𝐫)0`$. The intersection is perpendicular since $`\mathrm{cos}2(\varphi +\varphi _2)`$ vanishes along two perpendicular lines which intersect at $`𝐫`$. For the time being we shall continue the discussion assuming that $`\alpha _2(𝐫)0`$. The more general case will be commented on at the end of this section.

An avoided crossing occurs at $`𝐫`$ when $`\alpha _1(𝐫)=0`$ and $`|\alpha _2(𝐫)|>|\beta _0(𝐫)|>0`$. In other words, when $`𝐫`$ is a saddle point of the wave function. This can be easily seen by writing the equation of the zero set of (II) in terms of the local coordinates $`𝝆=(\xi ,\eta )`$,

$$1=\xi ^2\left(\frac{\beta _0(𝐫)|\alpha _2(𝐫)|}{\beta _0(𝐫)}\right)+\eta ^2\left(\frac{\beta _0(𝐫)+|\alpha _2(𝐫)|}{\beta _0(𝐫)}\right).$$
(8)
This is a hyperbola (ellipse) if $`|\alpha _2(𝐫)|`$ is larger (smaller) than $`|\beta _0(𝐫)|`$. At an avoided crossing, the scaled distance between the two branches is

$$z(𝐫)kd(𝐫)=4\sqrt{\frac{|\beta _0(𝐫)|}{|\beta _0(𝐫)|+|\alpha _2(𝐫)|}}.$$
(9)
This is the avoidance range associated with the avoided crossing at $`𝐫`$.

A few comments are in order:

i. At a nodal intersection $`z=0`$.

ii. At a saddle point $`|\alpha _2(𝐫)|>|\beta _0(𝐫)|`$, hence $`z<2\sqrt{2}`$.

iii. An equivalent expression for $`z(𝐫)`$ in terms of $`\mathrm{\Psi }(𝐫)`$ and its Cartesian derivatives reads

$$z(𝐫)=4\sqrt{\frac{k^2|\mathrm{\Psi }|}{k^2|\mathrm{\Psi }|+\sqrt{4\mathrm{\Psi }_{xy}^2+(\mathrm{\Psi }_{xx}\mathrm{\Psi }_{yy})^2}}}|_𝐫.$$
(10)
iv. Consider an elliptic critical point of the wave function. In the quadratic approximation the area of the elliptic nodal domain is

$$𝒜=4\pi k^2\sqrt{\frac{\beta _0^2}{\beta _0^2|\alpha _2|^2}}|_𝐫.$$
(11)
This area is always larger than $`4\pi k^2`$. Krahn’s theorem krahn24  gives $`j_{0,1}^2\pi k^2`$ as the maximal lower bound to the area of any nodal domain, where $`j_{0,1}2.405`$ is the first zero of the Bessel function $`J_0(x)`$. The lower bound $`4\pi k^2`$ is smaller but not very far from Krahn’s exact value and thus consistent.

So far we considered the intersections of two nodal lines. However, higher order intersections may occur. In general, if the first non vanishing coefficient at $`𝐫`$ is $`\alpha _q`$, then $`𝐫`$ is a nodal point of order $`q`$, where $`q`$ nodal lines intersect at angles $`\frac{\pi }{q}`$. The higher $`q`$, the more rare are the intersections, since more conditions are to be satisfied by the coefficients. This explains Uhlenbeck’s theorem uhlenbeck76  mention above. From now on we shall discuss the most common intersections with $`q=2`$, and comment about the higher order intersections whenever necessary.

Up to now we discussed individual avoided crossings, and defined the associated avoidance ranges. In the next section we shall consider the distribution of the avoidance ranges of a wave function in the domain of its definition, and compute its mean for random waves ensembles.

## III Avoidance range distributions

The number of the critical points of $`\mathrm{\Psi }(𝐫)`$ in the domain $`\mathrm{\Omega }`$ is given by

$$N_C=_\mathrm{\Omega }rdrd\theta \delta (\mathrm{\Psi }_r(𝐫))\delta \left(\frac{1}{r}\mathrm{\Psi }_\theta (𝐫)\right)|𝒥(𝐫)|,$$
(12)
with $`𝒥(𝐫)=\frac{1}{r}(\mathrm{\Psi }_{rr}(𝐫)\mathrm{\Psi }_{\theta \theta }(𝐫)\mathrm{\Psi }_{r\theta }^2(𝐫))`$ the Jacobian. Using (II) we get

$$N_C=\frac{k^2}{4}_\mathrm{\Omega }rdrd\theta \delta (\beta _1(𝐫))\delta (\gamma _1(𝐫))\left||\alpha _2(𝐫)|^2\beta _0^2(𝐫)\right|.$$
(13)
To count the saddle points, we add the restriction $`|\alpha _2(𝐫)|^2\beta _0^2(𝐫)>0`$ and obtain

$$N_S=\frac{k^2}{4}_\mathrm{\Omega }rdrd\theta \delta (\beta _1(𝐫))\delta (\gamma _1(𝐫))(|\alpha _2(𝐫)|^2\beta _0^2(𝐫))\mathrm{\Theta }(|\alpha _2(𝐫)|^2\beta _0^2(𝐫)).$$
(14)
Combining (14) and (9), the number of avoided crossings with avoidance ranges less than $`z`$ is

$$\stackrel{~}{}(z)=\frac{k^2}{4}_\mathrm{\Omega }rdrd\theta \delta (\beta _1(𝐫))\delta (\gamma _1(𝐫))(|\alpha _2(𝐫)|^2\beta _0^2(𝐫))\mathrm{\Theta }(|\alpha _2(𝐫)|^2\beta _0^2(𝐫))\mathrm{\Theta }(zz(𝐫)),$$
(15)
and the fraction of the total number is $`(z)=\stackrel{~}{}(z)/N_S`$. This counting function, and its associated density $`𝒫(z)\frac{\mathrm{d}(z)}{\mathrm{d}z}`$ are the distributions which characterize the avoided crossings of the nodal set. The range of $`z`$ is $`[0,2\sqrt{2}]`$. $`(z)`$ is normalized such that it takes the value $`1`$ at $`z=2\sqrt{2}`$. One can easily check that an avoided crossing of order $`q`$ is counted in (15) with a multiplicity $`q1`$. The number (including multiplicity) of nodal crossings provides the value of $`(0)`$. The effective multiplicity of nodal crossings of the boundary is reduced by a factor one half.

Numerically computed $`(z)`$ are shown in figure 2 for the three wave functions plotted in figure 1. The computed $`(z)`$ do not take into account the nodal crossings of the boundary. For the separable billiard all the saddle points have zero avoidance range, and the function $`(z)`$ is trivially equal to one in all the $`z`$ range. In the triangular billiard there are still some saddle points with zero avoidance range (nodal crossings), but most of them have a finite avoidance range. For the chaotic billiard, all the calculated saddle points have a finite avoidance range.

### III.1 Avoidance range distribution for the random waves ensemble

One of the main goals of the present work is to show that the properties of the nodal set of chaotic billiards, as detected by the distribution of avoidance ranges, are reproduced by the distributions computed for the random wave ensembles. Because wave functions are subject to boundary conditions, it is expected that the predictions of the isotropic random wave ensemble used e.g. in BGS ; bogoschmidt , are relevant only to the bulk of the domain, and will do poorly in the $`\lambda =(\frac{2\pi }{k})`$ vicinity of the boundary. This was observed and discussed in BGS .

The isotropic random wave ensemble is the ensemble of wave functions (2) where the real parameters $`\beta _l`$ and $`\gamma _l`$ are independent, identically distributed random Gaussian variables with zero mean and unit variance for all $`|l|1`$. Because $`\gamma _0=0`$ the variance of $`\beta _0`$ is twice that of all the others. The ensemble average of a function $`f`$ will be denoted by $`f`$. The local coefficients $`\alpha _l(𝐫)=\beta _l(𝐫)+i\gamma _l(𝐫)`$ were derived from the original ones by a unitary transformation. Hence, they are also independent identically distributed Gaussian variables, and the ensemble averages of the number of critical points (12), the number of saddle points (14) and the mean distribution of avoidance ranges (15) can be computed by considering Gaussian integrations with respect to the variables $`\beta _0,\beta _1,\gamma _1,\beta _2`$ and $`\gamma _2`$. A straightforward integration gives,

$$N_C=\frac{k^2|\mathrm{\Omega }|}{2\pi \sqrt{3}};N_S=\frac{k^2|\mathrm{\Omega }|}{4\pi \sqrt{3}},$$
(16)
where $`|\mathrm{\Omega }|`$ is the area of the domain. The average number of saddle points is a half of the number of critical points. The other half are the points were the wave function has either a minimum or a maximum.

The mean number of saddle points with an avoidance range less than or equal to $`z`$ can be also computed,

$$\stackrel{~}{}(z)=\frac{k^2|\mathrm{\Omega }|}{4\pi }\frac{3z^2(16z^2)^2}{(51264z^2+3z^4)^{3/2}},0<z<2\sqrt{2}.$$
(17)
The normalization gives trivially the mean number of saddle points in the ensemble. Therefore we can define the probability

$$_{r.w.}(z)=\frac{\stackrel{~}{}(z)}{N_S}=\frac{3\sqrt{3}z^2(16z^2)^2}{(51264z^2+3z^4)^{3/2}},0<z<2\sqrt{2},$$
(18)
with a corresponding density

$$𝒫_{r.w.}(z)=\frac{6144\sqrt{3}z(8z^2)(16z^2)}{(51264z^2+3z^4)^{5/2}},0<z<2\sqrt{2}.$$
(19)
These results served to test the conjecture that the bulk properties of the nodal sets of chaotic wave functions are reproduced by the predictions of the isotropic random wave ensemble. For this purpose we computed numerically the first 2400 eigen-functions of the Sinai billiard shown in figure 1.c. For each wave-function the critical points were found by a numerical search. The saddles which correspond to boundary intersections were excluded. The avoidance range (10) was computed for each saddle. In figure 3 we plot the cumulative histogram of the avoidance ranges for the 2000th eigen-state in full line (1102 saddle points), compared with formula (17) in dashed line. Even for a single eigen-state the agreement is very good. Averaging the avoidance range distributions over a group of neighboring eigen-states the numerical histogram and the theoretical curve approximately coincide. The inset in figure 3 shows the differences between the random wave prediction and the mean avoidance range distributions computed for two groups of eigen-states. The difference is small, and shows no systematic deviations.

Another distribution which we compared to the prediction of the random waves ensemble is the number of saddle points. Normalizing the number of saddle points by the prediction of the isotropic random wave ensemble, we observe as $`k`$ increases that the numerical computation approach the predicted value (see figure 4).

The systematic deviation observed at finite $`k`$ is due to boundary effects. Following Berry02  we compute the effect of an infinite straight Dirichlet line, which is reproduced by the wave ensemble

$$\mathrm{\Psi }(𝐫)=2\underset{n=1}{\overset{\mathrm{}}{}}c_n\mathrm{sin}(n\theta )J_n(kr),$$
(20)
where the $`c_n`$ are real coefficients taken as independent random Gaussian distributed variables.

Not entering into the details of the calculation, the computed mean density of saddle points approach the bulk expression (16) as the distance from the Dirichlet line increases. The integrated density in the perpendicular direction shows a global deficiency of saddle points relative to the bulk value. It diverges logarithmically as a function of the distance from the Dirichlet line Berry02 .

We compared the results of the ensemble (20) to the number of saddle points counted for highly exited eigen-states of a chaotic billiard with sufficiently smooth boundaries. The integral along the boundary multiplies the density by the perimeter of the billiard $``$. The perpendicular integral must be truncated in view of the logarithmic divergence mentioned above. A sensible choice of the truncation distance is $`R=\sqrt{|\mathrm{\Omega }|}/2`$. The resulting estimate for the mean number of saddle points is

$$N_S\frac{k^2|\mathrm{\Omega }|k(\sigma _1\mathrm{log}(kR)+\sigma _2)}{4\pi \sqrt{3}}.$$
(21)
with $`\sigma _10.014`$ and $`\sigma _22.0`$. The deficiency of saddle points is explained by the effect of the Dirichlet boundaries that affects the statistics. The dashed line in figure 4 represents the equation (21) and the agreement is definitely improving. However, the domain of low values $`k<500`$ is not in complete agreement. This deficiency can be associated with the corners of the billiard, to its curvature, or to the finiteness of the Dirichlet line that the ensemble (20) can not reproduce. These possible causes will be studied elsewhere.

In conclusion we can say that the properties of the nodal set of chaotic wave functions which were investigated in this work, are very well reproduced by the isotropic random waves ensemble in the semiclassical domain. These findings are consistent with previous works on the subject, and add support to the random waves conjecture.

## IV Acknowledgments

This work was supported by the Minerva Center for non-linear Physics at the Weizmann Institute and by an ISF research grant. AGM acknowledges a post doctoral fellowship from the European Network on Mathematical aspects of Quantum Chaos which support his stay at the Weizmann Institute.