text:
# Spins, charges and currents at Domain Walls in a Quantum Hall Ising Ferromagnet.

## Abstract

We study spin textures in a quantum Hall Ising ferromagnet. Domain walls between ferro and unpolarized states at $`\nu =2`$ are analyzed with a functional theory supported by a microscopic calculation. In a neutral wall, Hartree repulsion prevents the appearance of a fan phase provoked by a negative stiffness. For a charged system, electrons become trapped as solitons at the domain wall. The size and energy of the solitons are determined by both Hartree and spin-orbit interactions. Finally, we discuss how electrical transport takes place through the domain wall.

Recently there is a great interest on the study of spin properties of quantum Hall states. For some filling factors, $`\nu `$, in the integer quantum Hall effect (QHE), a transition from a ferromagnetic (F) to an unpolarized (U) ground state (GS) can be achieved by changing the ratio between the cyclotron ($`\mathrm{}\omega _c`$) and the Zeeman ($`E_Z`$) energiesGiuliani . Experimental evidence of this transition has been addressed recentlyPiazza ; Poortere ; Jungwirth0 ; Daneshvar . In the fractional QHE, transition between F and U states at filling factors $`\nu =2/3`$ and $`2/5`$ can be tuned by varying $`E_z`$ with respect the electron-electron interaction energy. Experimental indication of this transition has been reportedEisenstein ; Engel ; Kukushkin ; Smet ; Freytag ; Eom ; Kronmuler ; Hashimoto  although the order of the transition is not clearChakraborty ; Apalkov ; Murthy1 .

In this work we study a domain wall (DW) separating the U state from the F state, at $`\nu `$=2. At this filling factor, the F state has electrons occupying the Landau levels $`n=0,1`$ with spin up. In the U state, the electrons occupy the Landau level $`n=0`$ with the two spins orientations, i.e., the U state coincides with the singlet (S) state.

At integer $`\nu `$, the phase transition between the uniform GS’s can be described by the functional $`=\alpha m_z+\beta m_z^2`$, where $`𝐦=(𝐦_{},m_z)`$ is a unitary vector field, parallel to an isospin variable which points to the positive (negative) $`z`$ direction when the GS is F (U)Jungwirth01 . The phase transition occurs when $`\alpha `$=0. For $`\nu `$ being an integer greater than or equal to 2, $`\beta `$ is negative, indicating a first order phase transition between the $`m_z`$=1 and the $`m_z`$=-1 states: for this reason these systems are called quantum Hall Ising ferromagnets. The existence of hysteresis in transport experiments is the smoke signal of the occurrence of a first order phase transitionPiazza ; Poortere ; Jungwirth0 ; Daneshvar .

In this work we present the following results:

i) We obtain a functional for describing isospin textures in the system at $`\nu `$=2. Due to the odd parity of the product of the Landau levels wave-functions participating in the F and S phases, we find a negative stiffness for distortions of $`𝐦_{}`$. This stiffness is not able to change the order of the phase transition at $`\nu `$=2.

ii) By integrating out the transverse coordinate, we obtain a one-dimensional functional for describing spin textures at the DW. The adequacy of the functional is established by microscopic calculations. This functional also has a negative stiffness, which could produce a fan spin texture in the transverse isospin component along the DW but, for GaAs quantum wells having widths of a few hundred $`\AA `$ and electron densities from $`10^{11}`$ to $`10^{12}cm^2`$, Hartree repulsion prevents the formation of such topological structure.

iii) When the system is charged, extra electrons get trapped at the DW as topological excitations, solitons, with size and energy controlled by both the Hartree and spin-orbit (SO) interactions. The energy of the soliton controls the transport properties through the DW, the conductance being non-zero only for finite SO coupling. This result solves the problem of the spin conservation; transport through a DW implies a carrier spin flip something that can occur in presence of SO interaction.

Energy functional for isospin textures. The electron states of a two-dimensional electron gas confined in the $`xy`$ plane and a magnetic field applied in the $`z`$-direction, are characterized by the Landau level index $`n`$, the degeneracy index $`X`$ and the spin $`\sigma `$. In the Landau gauge, $`X`$ is the momentum in the $`y`$ direction as well as the orbit center of the $`x`$-part of the wave-function. In this work, the magnetic length $`\mathrm{}`$ and the interaction $`e^2/ϵ\mathrm{}`$ are the units of length and energy respectively. In both the F and the S states, all the $`|n=0,X,\sigma =`$ states are occupied and we consider them as electrically inert, being included in the vacuum. The $`\nu =2`$ states are described by

$$\mathrm{\Psi }=\underset{X}{}(\mathrm{cos}\theta (X)c_{X,}^{}+\mathrm{sin}\theta (X)e^{i\psi (X)}c_{XG,}^{})|0,$$
(1)
where $`|0`$ is the vacuum, $`c^{}`$ are creation operators, $`X`$ runs over all possible states and the isospins $``$ and $``$ represent the states $`n=1,\sigma =`$ and $`n=0,\sigma =`$ respectively. In Eq.(1) we only mix two isospins, since we suppose that $`\mathrm{}\omega _c`$ is large enough for not producing Landau level mixing in the S and F phasesMurthy2 ; Park . Assuming that $`\theta (X)`$ and $`\psi (X)`$ change slowly and $`G`$ is small, the unitary vector field corresponding to the state (1) has the form,

$`m_z(x)`$ $`=`$ $`\mathrm{cos}2\theta (x)`$
$`m_x(x,y)+im_y(x,y)`$ $`=`$ $`\mathrm{sin}2\theta (x)e^{i(\psi (x)+Gy)}.`$ (2)
By computing the expectation value of the energy for the wave-function Eq.(1), we obtain the following energy functional for isospin textures:

$`_{2D}=\alpha {\displaystyle 𝑑𝐫m_z(𝐫)}+\beta {\displaystyle 𝑑𝐫m_z^2(𝐫)}+{\displaystyle \frac{\rho _{}}{2}}{\displaystyle 𝑑𝐫\left(_\mu m_z(𝐫)\right)^2}`$
$`+{\displaystyle \frac{\rho _{}}{2}}{\displaystyle 𝑑𝐫\left(_\mu 𝐦_{}(𝐫)\right)^2}+\widehat{\rho }{\displaystyle 𝑑𝐫\left(_{\mu }^{}{}_{}{}^{2}𝐦_{}(𝐫)\right)^2}+V_H.`$ (3)
The coefficients are,

$`\alpha `$ $`=`$ $`\alpha _1+\left(E_Z{\displaystyle \frac{\mathrm{}\omega _c}{2}}\right){\displaystyle \frac{1}{2\pi }}`$
$`\alpha _1`$ $`=`$ $`{\displaystyle \frac{1}{8\pi }}\left(\mathrm{\Sigma }_{0,0,0,0}\mathrm{\Sigma }_{1,1,1,1}\mathrm{\Sigma }_{1,0,1,0}\right)`$
$`\beta `$ $`=`$ $`{\displaystyle \frac{1}{16\pi }}\left(\mathrm{\Sigma }_{0,0,0,0}+\mathrm{\Sigma }_{1,1,1,1}2\mathrm{\Sigma }_{1,1,0,0}\right)`$
$`\rho _{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\rho ^{0,0}+\rho ^{1,1}\right),\rho _{}=\rho ^{1,0}`$ (4)
with

$`\mathrm{\Sigma }_{n,n_1,n_2,n_3}={\displaystyle \frac{1}{S}}{\displaystyle \underset{𝐪}{}}v(𝐪)F_{n,n_1}(𝐪)F_{n_2,n_3}(𝐪)`$
$`\rho ^{n,n_1}={\displaystyle \frac{1}{2\pi L}}{\displaystyle \underset{𝐪}{}}{\displaystyle \frac{q^2}{4}}v(𝐪)F_{n,n}(𝐪)F_{n_1,n_1}(𝐪)`$
$`\widehat{\rho }={\displaystyle \frac{1}{48\pi L}}{\displaystyle \underset{𝐪}{}}{\displaystyle \frac{q^4}{4}}v(𝐪)F_{n,n}(𝐪)F_{n_1,n_1}(𝐪)`$ (5)
with $`S`$ and $`L`$ being the area and length (along the $`y`$-direction) of the system and $`v(𝐪)`$ the Fourier component of Coulomb interaction. For a strictly two dimensional system, the form factors are: $`F_{0,0}(𝐪)=e^{q^2/4}`$, $`F_{1,1}(𝐪)=(1q^2/2)e^{q^2/4}`$ and $`F_{1,0}(𝐪)=(q_y+iq_x)e^{q^2/4}/\sqrt{2}`$. The coefficients become: $`\alpha _1=3/32\sqrt{2\pi }`$, $`\beta =3/64\sqrt{2\pi }`$, $`\rho _{}=11/128\sqrt{2\pi }`$, $`\rho _{}=1/32\sqrt{2\pi }`$ and $`\widehat{\rho }=0.0035`$.

The novelty in this functional is the negative value of the transversal stiffness, $`\rho _{}`$. In order to control the spatial variation of $`𝐦_{}`$, it is necessary to include in the expansion a higher derivative of $`𝐦_{}`$. $`\rho _{}<0`$ due to the different parity of the $`n=0`$ and the $`n=1`$ Landau level wave-functions; in this way $`\rho ^{1,0}`$ and $`\rho ^{2,1}`$ are negative whereas $`\rho ^{2,0}`$ and $`\rho ^{3,1}`$ are positive. $`\rho _{}<0`$ could produce intermediate helical phases between the F and the S states, however, at $`\nu =2`$ the magnitude of $`\rho _{}`$ is not big enough for this occurrence.

Quantum Hall ferromagnets have the unique property that the topological charge is directly related to the electrical chargeSondhi ; Fertig . Therefore, we include in the functional Eq.(3) a Hartree term, $`V_H`$, representing the interaction between the charge densities $`q(𝐫)=\epsilon _{\mu ,\nu }𝐦(_\mu 𝐦\times _\nu 𝐦)/8\pi `$ associated to the isospin texture. We use a standardSondhi ; Fertig  expresion for $`V_H`$ including the semiconductor dielectric constant and finite width of the quantum well. The F-S degeneracy occurs when $`\alpha `$=0, and the negative sign of $`\beta `$ indicates the first order character of the transition. For $`E_Z`$=0 the phase transition occurs at $`\mathrm{}\omega _c=0.472`$, which corresponds to an electron separation $`r_s=2.12`$. This justifies the use, in Eq.(1), of just the $`n=0`$ and the $`n=1`$ Landau levelsMurthy2 .

Domain wall structure When $`\alpha `$=0, the S and the F states are degenerated, and disorder or finite temperature can produce DW’s separating these GS’s. For studying the structure of a DW, we assume $`\alpha `$=0 and impose to the functional (3) the boundary conditions $`m_z=\pm 1`$ at $`x=\pm \mathrm{}`$. By doing that we obtain a DW thickness, $`W_X`$, of the order of $`\mathrm{}`$. In order to get a functional to describe a DW, we write $`𝐦_{}=\mathrm{sin}2\theta (x)(\mathrm{cos}\varphi (y),\mathrm{sin}\varphi (y))`$, and integrate in Eq.(3) over $`x`$ using a simple model in which $`\theta (x)`$ varies linearly through the DW, obtaining:

$`\mathrm{\Delta }_{DW}(\varphi (y))=\rho /2{\displaystyle 𝑑y(_y\varphi (y))^2}`$
$`+B{\displaystyle 𝑑y[(_y\varphi (y))^4+(_y^2\varphi (y))^2]}+V_H+\mathrm{\Delta }_{SO}`$ (6)
where $`\mathrm{\Delta }_{SO}`$ is a SO term that will be essential in the discussion below. Using the simple model $`\theta (x)=\pi x/2W_X`$ for $`x<W_X`$ and zero otherwise, the parameters in (6) are $`\rho =W_X\rho _{}/2`$ and $`B=W_X\widehat{\rho }/2`$. However, the rapid change of $`m_z(x)`$ over a magnetic length, raises some doubts on the validity of the functional (3) as a good starting point to obtain (6). Therefore, we have taken the alternative of performing a microscopic Hartree-Fock (HF) calculationJungwirth  for describing DW’s. In Fig. 1 we plot the HF quasiparticle energies as a function of the orbit guiding center. The chemical potential is located at the gap energy. The reduction of the energy gap at the DW is an indication of the loss of coherence of the wave-function. We find that $`W_X`$ is roughly $`2\mathrm{}`$ and the energy per magnetic length of the DW is $`0.0448`$. In the inset of Fig. 1 we plot the $`z`$-component of the unitary vector field, $`m_z`$, isospin as a function of the position. At the center of the DW $`m_z=0`$ and $`m_{}`$ should be the unity. In absence of spin-orbit coupling, the system has $`U(1)`$ symmetry and the energy of the DW does not depend on global rotations of $`𝐦_{}`$. ¿From the HF results, we find that the functional (6) is adequate for describing textures of $`𝐦_{}`$ along the DW. The coefficients $`\rho `$ and $`B`$ for the terms with derivatives can be obtained from a fitting to the HF results. The ratio $`\rho /B`$ is the same than for the simple model above, but each coefficient has increased in a factor 3.7.

Since $`\rho <0`$, the first term in (6) tends to produce a rotation of the isospin along the DW. Although this rotation is limited by the second term, a fan phase could appear if one neglects any Hartree contribution as it is usually doneFalko ; Mitra . However, a rotation implies the existence of electrical dipoles associated to oscillations of the topological charge (but with zero total topological charge)

$`q(𝐫)={\displaystyle \frac{1}{4\pi }}_y\varphi (y)_x(m_z(x)).`$ (7)
The $`V_H`$ prevents the appearance of a fan phase induced by the negative stiffness. The Hartree repulsion of the charge density associated to the texture keeps the spin direction constant along the DW. It must be stressed that the Coulomb term in Eq.(6), has the same dependence in derivatives of the field $`\varphi `$ that the elastic term, and therefore it can not be neglected in the study of DW’s.

We have also included in Eq.(6) a SO term $`_{SO}`$. The SO interaction couples directly a state $`0,X,`$ with a state $`1,X,`$ Hanna  producing a Zeeman-like coupling to the isospin and an effective in plane magnetic field. Therefore, in our functional, SO is described by a term,

$$\mathrm{\Delta }_{SO}=\lambda _{SO}𝑑y\left(cos\varphi (y)1\right)$$
(8)
with $`\lambda _{SO}=W_X\beta _{SO}/2^{3/2}\pi ^2`$ where $`\beta _{SO}`$ is the bulk spin-orbit coupling Hanna .

Charged domain wall The solutions of Eq.(6) can be characterized by integers which correspond to the total topological charge $`Q_T`$ of the solution. $`Q_T`$ is the increase, in units of $`2\pi `$, of the phase $`\varphi (y)`$ when going from $`\mathrm{}`$ to $`+\mathrm{}`$. Hitherto, we have just considered solutions with $`Q_T=0`$. Let us now consider the solutions for $`Q_T>0`$.

Solutions of Eq.(6) in the sector $`Q_T=1`$, are very important since the equivalence between topological and electrical charge allows the isospin textures to be the relevant charged excitations in the systemSondhi ; Fertig . In the presence of domains, charge excitations can be trapped in the walls forming confined isospin textures, which are solitons in the phase $`\varphi (y)`$Falko . Analytical expression for the soliton have been obtained, neglecting the Hartree interaction, in the case of positive stiffnessFalko . In our case, the Hartree term is essential and we have not been able to obtain an analytic solution. For $`Q_T=1`$ we take a simplified shape for the soliton. In the sector of $`Q_T=1`$, we look for solitons of size $`\xi `$ having a simple form $`\varphi (y)=2\pi y/\xi `$, within an interval of length $`\xi `$ and zero out of that interval. The spin texture and the charge density of this soliton is shown schematically in Fig. 2. The size $`\xi `$ of the soliton is determined by the competition between the different terms of the functional; the Hartree and the quartic term ($`B>0`$) tend to make the texture large whereas the SO and the quadratic term ($`\rho <0`$) try to make it small. Fig. 3 shows the energy and size of the charged wall as a function of $`\lambda _{SO}`$. For $`\lambda _{SO}`$=0, the functional (6) has $`U(1)`$ symmetry and the soliton has zero energy being extended to the whole wall, i. e. $`\xi =L_y`$. SO interaction ($`\lambda _{SO}0`$) reduces the soliton size. $`\xi `$ takes a value much smaller than $`L_y`$ and the energy of the excitation becomes finite. The energy of the soliton is the energy cost to add an electron to the DW. Since this value is smaller than the energy gap in the S and F phases, $`e^2/ϵ\mathrm{}`$, we expect that extra charges in the system will become located at the DW. The coupling $`\lambda _{SO}`$ depends of the system characteristics as the DW width $`W_X`$. Typical values Hanna ; Falko2  vary from $`2\times 10^4`$ to $`8\times 10^4`$. In this range, the energy of the soliton is much smaller than the DW gap obtained in HF calculations (Fig. 1), $`e^2/ϵ\mathrm{}`$. The HF gap is dominated by exchange Coulomb interactions, and represents the excitation gap when the isospin order parameter is held fixed. The actual low-energy charge excitations come from fluctuations of the order parameter field. Once again, $`V_H`$ has been essential in the properties of the (in this case charged) DW.

Transport properties of the DW. Let us analyze the transport through a DW. If the chemical potential of the system, fixed by impurities or edge states, is located at the energy gap of the DW, no current can flow parallel to the DWMitra2 . The only possibility for the carriers is to pass across the DW. On the contrary, when the chemical potential resides in a band, there is a perfect unity transmission along the direction parallel to the wall and no carriers are passing through the DW. The same argument is valid for the charged excitations gap instead of that of the uncharged DW. The current through a DW separating a F from a U phase is different from zero if and only if the chemical potential lies on the charged excitation gap of the DW.

In the absence of SO coupling, there is not a gap for the charged excitations and, consequently, transport across the DW is not possible. This is in agreement with spin conservation arguments; when $`\lambda _{SO}=0`$, the spin is a good quantum number and no transport of charge through the DW is possible unless some other scattering mechanism is able to flip an electron spin. The hyperfine coupling to nuclear spins has been sometimes invokedSmet , but a non zero SO coupling is much more efficient to flip electron spins. Due to SO, the solution of the functional (6) changes smoothly its isospin when going from one side to the other of the wall. One electron with a given (real) spin can pass across the wall smoothly flipping its spin. The finite, due to SO coupling, energy of the soliton is rather small which means that very few electrons pass across the wall flipping their spins because a small gap reflects coupling between very few states at the two sides of the barrierBrey . In other words, there is a small current passing across a domain wall with a large resistance. This explains the large resistance observed in different systems where domains existPiazza ; Poortere ; Smet ; Eom ; Kronmuler ; Hashimoto .

A final question to comment on is the role played by nuclear spins. Apart from the possible role played in the process of domain formation, nuclear spins are not needed, in our picture, in the process of carrier transport. However, due to the hyperfine interaction, nuclear spins will suffer a dynamic nuclear spin polarization within the electronic domains. This is very important because, if current is turned off for a while, as done in some experimentsHashimoto ; Smet2 , the electrons in different domains immediately lose memory of theirs spins. However, nuclear spins relax so slowly in time that they serve as memory reservoirs of spin states and, if electronic current is reestablished after a while, the domains will reappear in exactly the same position they had before.

In summary, we study DW in a quantum Hall Ising ferromagnet at $`\nu =2`$ by means of a functional theory supported by a HF calculation. In a neutral DW, Hartree repulsion prevents the appearance of a fan phase provoked by a negative stiffness. When the system is charged, electrons are trapped as solitons at the DW. Hartree and SO interactions determine the energy and size of these solitons. Finally, a discussion of transport through the DW is presented.

We are indebted to L. Martin-Moreno for helpful discussions. Work supported in part by MEC of Spain under contract No. PB96-0085.


format: text

text:
# Achieving Extreme Resolution in Numerical Cosmology Using Adaptive Mesh Refinement: Resolving Primordial Star Formation11footnote 1 Further information and visualizations can be found at http://www.TomAbel.com/GB/

## 1 Introduction

Cosmic structure is formed by the gravitational amplification of initially small density fluctuations present in the early universe. A fluctuation containing the mass of the Milky Way galaxy will collapse by a factor of $`10^3`$ in size before it comes into dynamical equilibrium. In order to adequately resolve its internal structure, another two orders of magnitude of spatial resolution per dimension are needed, at a minimum, for a total spatial dynamic range $`SDR=10^5`$. Resolving the formation of individual stars in an entire galaxy would require vastly more resolution: $`SDR10^{20}`$—a seemingly unreachable goal.

N-body tree codes are widely used in numerical cosmology to achieve high spatial dynamic ranges in general 3-d evolutions, and some parallel implementations have won past Gordon Bell Awards . The largest cosmological N-body simulation used just over $`10^9`$ particles and achieved a $`SDR10^4`$ . The highest dynamic range N-body simulation achieved $`SDR=2\times 10^5`$ with substantially fewer particles . These calculations simulate only the collisionless cold dark matter (CDM) which dominates the gravitational dynamics of structure formation.

We are interested in simulating the formation of cosmic structures including the all–important baryonic gas which forms the visible galaxies. Our original design goal was to achieve $`SDR=10^4`$, a mark which we have far surpassed. In this paper we describe the Enzo cosmological adaptive mesh refinement (AMR) code we have developed for parallel computers and its application to a simulation of the formation of the first stars in the universe. The time-dependent calculation is carried out in full 3-d on a structured adaptive grid hierarchy which follows the collapsing protogalaxy and subsequent protostellar cloud to near stellar density starting from primordial fluctuations a few million years after the big bang. More than 8000 subgrids at 34 levels of refinement are generated automatically to achieve a final spatial and temporal dynamic range of $`10^{12}`$ (for comparison, $`10^{12}`$ is roughly the ratio of the diameter of the earth to the size of a human cell). It is the highest dynamic range 3-d simulation ever carried out in astrophysics.

The calculation incorporates all known relevant dynamical, chemical, and thermodynamic processes and is carried out in a proper expanding cosmological background spacetime. Enzo combines an Euler solver for the primordial gas , an N-body solver for the collisionless dark matter , a Poisson solver for the gravitational field , and a 12-species stiff reaction flow solver for the primordial gas chemistry . The latter is needed to determine the nonequilibrium abundance of molecular hydrogen which dominates the radiative cooling of the gas. This system of equations is solved on every level of the grid hierarchy with full two-way coupling, making it one of the most complex AMR simulations ever carried out. Adaptive mesh refinement allows us to resolve locally all the important length and timescales everywhere within the gravitationally collapsing cloud at all times, giving us confidence that our results are accurate. The calculation is effectively an ab initio simulation of star formation which connects initial conditions to the final state in full generality.

In this paper, we concentrate on the technical and performance aspects of the AMR code, but also discuss some of the exciting astrophysical results we have obtained. As we will see, achieving extreme resolution is more a matter of high performance data structures and extended precision arithmetic than raw gigaflops, although the latter is certainly needed.

## 2 Physical Model

The physics of present-day star formation is complicated and fraught with uncertainties. This is due to the relative complexity of the interstellar medium: it is composed primarily of hydrogen and helium but also contains heavier elements which are important for radiative cooling and heating. In addition, there are dynamically important magnetic fields, radiative backgrounds and uncertain initial conditions.

In contrast, the formation of the first star takes place in a much simpler environment: the gas composition is essentially just hydrogen and helium (since heavier elements are formed inside stars), there is no important radiative source, the magnetic fields are most likely insignificant and — most importantly — the initial conditions are precisely specified by cosmological models.

This makes primordial star formation a clean initial value problem in which we can include, from first principles, all of the relevant physical processes. Moreover, the problem is extremely important – it provides the starting point for the formation for all other structure in the universe, from galaxies to superclusters. The combination of all of these elements is rare in physics, and efforts to solve the problem date back over three decades. However, progress has been difficult because of the extreme multi-scale nature of the problem, as discussed in the introduction. In the following, we review the physical ingredients which are required to solve this problem, emphasizing the numerical implications.

### 2.1 Initial Conditions

In modern cosmological models, cold dark matter (CDM) dominates the dynamics of structure formation because its mass density exceeds that of the ordinary baryonic material by a factor of ten or more. Although the composition of CDM is unknown, its cosmologically relevant properties—its mean density and the functional form of its power spectrum of density fluctuations $`P(k)`$—are known or are calculable once a Friedmann “world” model is specified. We simulate a “standard” CDM model which means that the cosmological parameters and the amplitude of the power spectrum are those that have been shown to reproduce the statistical properties of galaxies and clusters of galaxies in the present universe .

The CDM power spectrum has the interesting property that rms density fluctuations are logarithmically divergent on small mass scales, implying that structures form hierarchically from the bottom-up: small gravitational condensations form first which then merge to form larger, more massive, objects later on. This merging behavior is captured numerically by combining an N-body technique, to follow the motion of the collisionless dark matter particles, with a Poisson solver, to compute the self-gravity of the dark matter (and the gas). But when does collapse to high density occur? How does it proceed? And what kind of astrophysical object is produced?

### 2.2 Gas Physics and Chemistry

In order to answer these questions, we need an accurate description of the thermal state of the gas. The fluid equations plus the ideal gas law (and Poisson’s equation) describe the evolution of the gas density, velocity and energy, but require knowledge of the radiative heating and cooling rates. This, in turn, depends on the chemical composition.

The primordial gas contains about 76% hydrogen and 24% of helium (trace amounts of deuterium and lithium play no important role here). For such a gas, and for temperatures below $`10^4`$ K, the primary cooling agent is molecular hydrogen, which has rotational and vibrational excited states than can be collisionally excited. Because the cosmological background density of baryons is small, chemical reactions in the smooth background gas occur on long timescales. As a consequence chemical equilibrium is rarely an appropriate assumption. The number of possible chemical reactions involving H2 is large, even in the simple case of metal free primordial gas. We have tediously selected the dominant reactions and collected the most accurate reaction rates available . We solve the time dependent chemical reaction network involving twelve species (including deuterium and helium). A fast numerical method to solve this set of stiff ordinary differential equations has been developed by some of us .

Once the chemical state is known, we must determine the radiative cooling rate. For densities encountered in the collapsing cloud, it is accurate to assume the gas is optically thin and only consider excitations from atomic and molecular ground states. We include all known radiative loss terms due to atoms, ions, and molecules that are appropriate for our primordial gas. Also the energy exchange between the cosmic microwave background and free electrons (Compton heating and cooling) is included.

## 3 The Enzo Cosmological Adaptive Mesh Refinement Code

To solve the star formation problem using a traditional static, uniform mesh is simply not feasible. Even assuming a continuation of Moore’s law in its present form, it would not be until about 2200 that a problem of this dynamic range could even fit into memory of the largest systems. The solution method that we adopt is structured adaptive mesh refinement (SAMR: ), because it: (1) is spatially- and time-adaptive, (2) uses accurate and well-tested grid-based methods for solving the hydrodynamics equations, and (3) can be well optimized and parallelized. In this section, we describe our implementation; more details are available elsewhere .

### 3.1 Structured Adaptive Mesh Refinement

The central idea behind SAMR is deceptively simple. While solving the desired set of equations on a uniform grid, monitor the quality of the solution; when necessary, add an additional, finer mesh over the region that requires enhanced resolution. This finer (child) mesh obtains its boundary conditions from the coarser (parent) grid or from other neighbouring (sibling) grids with the same mesh spacing. The finer grid is also used to improve the solution on its parent. As the evolution continues, it may be necessary to move, resize or remove the finer mesh. Even finer meshes may be required, producing a tree structure that may continue to any depth.

We show an example of this for the two-dimensional case in Figure 1. Each level up the hierarchy decreases the mesh spacing relative to its parent by a factor $`r`$ (the refinement factor). Therefore, if we define the resolution (or SDR) as the ratio of the length of the entire region to the cell spacing, the resolution at a level $`l`$ is given by $`nr^l`$, where $`n`$ is the number of cells in the root grid.

There are a number of things to note about this algorithm. First, in order to be efficient, the covering fraction of the refined meshes must be relatively small. For our problem the covering fraction of the most refined grids is extremely small, approximately $`10^{30}`$ (nevertheless, this is by far the most interesting region!). Second, the refinement factor is constrained to be an integer so that meshes can be aligned, making communication between grids easier and more accurate. In addition, the subgrids must be rectangular and completely contained within their parents. Third, the finer mesh is in addition to, not instead of, the coarse mesh it covers. At first this seems wasteful; however, the overhead is quite small and there are significant advantages. The individual grids remain uniform, so we can use off-the-shelf solvers. Also, since the timestep constraint depends on the mesh spacing, the uniform (coarse) mesh has a uniform (large) timestep, thus naturally permitting multiple timesteps at different levels.

### 3.2 Advancing the Grid Hierarchy

To advance our system of coupled equations in time as represented on this hierarchy requires some care. The overall control algorithm is quite similar to that used for a single grid. The heart of the algorithm is the recursive EvolveLevel routine which is passed the level it is to work on and the new time:

```
EvolveLevel(level, ParentTime)
begin
  SetBoundaryValues(all grids)
  while (Time < ParentTime)
  begin
    dt = ComputeTimeStep(all grids)
    SolveHydroEquations(allgrids, dt)
    Time += dt
    SetBoundaryValues(all grids)
    EvolveLevel(level+1, Time)
    FluxCorrection
    Projection
    RebuildHierarchy(level+1)
  end
end

```

The interior loop is used to advance the grids on this level from Time (the time at their current state) to ParentTime (the time to which the parent grids have been advanced). Inside this loop there is a recursive call so that all finer levels are advanced as well. The resulting order of timesteps is like the multigrid W cycle, as shown in Figure 2.

#### 3.2.1 Solving the equations of hydrodynamics

We treat each grid as its own self-contained unit, so that each routine in the pseudo-code above is applied to all the grids on a given level. As with any hyperbolic equation, we first need to set the boundary conditions on the grids. For the root grid, we use a predefined boundary condition (periodic in this case), while for finer grids, a two step procedure is used:

1. All boundary values are first interpolated from the grid’s parent.
2. Grids which border other grids on the same level (i.e. siblings) use the solution from the sibling grid.

This ensures that all boundary values are set using the highest resolution solution available.

Once the boundary values have been set, it’s straightforward to solve the equations of hydrodynamics on a grid. Because it is a Cartesian rectangular grid, we can use highly-optimized solvers. We have implemented two: the piecewise parabolic method (PPM)  modified for use in cosmological hydrodynamics , as well as a robust finite difference technique . This allows us a double check on any result.

A disadvantage to such a hierarchy is that we must ensure consistency in regions which are represented by both coarse and fine meshes. This involves two steps: the first is to correct the coarse fluxes (of conserved quantities) at subgrid boundaries to reflect the improved flux estimates from the subgrid. This is required to ensure mass, momentum and energy conservation as material flows into and out of a refined region. The second step, termed projection, updates the solution on the coarse mesh points which are covered by finer meshes. Taken together, these two steps represent one side of the two-way communication between parent and child grids (the other side being boundary interpolation from parent to grid).

#### 3.2.2 Rebuilding the Grid Hierarchy

The final part of the EvolveLevel routine adapts the grid hierarchy to the new solution. The RebuildHierarchy procedure modifies (“rebuilds”) the grids on the specified level and higher. This involves three steps:

1. A refinement test is applied to the parent grids of the current level, resulting in a boolean field indicating if this mesh point needs to be refined.
2. Rectangular regions are chosen which cover all of the refined regions, while attempting to minimize the number of unnecessarily refined points. This is done with an edge-detection algorithm from machine vision studies .
3. Finally, the new grids are created and their values are copied from the old grids on the same level where possible, or are otherwise interpolated from the parent grids. The old grids are then deleted, freeing their memory.

This process is repeated on the next refined level, until the entire hierarchy has been rebuilt. This is both easier and more effective than moving and resizing the grids, and still only takes about 10% of the overall cpu time.

#### 3.2.3 Refinement Criteria

In this application, there are three criteria used to determine when refinement is required. Since gravitational collapse causes mass to flow into a small number of cells, the first two are designed to preserve a given mass resolution in the solution. We define a threshold $`M_{}`$ such that whenever a cell accumulates at least this much mass, it is refined. This is done separately for the gas and the dark matter component.

The third refinement criterion ensures that all possible collapsing fragments are resolved. More specifically, we require that the cell width be less than some fraction of the local Jeans length ($`\mathrm{\Delta }x<L_J/N_J`$), the smallest perturbation which is unstable to gravitational collapse. We have varied $`N_J`$, the number of cells across the local Jeans length, from 4 to 64 without seeing a significant difference in the results.

### 3.3 Adding Physics: Gravity, Dark Matter and Chemistry

To simulate a system such as a collapsing star, we must also solve the elliptic Poisson equation on the hierarchy. On the root grid, this is done with an FFT which naturally provides the periodic boundary conditions required. On subgrids, we interpolate the gravitational potential field and then solve the Poisson equation using a traditional multi-grid relaxation technique. In order to produce a solution that is consistent across the boundaries of sibling grids, we use an iterative method: first solving each grid separately, exchanging boundary conditions, and then solving again.

In contrast to the familiar baryonic fluid, the dark matter is pressureless and collisionless, only interacting via gravity. Therefore, it is not appropriate to use the fluid equations; rather, we solve for the individual trajectories of a representative sample of particles. This we do using particle-mesh techniques specially tailored to adaptive mesh hierarchies. This adds another layer of complexity to the software, the details of which are discussed elsewhere .

Finally, as discussed earlier, the cooling of the primordial star is controlled by the formation of molecular hydrogen. Therefore, we need to follow the chemical evolution of various ionized states of hydrogen, helium, deuterium, and molecular hydrogen, twelve in all. This coupled set of ordinary differential equations can be solved independently for each grid point. Because the equations are stiff, we use a backward finite-difference technique for stability, sub-cycling within a fluid timestep for additional accuracy .

### 3.4 Parallelization Strategy

The code is implemented in C++, with compute-intensive kernels in FORTRAN 77. An object oriented approach provides a number of benefits. The first is encapsulation: a grid represents the basic building block of AMR. We have defined a series of atomic and binary operations that perform functions such as solving a certain set of equations on a grid or interpolating from one grid to another. This greatly simplifies the control algorithms. The second is extensibility: adding new equations is generally a matter of adding a few methods to a grid class and a few lines to the control algorithm.

Efficiently parallelizing SAMR is difficult, particularly for distributed memory systems. Grids have a relatively short life, so information must be updated frequently. Moreover, load balancing becomes a serious headache since small regions of the original grid eventually dominate the computational requirements.

As a first attempt, we tried compiler-based parallelization (a predecessor to OpenMP). Using pragmas, we identified large-grained loops that performed operations over entire grids. This was successful for a relatively small number of processors (8-16), but could not be extended to distributed memory systems, such as clusters of PCs. To adapt to these machines, we developed an MPI version that exploited the following optimization techniques:

* Distributed objects. We leveraged the object-oriented design by distributing the objects over the processors, rather than attempting to distribute an individual grid. This makes sense because the grids are generally small ($`20^3`$) and numerous (sometimes in excess of 50,000 grids). This had the added benefit of ease of coding: we simply added a small number of communication primitives to the class methods.
* Sterile objects. Although distributing the objects results in good load balancing, it can greatly increase the amount of communication since each processor has to probe other processors to find potential neighbours. We solved this problem by creating a type of object which contained information about the location and size of a grid, but did not contain the actual solution. These sterile objects are small and so each processor can hold the entire hierarchy. Only those grids which are local to that processor are non-sterile. This means that almost all messages are direct data sends; very few probes are required.
* Pipeline communications. One result of distribution is that all binary operations (e.g. obtaining boundary values) are potentially non-local. However, a given stage of the calculation generally involves many sends. We optimize this by dividing each stage into two steps. First, all of the data (such as boundary values) are processed and sent. Since all processors have the location of all other grids locally (thanks to the sterile objects), we can order these sends such that the data that are required first are sent first. Then, in the receive stage, the data needed immediately have had a chance to propagate across the network while the rest of the sends were initiated. This was implemented through the use of MPI’s asynchronous communication mode, and resulted in a large decrease in wait times.

### 3.5 Extended Precision Arithmetic

Performing a simulation with such extreme dynamic range presents unique difficulties. One of these is simply the requirement of being able to accurately describe the positions of grids and particles within the problem domain. At a minimum, we must be able to distinguish between the two locations $`x+\mathrm{\Delta }x`$ and $`x`$, where $`\mathrm{\Delta }x`$ is the minimum grid spacing. This implies a precision of at least $`\mathrm{\Delta }x/x10^{12}`$ (in fact we require roughly 100 times larger than this because of various mathematical operations which are applied to this ratio). Therefore, we need 128 bit precision. Unfortunately, this is not always available and when available, we often suffer a large performance penalty. For example, the GNU FORTRAN compilers does not support 128 bit precision. The SGI compiler on the Origin2000 does provide 128 bit arithmetic, but is some 30 times slower than 64 bit. The RS/6000 compiler on the SP2 allows 128 bit precision, but only as an extra compiler option, so libraries (such as MPI) must be re-compiled.

One option for dealing with this issue is to develop routines which use native 64 bit arithmetic to perform 128 bit operations, circumventing patchy native support. Recent work  in this direction appears to be quite successful and we are examining the possibility of adopting this method.

Whichever approach is taken, it is clear that the use of 128 bit arithmetic should be kept to a minimum because of its impact on operation count, memory consumption and interprocessor message size. Throughout the code, we have identified only those operations which require high precision. These correspond to operations involving absolute position and time, as opposed to relative location (O($`\mathrm{\Delta }x`$)). Most grid atomic operations can be entirely written in terms of relative position and time, and hence low precision. This reduced the total high-precision operation count to $`5\%`$ of the total, resulting in considerable speed (and memory) improvements.

## 4 Astrophysical results

In the calculation discussed in this paper, we have set up a three dimensional volume that is 256 comoving kiloparsec (kpc) on a side (1 parsec is 3.26 light years, roughly the distance to the Sun’s nearest neighbour). We set up the initial conditions a few million years after the big bang, using an inflation-inspired cosmological model. We first run a low-resolution ($`64^3`$) simulation to determine where the first star will form and then restart the calculation including three additional levels of static meshes, covering the region from which the star will form (equivalent to $`512^3`$ initial conditions over the entire box). This allows us to increase mass resolution in this region by a large factor (512) and capture as many small-wavelength modes in the initial conditions as possible. We have experimented with using only two additional levels and find it has little effect on the overall evolution of the object. We have also carried out a number of experiments varying the refinement criteria and find the results described here are quite robust.

A series of earlier AMR simulations at increasing resolutions have revealed the following results. A simulation with $`SDR=2.5\times 10^5`$  showed that at redshifts of $`20`$ (approximately 150 million years after the big bang) primordial gas becomes gravitationally bound to dark matter condensations with characteristic mass of few $`10^5`$ M (solar masses). Non-equilibrium chemical reactions form $`H_2`$ from primordial hydrogen in the cores of these clumps in sufficient concentrations to cool the gas down to a few hundred K via molecular line radiation. A simulation with $`SDR=3\times 10^7`$  showed the onset of gravitational instability in this cool gas and the formation of a collapsing protostellar core. The mass of the collapsing cloud was a few hundred solar masses. Most recently, a simulation with $`SDR=10^{10}`$  resolved the supersonic accretion flow onto a central, fully molecular protostellar core, confirming the resulting star would be much more massive than the Sun.

The simulation reported here attempts to follow the protostar’s collapse to stellar density, which we estimate will require $`SDR=10^{15}`$. As of current writing, we have achieved $`SDR=10^{12}`$ We find a rapidly accreting protostar that will form a single massive star with a mass of order one hundred times that of our own sun. Such a large star will explode as a supernovae, generating the very first heavy elements in the universe. In Figure 3, we show a zoom into the forming star. At this stage, the protostar is still collapsing and not yet spherical since gravity still overwhelms internal pressure. We highly encourage readers to view the animated version of this zoom at www.TomAbel.com/GB/movies.html.

Although the cloud and protostar are not spherical, it is instructive to plot in Figure 4 radial profiles of mass-weighted spherical averages of various quantities for seven illustrative times of the evolution. At the first output time shown (marked z=19), about 700 $`M_{}`$ have cooled below 200 Kelvin at the center of the protogalactic halo of mass $`5.44\times 10^5M_{}`$. From the radial velocity profiles it is evident that cooling gas behind the accretion shock (at $`r100`$ pc) falls onto this cooled central material. This cold central region within $`r10`$ pc is the high redshift analog of a present day molecular cloud despite its minute molecular mass fraction of $`10^3`$ (see panel C of Figure 4).

As the evolution continues (the second output time shown in Figure 4) the total mass of the proto-galactic object has grown to $`7.4\times 10^5M_{}`$ at $`z=18.209`$. The primordial “molecular cloud” contains now $`1000M_{}`$. At the same time a central core of $`100M_{}`$ starts to become warmer ($`300`$ K) within the “molecular cloud” and shows small radial velocities. For this central core the radial velocities and temperatures are seen to increase monotonically as time goes on. The molecular mass fraction also starts to increase slightly. When the central densities exceed $`10^9`$ cm-3 three–body H2 formation shows its effect in binding atomic hydrogen to molecules. At central densities $`10^{11}`$ cm-3 atomic and molecular hydrogen exist in similar abundance. This enhanced molecular fraction which decreases the cooling time is the reason for the observed central temperature drop (see panel D of Figure 4). Interestingly, the inner $`1M_{}`$ shows similar thermal evolution and start to collapse more rapidly than its surrounding. In fact as can be seen in panel E of Figure 4 it eventually starts to contract faster than the local sound speed. Before this dynamical collapse all the cooling gas has been settling to the center with sound crossing times shorter than or equal to the dynamical times, and long cooling times. I.e. the contraction proceeds quasi–hydrostatically until 3 body H2 formation turns the gas fully molecular. During this slow contraction phase sound waves damp sub–Jeans density perturbations in the baryons yielding a rather smooth gas distribution prior to the onset of 3 body H2 formation. The time-scales become as short as weeks for the densest material followed in the simulations. We have not observed any signs of fragmentation by examining the three–dimensional density field of the final output time. Although the protostar has not yet reached stellar density, we conclude that a single massive star will ultimately form.

## 5 Code Performance Results

Our primary interest is not raw performance, but rather the complexity and difficulty of the problem and the algorithm. Therefore, we first examine some characteristics of the resulting grid hierarchy.

In the top two panels of Figure 5, we show how the grid hierarchy grows as time progresses. Note the slow increase in the number of grids as the proto-star condenses and the final, very sudden jump in the depth of the grid tree at the end when the core of the cloud collapses to high density. This demonstrates how the data structures themselves adapt to fit the physical solution. Note also the extremely large number of memory allocations and frees, since the entire grid hierarchy is rebuilt thousands of times. This kind of method represents a new class of scientific computing that place great strain on the operating system infrastructure. Total memory usage is also substantial, often reaching up to 20 GB. With outputs in the 2-4 GB range, we require at least 50-100 GB disk storage and much more mass storage space.

In the bottom two panels of this figure, we have chosen two representative times and plotted the distribution of grids per level. At early times, most of the grids are at moderate levels, representing the fact that relatively low resolution is sufficient to model the proto-stellar cloud. However, at late times, as the proto-star forms, a large investment is required at the very highest levels of resolution.

Finally, we turn to the raw performance of the code. Because of the adaptive nature of the algorithm it is difficult to determine an over-all flop rate for the entire computation. One solution to this problem would be to instrument each module to return its operation count (which depends on such factors as the grid size, particle number, and other parameters, as well as the solution itself in the case of some modules). Then we could simply add up the count returned by all the modules. However, the code is nearly 100,000 lines, so this remains a future project.

In the meantime, we have estimated the flop rate in the following way. We have used the hardware floating-point counter on the SGI Origin2000 to determine the total floating point operation count for a portion of the simulation (drawn from near the middle of the run, which should be representative). In doing so, we have been careful to use only 64 bit precision, rather than 128 bit (this somewhat affects the accuracy of the result but because it is only used for timing this is not a concern). Then, when we run the calculation on our production machine, the Blue Horizon SP2 system at the San Diego Supercomputer Center, we use full 128 bit precision and measure the wall clock time for this same portion of the simulation. This ensures that the operation count employs a widely-used standard (the hardware counter on the R10000 processor), while the timing reflects the use of 128 bit precision. This provides a conservative estimate of the flop rate since we treat 128 bit operations as a single operation, even though the hardware itself may not do so. Given this operation count coupled with the measured wall clock time running on 64 processors of the Blue Horizon SP2, produced a total speed of approximately 13 Gflop/s. Since most of the run was done on 64 processors and the section timed is representative of the entire calculation, this represents our estimate of the sustained performance.

As an exercise, we can also ask how long this calculation would have taken with a traditional static grid code and compute an effective or virtual flop rate. To do this, we assume a grid with $`10^{12}`$ cells on each side, and assume the entire calculation would have taken (quite conservatively) $`10^{10}`$ timesteps. This works out to approximately $`10^{50}`$ floating point operations. Since the entire calculation took of order $`10^6`$ seconds, this converts to a virtual flop rate of $`10^{44}`$ flop/s.

We stress that this result has been achieved by optimizing all stages of the calculation. Since we solve a number of different kinds of equations with various techniques, it is instructive to examine the relative time taken by each part. In the table below, we examine the fractional part of the compute time taken by each of the following science components: (a) the hydrodynamics solver, (b) the gravity solver, (c) the chemistry solver, (d) the particle solver for the dark matter, as well as a few major parts of the amr structure. We stress that these are fractions of the compute time which, for our 64 processor run, represents approximately 60% of the total wall clock time, the other 40% being consumed by communication and imperfect load balancing. More detailed analysis of this code (on other problems) has been performed elsewhere .

| component | usage |
| --- | --- |
| hydrodynamics | 36 % |
| Poisson solver | 17 % |
| chemistry & cooling | 11 % |
| N-body | 1 % |
| hierarchy rebuild | 9 % |
| boundary conditions | 15 % |
| other overhead | 11 % |

## 6 Analysis and Visualization Challenges

Beyond simply running the simulation, we have also faced major challenges in data analysis and visualization of these large dynamic range simulations. Our analysis routines have had to be implemented in parallel to hold the gigabytes of data in memory on available systems. These routines facilitate finding collapsed objects and other regions of interest. They range from computing direct hydrodynamical quantities, such as temperatures and densities, to derived quantities like cooling times, two–body relaxation times, X–ray luminosities and inertial tensors. To study flattened objects such as galactic or proto stellar disks versatile routines to find such objects and derive projections, surface densities and other useful diagnostic quantities were created or modified so that they understood the structure of the hierarchy.

Visualization has proven equally involved. To allow interactive exploration of the full data sets on available systems we developed Jacques, a GUI-based visualization tool which allows simultaneous interactive analysis of tens of thousands of grids of the AMR hierarchy on modest memory machines. Jacques is implemented in the Interactive Data Language (IDL), and was used to produce the images in Fig. 3. The hundred solar radii protostar in our computational volume is tiny compared to a needle in a haystack and navigation techniques had to be devised to simplify the identification of regions of interest and to visualize slices and projections of them (Jacques has a “zoom in by $`10^{10}`$ button”!). Other features include velocity fields, isosurfaces, and a preliminary volume renderer.

## 7 Conclusions

We have presented results on the highest-resolution 3-d computational fluid dynamics simulation yet performed, achieving a spatial dynamic range of $`10^{12}`$ in the region of interest. The simulation breaks new ground in the complexity of the method and the realism of the result. On the computational front, we solve a series of coupled partial and ordinary differential equations with a range of techniques, in parallel, on a hierarchy of 34 levels containing nearly 10,000 grids. From an astrophysical point of view, we have performed a long sought feat: the ab initio simulation of one of the first stars in the universe. Put together, we argue that this represents a watershed in scientific computing.

We acknowledge support from NSF grant AST-9803137. Support for GLB was provided by NASA through Hubble Fellowship grant HF-01104.01-98A from the Space Telescope Science Institute, operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS6-26555. We gratefully acknowledge the support of the San Diego Supercomputing Center.


text:
# 1 Introduction

## 1 Introduction

The inverse-square interacting particle system  and its spin generalizations  are important models of many-body systems, due to their exact solvability and intimate connection to spin chain systems , 2-dimensional Yang-Mills theories  etc. Reference  is a classic review, while  present newer different perspectives.

The prototype of these models is the spin-Calogero (‘rational’) scattering model of particles on the line carrying $`U(\mathrm{n})`$ spin and interacting with two-body inverse-square potentials with a $`U(\mathrm{n})`$-invariant spin coupling. Most other models can be obtained as appropriate reductions of this model, taking advantage of its discrete or continuous symmetries . In particular, generalizations involving $`U(\mathrm{n})`$ non-invariant interactions can be obtained this way, recovering the trigonometric models derived in  and extending them to the elliptic case .

An unrelated development has been the recent progress in noncommutative field theory and matrix models. Spatial noncommutativity can be traced back to Heisenberg and naturally arises in lowest Landau level physics . Its current manifestation originates in matrix, string and membrane theory  and came into focus with the work of Connes, Douglas and Schwartz .

So far these two fields remained unrelated. In this letter we show how they can be cross-fertilized by borrowing notions of noncommutative geometry and applying them in the reduction scheme of the Calogero model to obtain a new integrable elliptic model involving non-$`U(\mathrm{n})`$ invariant noncommutative spin interactions. Such a modification of the spin interaction may serve to test the ‘flavor stiffness’ of the original spin model, to stress the degeneracy structure of the energy spectrum and to identify universality features of this class of models.

## 2 The reduction scheme

The basic technique that we will use consists of reducing a system of infinitely many particles with spin to a finite system with generalized interactions. The reader should refer to , and especially to the elliptic case with spin, for a more detailed description of the method.

The starting point is the spin-Calogero system with classical $`U(\mathrm{n})`$ degrees of freedom. This system can be obtained, for instance, from the model in  (which can itself be obtained as a reduction of a hermitian matrix model  into nontrivial angular momentum sectors) by redistributing the global $`U(N)`$ degrees of freedom of this model into individual particle spins, or, alternatively, from the infinite-volume classical limit of the spin model derived and solved in . The hamiltonian for $`N`$ particles reads

$$H=\underset{i=1}{\overset{N}{}}\frac{1}{2}p_i^2+\frac{1}{2}\underset{ij}{}\frac{\mathrm{tr}(S_iS_j)}{x_{ij}^2}$$
(1)
$`x_i`$ and $`p_i`$ are one-dimensional canonical coordinates and momenta; $`S_i`$ are a set of independent classical $`U(\mathrm{n})`$ spins of rank one and length $`\mathrm{}`$, that is, $`\mathrm{n}\times \mathrm{n}`$ rank-one hermitian matrices satisfying

$$\mathrm{tr}(S_i)^2=\mathrm{}^2$$
(2)
and with Poisson brackets

$$\{(S_i)_{ab},(S_j)_{cd}\}=i\delta _{ij}\left[(S_i)_{ad}\delta _{cb}\delta _{ad}(S_i)_{cb}\right]$$
(3)
Such spins can be realized in terms of oscillators :

$$(S_i)_{ab}=\overline{A}_i^aA_i^b,a,b=0\mathrm{}\mathrm{n}1$$
(4)
where $`(A_i^a,\overline{A}_i^a)`$ are a set of $`\mathrm{n}N`$ independent classical harmonic oscillator canonical pairs with Poisson brackets:

$$\{A_i^a,\overline{A}_j^b\}=i\delta _{ij}\delta _{ab}$$
(5)
and satisfying the constraint

$$\underset{a}{}\overline{A}_i^aA_i^a=\mathrm{}\mathrm{for}\mathrm{all}i$$
(6)
In the above model we can analytically continue the coordinates $`x_i`$ and momenta $`p_i=\dot{x}_i`$ to the complex plane. The integrability and solvability of this model trivially extends to the complex case. Such extensions will be useful, provided that we can identify a real subsystem, which will be the physical system of interest.

The hamiltonian $`H`$ is invariant under particle permutations, coordinate space translations and global spin rotations:

$$ip__N(i);x_ix_i+\mathrm{c};S_iUS_iU^1$$
(7)
where $`p__N`$ is an element of the permutation group $`S_N`$, $`\mathrm{c}`$ is a constant complex parameter and $`U`$ is a constant $`U(\mathrm{n})`$ matrix. We can arrange, therefore, for special initial configurations that are invariant under some combination of the above symmetries and be sure that these conditions will be preserved in time.

We shall choose the configuration to be a replication of $`N`$ real coordinates over an infinite complex lattice with periods $`\mathrm{c}_1`$ and $`\mathrm{c}_2`$. The $`N`$ particles effectively live on the real coordinate of a complex torus. The total number of particles on the covering space is, thus, infinite and we can parametrize them with the triplet of indices $`(i;m,n)`$, where $`i=1,\mathrm{}N`$ labels the particles in each cell and $`m,n𝐙`$ parametrize the cell. Shifts in $`m`$ and $`n`$ are elements of the (infinite) permutation symmetry of the system. The kinematical variables are chosen to obey

$`x_{i;m+1,n}=x_{i;m,n}+\mathrm{c}_1,p_{i;m+1,n}=p_{i;m,n}`$ (8)
$`x_{i;m,n+1}=x_{i;m,n}+\mathrm{c}_2,p_{i;m,n+1}=p_{i;m,n}`$ (9)
ensuring that we are dealing with lattice copies. This means

$$x_{i;m,n}=x_i+m\mathrm{c}_1+n\mathrm{c}_2,p_{i;m,n}=p_i$$
(10)
The above conditions, being an invariance under combined permutations and translations, are dynamically preserved. To also preserve the condition that $`x_i`$ are real, we should impose invariance under the imaginary parity transformation $`xx^{}`$. Choosing $`\mathrm{c}_1`$ real, the only possibilities for $`\mathrm{c}_2`$ are

$$\mathrm{c}_2+\mathrm{c}_2^{}=0\mathrm{or}\mathrm{c}_2+\mathrm{c}_2^{}=\mathrm{c}_1$$
(11)
The first choice ($`\mathrm{c}_2`$ imaginary) leads to an orthogonal lattice, while the second choice leads to a rhombic lattice.

To ensure full preservation of the lattice structure, we should also impose appropriate periodic conditions for the spins. In this, we can take advantage of the global spin symmetry of the system and impose

$`S_{i;m+1,n}`$ $`=`$ $`US_{i;m,n}U^1`$
$`S_{i;m,n+1}`$ $`=`$ $`VS_{i;m,n}V^1`$ (12)
with $`U`$ and $`V`$ two constant matrices. That is, spins can pick up $`U(\mathrm{n})`$ transformations as they move around the cycles $`\mathrm{c}_1`$ or $`\mathrm{c}_2`$. Consistency requires that $`S_{i;m+1,n+1}`$ be uniquely determined irrespective of the order of increase of the indices $`m`$ and $`n`$. That is,

$$UVS_{i;m,n}V^1U^1=VUS_{i;m,n}U^1V^1$$
(13)
which implies

$$[U^1V^1UV,S_{i;m,n}]=0$$
(14)
For this to hold for all generic $`S_i`$ we must require $`U^1V^1UV\omega `$ to be proportional to the identity matrix. Clearly $`\omega `$ satisfies $`det(\omega )=\omega ^\mathrm{n}=1`$, so we obtain

$$UV=\omega VU,\omega =e^{i2\pi \frac{\nu }{\mathrm{n}}}$$
(15)
with $`\nu `$ an integer $`0\nu <\mathrm{n}`$. $`U`$ and $`V`$ then satisfy Weyl’s braiding condition which characterizes a noncommutative (‘quantum’) torus .

The spin matrices $`S_{i;m,n}`$ are now expressed as

$$S_{i;m,n}=U^mV^nS_iV^nU^m=V^nU^mS_iU^mV^n$$
(16)
Inserting the expressions (10,16) in the hamiltonian (1) we obtain the hamiltonian of the reduced system. As usual, the resulting hamiltonian has an infinite factor, due to the summation of the hamiltonians of the infinitely many identical cells over the complex plane. Dropping this trivial infinity, the reduced hamiltonian includes the kinetic terms of the fundamental cell and the interaction potential of particles in this cell with all other particles in all cells:

$$H=\underset{i=1}{\overset{N}{}}\frac{1}{2}p_i^2+\frac{1}{2}\underset{i,j}{}\underset{m,n=\mathrm{}}{\overset{\mathrm{}}{}}\frac{\mathrm{tr}(U^mV^nS_iV^nU^mS_j)}{(x_{ij}+m\mathrm{c}_1+n\mathrm{c}_2)^2}$$
(17)
where we adopted the notation $`x_{ij}=x_ix_j`$.

## 3 Noncommutative spin interaction potentials

To proceed, we must identify the possible forms of $`U,V`$. We need the irreducible representations of the relation (15). Call $`k`$ the greatest common divisor of $`\nu `$ and $`\mathrm{n}`$. Then $`\mathrm{n}=k\mathrm{m}`$ and $`\nu =k\mu `$, for relatively prime $`\mathrm{m},\mu `$. The irreducible representations for $`U,V`$ are $`\mathrm{m}`$-dimensional ‘clock’ and ‘shift’ matrices. By a global $`U(\mathrm{n})`$ spin transformation we can diagonalize either of $`U,V`$. Choosing $`U`$ diagonal, the general form of $`U`$ and $`V`$ will be the direct sum of $`k`$ of the above irreducible representations:

$$U=diag\{e^{i\varphi _0},\mathrm{}e^{i\varphi _{k1}}\}u,V=diag\{e^{i\theta _0},\mathrm{}e^{i\theta _{k1}}\}v$$
(18)
where $`\varphi _q,\theta _q`$ are arbitrary phases, determining the Casimirs $`U^\mathrm{m}`$ and $`V^\mathrm{m}`$, and $`u,v`$ are the $`\mathrm{m}`$-dimensional clock and shift matrices

$$u_{\alpha \beta }=\omega ^\alpha \delta _{\alpha \beta },v_{\alpha \beta }=\delta _{\alpha +1,\beta }(mod\mathrm{m}),\alpha ,\beta =0,\mathrm{}\mathrm{m}1$$
(19)
So the acceptable $`U`$ and $`V`$ depend on $`2k`$ arbitrary parameters.

To take advantage of the form (18) for $`U,V`$ we partition $`S_i`$ into $`k^2`$ blocks of dimension $`\mathrm{m}\times \mathrm{m}`$ each by using the double index notation

$$(S_i)_{ab}=(S_i)_{\alpha \beta }^{pq},a=p\mathrm{m}+\alpha ,b=q\mathrm{m}+\beta $$
(20)
The $`U(\mathrm{n})`$ Poisson brackets in this notation are

$$\{(S_i)_{\alpha \beta }^{pq},(S_j)_{\gamma \delta }^{rs}\}=i\delta _{ij}\left[(S_i)_{\alpha \delta }^{ps}\delta _{\gamma \beta }\delta _{rq}\delta _{\alpha \delta }\delta _{rq}(S_i)_{\gamma \beta }^{rq}\right]$$
(21)
The $`m,n`$-sums that appear in (17) then become

$$\underset{m,n;\alpha ,\beta ;p,q}{}(S_i)_{\alpha +n,\beta +n}^{pq}(S_j)_{\beta ,\alpha }^{qp}\frac{e^{im\varphi _{pq}in\theta _{pq}}\omega ^{m(\alpha \beta )}}{(x_{ij}+m\mathrm{c}_1+n\mathrm{c}_2)^2}$$
(22)
where the term $`m=n=0`$ is omitted if $`i=j`$.

The above gives a potential interaction between particles $`i`$ and $`j`$ in the form of a modular function in $`x_{ij}`$ which depends on the spin components of particles $`i`$ and $`j`$. To make the noncommutative character of the spin interaction manifest, we perform a change of basis in the spin states and define

$$(\stackrel{~}{S}_i)_{\alpha \beta }^{pq}=\underset{\sigma }{}\omega ^{\left(\sigma +\frac{\alpha }{2}\right)\beta }(S_i)_{\alpha +\sigma ,\sigma }^{pq}$$
(23)
This is essentially a discrete Fourier transform in the sum of the $`\alpha ,\beta `$ indices of $`S_{\alpha \beta }^{pq}`$. (Note that, for $`\mathrm{m}`$ odd, $`\stackrel{~}{S}_{\alpha \beta }^{pq}`$ is actually antiperiodic in the index $`\alpha `$ if $`\beta `$ is odd, and vice versa. Although we could have defined a properly periodic matrix, we prefer this slight inconvenience in order to make the ensuing formulae more symmetric.) In fact, it will be convenient to assemble the double indices $`(\alpha ,\beta )`$ and $`(m,n)`$ into vectors $`\stackrel{}{\alpha }`$ and $`\stackrel{}{m}`$. Similarly, we define $`\stackrel{}{\mathrm{c}}=(\mathrm{c}_1,\mathrm{c}_2)`$ and $`\stackrel{}{\varphi }_p=(\varphi _p,\theta _p)`$.

The Poisson brackets of the $`\stackrel{~}{S}_i`$ are found from (21)

$$\{(\stackrel{~}{S}_i)_\stackrel{}{\alpha }^{pq},(\stackrel{~}{S}_j)_\stackrel{}{\beta }^{rs}\}=i\delta _{ij}\left[\omega ^{\frac{\stackrel{}{\alpha }\times \stackrel{}{\beta }}{2}}\delta _{ps}(\stackrel{~}{S}_i)_{\stackrel{}{\alpha }+\stackrel{}{\beta }}^{rq}\omega ^{\frac{\stackrel{}{\alpha }\times \stackrel{}{\beta }}{2}}\delta _{rq}(\stackrel{~}{S}_i)_{\stackrel{}{\alpha }+\stackrel{}{\beta }}^{ps}\right]$$
(24)
This is a structure extending the Moyal (star-commutator) algebra, the exponent of $`\omega `$ being the cross product of the discrete ‘momenta’ $`\stackrel{}{\alpha }`$ and $`\stackrel{}{\beta }`$. For $`(rs)=(pq)`$, in particular, it becomes the torus Fourier transform of the Moyal bracket

$$\{(\stackrel{~}{S}_i)_\stackrel{}{\alpha }^{pq},(\stackrel{~}{S}_j)_\stackrel{}{\beta }^{pq}\}=i\delta _{ij}(\omega ^{\frac{1}{2}}\omega ^{\frac{1}{2}})\left[\stackrel{}{\alpha }\times \stackrel{}{\beta }\right]_\omega (\stackrel{~}{S}_i)_{\stackrel{}{\alpha }+\stackrel{}{\beta }}^{pq}$$
(25)
where

$$[x]_\omega =\frac{\omega ^{\frac{x}{2}}\omega ^{\frac{x}{2}}}{\omega ^{\frac{1}{2}}\omega ^{\frac{1}{2}}}$$
(26)
is the $`\omega `$-deformation of $`x`$. This is the so-called trigonometric algebra with periodic discrete indices .

Finally, by inverting (23) and substituting in (22), the potential energy $`W`$ in terms of the $`\stackrel{~}{S}_i`$ acquires the form

$$W=\underset{i,j}{}\underset{\stackrel{}{\alpha };p,q}{}(\stackrel{~}{S}_i)_\stackrel{}{\alpha }^{pq}(\stackrel{~}{S}_j)_\stackrel{}{\alpha }^{qp}W_\stackrel{}{\alpha }^{pq}(x_{ij})$$
(27)
The above includes two-body interactions, for $`ij`$, as well as spin self-couplings, for $`i=j`$, arising from the interaction of each particle with its own images in different cells. The two-body potential $`W_\stackrel{}{\alpha }^{pq}(x)`$ is

$$W_\stackrel{}{\alpha }^{pq}(x)=\frac{1}{\mathrm{m}}\underset{\stackrel{}{m}}{}\frac{\omega ^{\stackrel{}{\alpha }\times \stackrel{}{m}}e^{i\stackrel{}{\varphi }_{pq}\stackrel{}{m}}}{(x+\stackrel{}{\mathrm{c}}\stackrel{}{m})^2}$$
(28)
while the spin self-coupling $`\stackrel{~}{W}_\stackrel{}{\alpha }^{pq}`$ is

$$\stackrel{~}{W}_\stackrel{}{\alpha }^{pq}=\frac{1}{\mathrm{m}}\underset{\stackrel{}{m}\stackrel{}{0}}{}\frac{\omega ^{\stackrel{}{\alpha }\times \stackrel{}{m}}e^{i\stackrel{}{\varphi }_{pq}\stackrel{}{m}}}{(\stackrel{}{\mathrm{c}}\stackrel{}{m})^2}$$
(29)
If the above potentials were independent of the $`U(\mathrm{n})`$ indices $`\stackrel{}{\alpha }`$ and $`p,q`$, the sum over $`U(\mathrm{n})`$ indices in the potential energy expression (27) would simply be a $`U(\mathrm{n})`$ trace and would give the $`U(\mathrm{n})`$-invariant coupling between the spins of particles $`i`$ and $`j`$ multiplying the standard Weierstrass potential of the elliptic Calogero model. In the present case, however, the above potential is spin-dependent and breaks $`U(\mathrm{n})`$ invariance, introducing a star-product twist in the indices $`\stackrel{}{\alpha }`$ and phase shifts $`\stackrel{}{\varphi }_p`$ in the indices $`p,q`$. Generically, the $`U(\mathrm{n})`$ invariance of the original model is broken down to an abelian $`U(1)^k`$, amounting to the transformation

$$(S_i)_{\alpha \beta }^{pq}e^{i\phi _p}(S_i)_{\alpha \beta }^{pq}e^{i\phi _q}$$
(30)
If $`\stackrel{}{\varphi }_p`$ are equal for $`k^{}`$ values of $`p`$, the remaining symmetry $`U(1)^k^{}`$ is enhanced to $`U(k^{})`$, corresponding to mixing the corresponding $`p`$-components.

The case $`\omega =1`$, $`\stackrel{}{\varphi }_p=0`$ reduces to the standard spin-elliptic Calogero-Moser model. The case $`\mathrm{m}=1`$ (and thus $`\omega =1`$) reproduces the $`U(\mathrm{n})`$-noninvariant model introduced in . The general case with $`\omega 1`$ is a new classical integrable model of the spin-Calogero type with a spin-dependent potential which is a modular function of the two-body distance $`x_{ij}`$.

The sums appearing in (28) and (29) could in principle have ambiguities due to the logarithmic divergence of the summation over the radial coordinate on the complex plane. This is, indeed, the case for the standard Weierstrass function and a specific prescription is needed to regularize it. Different prescriptions lead to different additive constants in the result. In our case, however, the presence of the extra phases renders the sums convergent and there is no regularization ambiguity.

The potentials can be expressed in terms of theta-functions. $`W_\stackrel{}{\alpha }^{pq}(x)`$ is a modular function on the complex torus $`(\mathrm{c}_1,\mathrm{c}_2)`$ with quasiperiodicity

$`W_\stackrel{}{\alpha }^{pq}(x+\mathrm{c}_1)`$ $`=`$ $`e^{i\varphi _{pq}+\frac{2\pi \mu }{m}\alpha _2}W_\stackrel{}{\alpha }^{pq}(x)`$
$`W_\stackrel{}{\alpha }^{pq}(x+\mathrm{c}_2)`$ $`=`$ $`e^{i\theta _{pq}\frac{2\pi \mu }{m}\alpha _1}W_\stackrel{}{\alpha }^{pq}(x)`$ (31)
It has a double pole at $`x=0`$, with principal part

$$W_\stackrel{}{\alpha }^{pq}(x)=\frac{1}{\mathrm{m}x^2}+O(x^0)$$
(32)
and no other poles in each cell. These properties uniquely define $`W_\stackrel{}{\alpha }^{pq}(x)`$ and allow for an expression in terms of theta-functions. We put

$$W_\stackrel{}{\alpha }^{pq}(x)=A\omega ^{i\frac{x}{\mathrm{c}_1}}e^{i\frac{x}{\mathrm{c}_1}\varphi _{pq}}\frac{\vartheta _1\left(\frac{\pi }{\mathrm{c}_1}(xQ_1)\right)\vartheta _1\left(\frac{\pi }{\mathrm{c}_1}(xQ_2)\right)}{\vartheta _1\left(\frac{\pi }{\mathrm{c}_1}x\right)^2}$$
(33)
where $`Q_{1,2}`$ are the as yet unknown zeros of $`W_\stackrel{}{\alpha }^{pq}(x)`$ and the theta-functions appearing above have complex period $`T=\mathrm{c}_2/\mathrm{c}_1`$. This has the right quasiperiodicity under $`xx+\mathrm{c}_1`$. In order to also have the right quasiperiodicity under $`xx+\mathrm{c}_2`$, $`Q_{1,2}`$ must satisfy

$$Q_1+Q_2=\frac{\stackrel{}{\varphi }_{ab}\times \stackrel{}{\mathrm{c}}}{2\pi }+\frac{\mu }{\mathrm{m}}\stackrel{}{\alpha }\stackrel{}{\mathrm{c}}$$
(34)
and to have the right behavior around $`x=0`$ we must further have

$$\frac{\vartheta _1^{}\left(\frac{\pi }{\mathrm{c}_1}Q_1\right)}{\vartheta _1\left(\frac{\pi }{\mathrm{c}_1}Q_1\right)}+\frac{\vartheta _1^{}\left(\frac{\pi }{\mathrm{c}_1}Q_2\right)}{\vartheta _1\left(\frac{\pi }{\mathrm{c}_1}Q_2\right)}=i\frac{\varphi _{ab}}{\pi }2i\frac{\mu }{\mathrm{m}}\alpha _1$$
(35)
$$A=\frac{\pi ^2\vartheta _1^{}(0)^2}{\mathrm{m}\mathrm{c}_1^2\vartheta _1\left(\frac{\pi Q_1}{\mathrm{c}_1}\right)\vartheta _1\left(\frac{\pi Q_2}{\mathrm{c}_1}\right)}$$
(36)
The equations (34) and (35) above determine $`Q_1`$ and $`Q_2`$, while (36) in turn determines $`A`$. The self-coupling $`\stackrel{~}{W}_\stackrel{}{\alpha }^{pq}`$ can then be extracted from $`W_\stackrel{}{\alpha }^{pq}(x)`$ as

$$\stackrel{~}{W}_\stackrel{}{\alpha }^{pq}=\underset{x0}{lim}\left(W_\stackrel{}{\alpha }^{pq}(x)\frac{1}{\mathrm{m}x^2}\right)$$
(37)
The sums appearing in (28) and (29) are in general convergent, due to the presence of the phases. For $`\omega =1`$, however, the phases are absent and terms with $`p=q`$ have an additive ambiguity due to the need for regularization for the expression (28). In the theta-function expression this manifests in the fact that the equations for $`Q_{1,2}`$ (34,35) are satisfied for any $`Q_1=Q_2`$. By applying the addition formula

$$\vartheta _1(x+Q)\vartheta _1(xQ)\vartheta _4(0)^2=\vartheta _1(x)^2\vartheta _4(Q)^2\vartheta _4(x)^2\vartheta _1(Q)^2$$
(38)
this is seen indeed to amount to an arbitrary additive constant to the expression for $`W^{pp}(x)`$. The same holds for terms $`p,q`$ for which $`\stackrel{}{\varphi }_{pq}=0`$. Such arbitrariness, however, corresponds to trivial redefinitions of the model by addition of constants of motion. This is explained in  and will not be elaborated here.

## 4 Conclusions and open questions

In conclusion, we identified an integrable generalization of the elliptic spin model which breaks the spin $`U(n)`$ invariance and promotes the potential to a modular function introducing noncommutative spin twists. (Quantum generalizations of elliptic spin models have appeared recently, but with $`U(\mathrm{n})`$ invariant interactions .)

There are clearly many issues that deserve further study. The conserved quantities and Lax matrix of this model can, in principle, be obtained from the corresponding quantities of the unreduced model; this was done, e.g., in  for a specific case. A derivation in the present case would be useful. Further, the modular potential involves implicitly defined $`Q_{1,2}`$; a more explicit and symmetric expression would be desirable.

The properties of the spin interaction should also be clarified. In particular, it would be interesting to see if some deformation of $`U(\mathrm{m})`$ can be identified as a dynamical symmetry.

Finally, the quantum mechanical extension of the model is, perhaps, the most interesting and pressing question. This is not trivially accessible by the method used here since, in general, the constraints implied by the phase space restrictions are second class and we cannot carry over the solution of the unrestricted quantum model and apply the constraints as operator relations on the Hilbert space. The above issues remain interesting topics for future investigation.

Acknowledgement: I would like to thank C. Zachos for useful comments on the manuscript.


text:
# General Flavour Blind MSSM and CP Violation

## Abstract:

We study FCNC and $`CP`$ violating processes in the MSSM without a new flavour structure (flavour blind MSSM). The low energy parameters are determined by the running of the soft breaking terms from the GUT scale with SUSY phases consistent with the EDM constraints. We find that the $`CP`$ asymmetry in $`bs\gamma `$ can reach large values potentially measurable at B factories, especially in the low $`BR(bs\gamma )`$ region. We analyze the SUSY contributions to the anomalous magnetic moment of the muon pointing out its impact on the $`bs\gamma `$ $`CP`$ asymmetry and on the SUSY spectrum.

conference: International Europhysics Conference on HEP

In this contribution we study in a systematic way the restrictions on the supersymmetry (SUSY) parameters and complex phases which can be derived from the experimental information on flavour changing neutral currents (FCNC) and on $`CP`$ violation. We choose the Minimal Supersymmetric Standard Model (MSSM) as our theoretical framework and focus on a class of SUSY models that we call flavour blind MSSM. With this term we refer to a model where the soft breaking terms at the grand unification (GUT) scale do not introduce any new flavour structure beyond the usual Yukawa matrices. In this restricted class of models, the number of parameters is largely reduced and it is therefore possible to perform a complete phenomenological analysis. Here we present the main results of our analysis. For details we refer to .

We consider two cases, which are specified by the structure of the soft breaking parameters at the GUT scale. The first case is the simplest version of the constrained MSSM, where we take the following independent parameters:
(I) $`M_{1/2}`$, $`M_0^2`$, $`|A_0|`$, $`\varphi _{A_0}`$, $`\varphi _\mu `$, and $`\mathrm{tan}\beta `$,
which are the universal gaugino and scalar mass parameter, the absolute value and the phase of the trilinear scalar coupling, the phase of $`\mu `$, and the ratio of the vacuum expectation values of the Higgs fields. The second case refers to the SUSY $`SU(5)`$ model, where we take the following set of parameters:
(II) $`M_{1/2}`$, $`M_5^2`$, $`M_{10}^2`$, $`M_{H_1}^2`$, $`M_{H_2}^2`$, $`|A_u|`$, $`|A_d|`$, $`\varphi _{A_u}`$, $`\varphi _{A_d}`$, $`\varphi _\mu `$, $`\mathrm{tan}\beta `$,
where now $`M_5^2`$, $`M_{10}^2`$, $`M_{H_1}^2`$ and $`M_{H_2}^2`$ are the scalar mass parameters of the $`\overline{5}`$ and $`10`$ sfermions and the two Higgs doublets, and $`A_u`$ and $`A_d`$ are the trilinear scalar couplings of the $`10`$ and $`\overline{5}`$.

With the SUSY parameters defined at the GUT scale we determine the soft SUSY breaking parameters at the weak scale by evolving them down with the renormalization group equations (RGE). In our analysis, we have used two–loop RGEs as given in  and one–loop masses as given in . We fix $`|\mu |^2`$ by demanding radiative breaking of the electroweak $`SU(2)_L\times U(1)`$ symmetry. At the weak scale we impose the constraints from direct SUSY and Higgs particle searches  and from the $`\rho `$–parameter, as well as the requirements of colour and electric charge conservation and the lightest SUSY particle (LSP) to be neutral. With these sets of soft SUSY parameters we calculate the electric dipole moment (EDM) of the electron and the branching ratio of $`bs\gamma `$ and compare them with the experimental data, $`|d^e|4.0\times 10^{27}e\text{cm}`$ and $`2\times 10^4BR(bs\gamma )4.5\times 10^4`$. The sets in agreement with the experimental constraints are used to calculate the $`CP`$ asymmetry of $`bs\gamma `$ and the SUSY contributions to $`ϵ_K`$, $`\mathrm{\Delta }M_{B_d}`$ and $`\mathrm{\Delta }M_{B_s}`$. We also calculate the SUSY contributions to the muon anomalous magnetic moment $`a_{\mu ^+}`$ in order to quantify the effect of the recent experimental data on the observables we are interested in. We assume that the SUSY contribution accounts for the difference between the experimental world average for this quantity  and the corresponding SM prediction , $`\delta a_{\mu ^+}=+43(16)\times 10^{10}`$.

We first study the restrictions on the phases from the electron EDM constraint. Due to cancellations between different contributions these restrictions are less stringent and are very different from the case that the EDM constraint is required for each contribution separately (see  and references therein). In our numerical calculations we scan the scalar and gaugino masses at $`M_{GUT}`$ in the range $`100\text{GeV}<M_i<1\text{TeV}`$, the trilinear terms $`0<|A_d|^2<M_{10}^2+M_5^2+M_{H_1}^2`$, $`0<|A_u|^2<2M_{10}^2+M_{H_2}^2`$, taking their phases arbitrary, and $`4<\mathrm{tan}\beta <50`$. We find that the EDM constraint implies correlations between $`\varphi _\mu `$ and $`\varphi _{A_0}`$ and between $`\varphi _\mu `$ and $`M_0^2`$. While it is possible to find any value for $`\varphi _{A_0}`$, the values of $`\varphi _\mu `$ are more constrained: $`\varphi _\mu <0.1`$ for lower values of $`M_0^2`$, whereas values up to $`\varphi _\mu =0.4`$ are allowed for higher $`M_0^2`$.

Next we discuss the main features of the SUSY particle spectrum relevant for $`CP`$ violating and FCNC observables. In Fig. 1 we show the scatter plot of the lightest chargino versus the lightest stop masses for the SUSY $`SU(5)`$ model (set (II)). All points fulfill the EDM and $`BR(bs\gamma )`$ constraints, whereas the black dots are the points of the parameter space that also fulfill the $`\delta a_{\mu ^+}`$ constraint. We confirm the presence of an upper bound on the chargino mass of about $`700\text{GeV}`$ for very large $`\mathrm{tan}\beta `$ (of order 50) , and lower for smaller values of $`\mathrm{tan}\beta `$, which is a consequence of the $`\delta a_{\mu ^+}`$ constraint. It can be shown that this bound is essentially due to our assumption of gaugino mass unification. We would like to stress, however, that in any RGE evolved MSSM with gaugino mass unification the $`a_{\mu ^+}`$ constraint has important consequences on the complete MSSM spectrum. As can be seen in Fig. 1, there is a strong correlation among the stop and chargino masses due to the dominant gaugino RGE effects. Then, without any further restriction on the SUSY parameter space, this correlation implies the presence of an upper bound on the light stop mass of $`m_{\stackrel{~}{t}}1500\text{GeV}`$.

An especially interesting observable is the $`CP`$ asymmetry in the partial width,

$$A_{CP}^{bs\gamma }=\frac{BR(\overline{B}X_s\gamma )BR(BX_s\gamma )}{BR(\overline{B}X_s\gamma )+BR(BX_s\gamma )}.$$

The SM prediction for this asymmetry is exceedingly small, therefore, it is sensitive to the new SUSY phases. As $`\varphi _\mu `$ and $`\varphi _A`$ are associated with chirality changing operators, we can expect large effects in $`bs\gamma `$ . In Fig. 2 we show the correlation among $`\delta a_\mu `$ and $`\text{BR}(bs\gamma )`$. All points shown fulfill the EDM constraint, whereas the black dots represent the points of the parameter space that also reproduce the measured anomalous magnetic moment. In this plot we can see that the $`\delta a_\mu `$ constraint cuts out the region of high branching ratios $`BR(bs\gamma )`$. In our flavour blind scenario the main contributions to the $`bs\gamma `$ decay are due to charged Higgs, light stop and light chargino loops. The $`BR(bs\gamma )`$ constraint and the $`\delta a_\mu `$ constraint favour $`Re(\mu )>0`$. Imposing the EDM constraint reduces the $`CP`$ asymmetry to less than $`1\%`$ in the region of large branching ratio $`BR(bs\gamma )`$. However, for small $`BR(bs\gamma )`$ the $`CP`$ asymmetry can go up to $`5\%`$, making it potentially measurable at B factories.

Finally we calculated the SUSY contributions to the $`CP`$ violating quantities $`ϵ_K`$, $`\mathrm{\Delta }M_{B_d}`$, $`\mathrm{\Delta }M_{B_s}`$ and studied the modifications of the unitarity triangle. We obtain only small deviations from the SM predictions. This result is different from those obtained in the context of Minimal Flavour Violation MSSM or in the MSSM with non–universal parameters

Acknowledgements: This work was supported by the ‘Fonds zur Förderung der wissenschaftlichen Forschung’ of Austria, FWF Project No. P13139-PHY, by the EU TMR Network Contracts No. HPRN-CT-2000-00148, HPRN-CT-2000-00149 and HPRN-CT-2000-00152 and by Cooperazione scientifica e tecnologica Italia-Austria 1999-2000, Project No. 2; W.P. was supported by the Spanish Ministry of Education and Culture under the contract SB97-BU0475382 and by DGICYT grant PB98-0693; T.G. acknowledges financial support from the DOE grant DE-FG02-96ER40967; E.L. and O.V. thank SISSA for support during the first stages of this work. E.L. was supported by the Alexander Von Humboldt Foundation; O.V. acknowledges financial support from a Marie Curie E.C. grant (HPMF-CT-2000-00457) and partial support from Spanish CICYT AEN-99/0692.


text:
# Hydrodynamical Models of Pulsar Wind Bow-Shock Nebulae

## 1. Introduction

Pulsars moving supersonically with respect to the ambient medium are expected to give rise to a bow shock. In the case of interaction with the interstellar medium (ISM) the H$`\alpha `$ emission from the nebula may be detected. Only four bow shocks have been discovered so far: PSR 1957+20 (Kulkarni & Hester 1988), PSR 2224+65 (Cordes, Romani & Lundgren. 1993), PSR J0437+4715 (Bell 1995), and PSR 0740–28 (Jones, Stapper & Gaensler 2001). These nebulae may be detected in optical Balmer lines (mostly H$`\alpha `$ ) as a signature of a non-radiative shock moving through a partially neutral medium (Chevalier & Raymond 1980).

A puzzling point is that two out of the four known nebulae show a peculiar shape with a conical tail. Various hypothesis have been put forward to justify the conical tail shape: a peculiar ISM density distribution in the surroundings of the PSR 2224+65; effects like mass loading due to neutral atoms which might penetrate the external layers of shocked material (Bucciantini & Bandiera 2001, hereafter Paper I) . The nebula near PSR 0740–28, seems to represent an intermediate case, with a “standard” head and a conical tail. The ISM is seen by the pulsar as a plane-parallel flow, likely with a constant density, and lasting long enough to produce a steady-state regime. Moreover the ISM is typically partly neutral and the neutral atoms may have collisional mean free paths comparable with the scale length of the system, then invalidating a purely fluid treatment. However as demonstrated in a previous paper (Paper I), if we suppose that the H atoms can be ionized only via collisions, for a large amount of pulsars the presence of a neutral component in the ISM can be neglected .

The interaction of the relativistic magnetized pulsar wind with the ISM ionized component is due to the pulsar wind magnetic field compressed on the head of the nebula. The mean free path of particles is typically much smaller than the typical dimension of the nebula: $`d_o=\sqrt{/4\pi c\rho _oV_o^2}`$,where $``$ and $`V_o`$ are the pulsar luminosity and velocity, and $`\rho _o`$ the density of the dragged component of ambient medium (Paper I).

## 2. Numerical Simulation

We have used the multi-D code CLAWPACK (LeVeque 1994): a second order Godunov with piecewise linear reconstruction. The problem has an intrinsic cylindrical symmetry which allows us to use a two-dimensional grid. Both the stellar and the external winds are taken with high Mach numbers (hypersonic limit). In fact the pulsar wind has a thermal pressure much lower than its ram pressure and the same holds true for the ISM. We decided to use a different adiabatic coefficient for the ambient medium (with adiabatic coefficient $`5/3`$) and the relativistic ($`4/3`$) fluid coming from the star. The physical conditions of such a bow-shock guarantee conservation of energy, and little mixing between the two media. However we added a little diffusion to avoid the growth of numerical instabilities. Ul is the unit length we have used. For the simulations we used a grid of 580$`\times `$400 cells and the star in (420,0). The stellar wind has a velocity of 1000 km s-1, a density of 0.04 cm-3 and a Mach number about 13 upstream of the termination shock along the bow shock axis. We have performed various simulations using different velocities for the external flow (400, 300, 200, 150 km s-1). The time needed to reach the steady state condition is of order 4-5 times the external flow crossing time in the computational box.

## 3. Comparison with the analytic model

The geometry and internal structure of bow-shock nebula resulting from our simulations are shown in Fig. 1. The positions of the termination shock, bow shock and contact discontinuity do not change among all the above mentioned simulations, and also the position of the sonic point remains essentially the same. This is because the solution, in the hypersonic limit, depends only on the ratio between ram pressures of the stellar and ambient medium, but not separately on their densities or velocities. The dashed lines represent two analytic solutions (Wilkin 1996), one scaled to match the contact discontinuity near the head, and the other to match the bow-shock near the head. Fig. 1 shows that the region of unperturbed stellar wind extends in the tail up to a “spherical shock”, whose position is strictly connected with the value of the thermal pressure in the ambient medium. The contact discontinuity in the tail (Fig. 1) tends to reach a limiting distance $`Z`$ from the axis, with radius of the same order of the stagnation point distance. In order the stellar wind material to expand sideway beyond $`Z`$, some process is required, which reduces the velocity in the layer or increases its pressure. This can be achieved via processes like mass accretion, or if some further shock develops. The distance of the bow-shock from the star given by the analytic model is 0.175 Ul while the two limiting analytic solutions of Fig. 1 (dashed lines) are at 0.22 and 0.29 Ul. The behaviour of the fluid in regions A and B may affect the synchrotron emission from the relativistic electrons and positrons coming from the pulsar and shocked in the termination shock. If particles have relativistic motions their synchrotron emission can be detected only if their pitch angle is close to the observer direction. In the tail the relativistic bulk motion produce a beaming of the emission in the backward direction.

To verify if the two-thin layers model could be taken as a reasonable approximation of the true structure of the nebula, at least as far as the external layer is concerned, and to evaluate how far in the tail it could be used, we have extracted from our simulations the surface density and the tangential velocity in the external layer, between the bow shock and the contact discontinuity. Fig. 2-left compares the surface density in the external layer, extracted from the simulation, with that evaluated by a two-thin layers model. The analytic curve that best reproduces the numerical points gives an external density of 0.26 cm-3 (instead of the true value 0.25 cm-3) and a stagnation point distance of 0.235 Ul, compatible with the range of the two limit solutions.

In the same way we have evaluated the tangential velocity in the layer between the bow shock and the contact discontinuity. Fig. 2-right shows the numerical values of tangential velocity in the external layer, normalized to the external wind value, and compares it to the two-thin layers model. The analytic curve that reproduces the points has been obtained using a values of 0.25 Ul for the distance of the stagnation point, a value in agreement with the two limit solutions.

## 4. Magnetic field and penetration thickness for neutral hydrogen

Using the numerical hydrodynamic solutions described above we have also evaluated the behaviour of a toroidal magnetic field wounded around the symmetry axis and passively advected by the fluid. We considered the case of a pulsar wind with the energy associated to the ordered magnetic field much smaller than the ram pressure of the relativistic particles, and with the ratio between the two being independent of the direction. The value of such a field increases in the head, due to the compression in the termination shock, ranging from 7 times its upstream value, suddenly after the shock, up to a factor 30-35 with respect to its upstream value, near the contact discontinuity. Even if initially the magnetic field’s energy is below equipartition, it can reach equipartition if the external wind Mach number is high enough.

As we have demonstrated in a previous paper (Paper I) the penetration thickness $`\tau `$ for neutral hydrogen in the external layer, depends on ambient density and pulsar velocity. The penetration thickness depends by the processes that could drag the atoms: essentially by charge-exchange and by ionization via the interaction with the protons confined in the external layer. In our previous paper we have evaluated the thickness using the two-thin layers model limited to a region near the axis. Our aim is to verify if such approximation is acceptable and how far in the tail it can be extended. For typical pulsar velocities, collisional ionization can be neglected, even in non equilibrium cases. The thickness decreases at higher velocities as expected: in fact it follows the behaviour of the cross section for charge-exchange, which is the dominant process. The thickness, at least for the head, increases slightly: only over a factor 2. This can be quite important, because, if the nebula is thin to the penetration of neutral hydrogen in the head, it can remains thin even in tail.

## References

Bell, J. F., Bailes, M., Manchester, R. N., Weisberg, J. M., & Lyne, A. G. 1995, ApJ, 440, L81

Buciantini, N., & Bandiera, R. 2001, A&A, 375, 1032

Chevalier, R. A., & Raymond, J. C. 1978, ApJ, 225, L27

Cordes, J. M., Romani, R. W., & Lundgren, S. C. 1993, Nature, 362, 133

Jones, H., Stappers, B., & Gaensler, B. 2001, The Messenger, 103, 27

Kulkarni, S. R., & Hester, J. J. 1988, Nature, 335, 801

LeVeque, R. J. 1994, Proceedings of the Fifth International Conference on Hyperbolic Problems.

Wilkin, F. P. 1996, ApJ, ApJ, L31