text:
# Diomega production in relativistic heavy ion collisions

## Abstract

Using a multiphase transport model, we study the production of a new strange dibaryon $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ in dense hadronic matter formed in relativistic heavy ion collisions. The (multi-)strange baryons ($`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$) are produced by strangeness-exchange reactions between antikaons and hyperons in the pure hadronic phase. The rescattering between the $`\mathrm{\Omega }`$s at midrapidity leads to a production probability of $`3\times 10^7`$ $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ per event at the RHIC energy of $`\sqrt{s}=130A`$ GeV. The production probability would be enhanced by two orders of magnitude if $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ and $`\mathrm{\Omega }`$ reach chemical equilibrium during heavy ion collisions. We further find that the yield of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ increases continuously from SPS to the highest RHIC energy.

PACS numbers: 25.75.-q, 24.10.Lx, 14.20.Pt

Since first proposed by Jaffe  that the H dibaryon with quark content ($`uuddss`$) could be a possible lightest strangelet $``$ droplet of bound strange quark matter, properties of dibaryons have been widely investigated in various models, such as the MIT bag model , the Skyrme model , and the constituent quark model . The consensus regarding the mass of H dibaryon is about 2.232 GeV . However, extensive experimental searches  have not identified any strangelets with small values of charge fraction, $`f_Z=|Z|/A<1`$, and strangeness fraction $`f_S=|S|/A<2`$.

Recently, the structure and properties of dibaryons with large strangeness are investigated in the chiral SU(3) model  that have been quite successful in reproducing several nuclear properties . In this model the H dibaryon is found to be only weakly bound . On the other hand, analysis of some six-quark cluster states with high strangeness fraction reveals that the diomega $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$, in particular, is rather deeply bound. Although the color magnetic interaction in the one gluon exchange term for this system exhibits repulsive feature, the large attraction stemming from the chiral quark coupling and from the symmetry property of the system lead to a rather large binding. In the Resonating Group Method calculation, the binding energy of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ is found to be as large as $`116`$ MeV, and the root mean square distance between the two $`\mathrm{\Omega }`$’s is 0.84 fm. Besides the large (negative) charge fraction $`f_Z=1`$ and strangeness fraction $`f_S=3`$, the new dibaryon $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ has quite a long mean lifetime of $`10^{10}`$ sec as it can undergo only weak decay.

Because of its large strangeness, $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ is not likely to be produced in proton-proton collisions. On the other hand, strangeness production is enhanced in heavy ion collisions and has been suggested as one of the possible signals of quark-gluon plasma (QGP) due to large gluon density and low energy threshold for $`s\overline{s}`$ formation . This may thus lead to the formation of exotic deeply bound objects composed of quarks or baryons with large strangeness. Therefore, in relativistic heavy ion collisions, especially at RHIC energies, the dibaryon $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ could be a new interesting candidate.

Indeed, recent measurements by the WA97 Collaboration  and the NA49 Collaboration  demonstrated substantial enhancement of the (anti-)hyperon yields ($`\mathrm{\Lambda }`$, $`\mathrm{\Xi }`$, and $`\mathrm{\Omega }`$) in $`158A`$ GeV Pb-Pb central collisions relative to p-Pb collisions. The enhancement pattern increases with the strangeness content of the (anti-)hyperon. Such a large enhancement for multistrange baryons at midrapidity was interpreted as a signal for quark-gluon plasma formation.

On the other hand, even without a phase transition (anti)strangeness can be abundantly produced by hadron rescatterings alone. In fact, within the microscopic transport model UrQMD, the WA97 data for multistrange baryon enhancement can be explained by reducing the constituent quark mass in the fragmentation of the initial strings in dense matter or by increasing the string tension . Based on the rate equation approach, the multimesonic reactions $`\overline{Y}+Nn\pi +n_YK`$ was demonstrated  to enhance antistrange hyperon $`\overline{Y}`$ production in a hadronic scenario. Using a multiphase transport (AMPT) model, we found that strangeness-exchange reactions between antikaons and hyperons in a pure hadronic stage also leads to a significant production of multistrange baryons at the SPS and RHIC energies .

The considerable production of (multi-)strange baryons in the dense hadronic stage, in absence of QGP formation, could then also result in the formation of exotic objects with enhanced strangeness. Properties of metastable multistrange baryonic objects consisting of nucleons, $`\mathrm{\Lambda }`$’s and $`\mathrm{\Xi }`$ are studied in the relativistic mean-field theory . These objects were found to have properties quite similar to those of strangelets $``$ their quark counterpart. Estimates on the production of such strange dibaryons were presented  at the RHIC energies by employing wave-function coalescence in the RQMD model.

The pronounced $`\mathrm{\Omega }`$ production obtained in the AMPT model from strangeness-exchange reactions suggests that multiple collisions between these omegas in the dense hadronic matter may lead to an appreciable production of the new strange dibaryon $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$. In this letter, we investigate $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ production in relativistic heavy ion collisions in the AMPT model.

The AMPT model  is a hybrid model that uses as input both the minijet partons from the hard processes and the strings from the soft processes in the HIJING model . The dynamical evolution of partons are modeled by the ZPC  parton cascade model, while the transition from the partonic matter to the hadronic matter is based on the Lund string fragmentation model . The final-state hadronic scatterings are then modeled by the ART model . The AMPT model has been very successful in describing the measured transverse momentum and rapidity distributions of charge particles  as well as the particle to antiparticle ratios . The multistrange baryon production is via the strangeness-exchange reactions $`\overline{K}\mathrm{\Lambda }\mathrm{\Xi }\pi `$, $`\overline{K}\mathrm{\Sigma }\mathrm{\Xi }\pi `$ and $`\overline{K}\mathrm{\Xi }\mathrm{\Omega }\pi `$. Since there is no experimental information on these cross sections, they are assumed to be the same as that for $`\overline{K}N\mathrm{\Sigma }\pi `$, which is the isospin averaged cross section for converting a nucleon to a sigma . This assumption is consistent with results obtained using SU(3) invariant hadronic Lagrangians . For diomega production from collisions between omegas, we consider both the the electromagnetic process $`\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })_{0+}+\gamma `$ and the strong interaction process $`\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })_{0+}+\eta `$. Their cross sections have been evaluated in Ref.  using an effective Hamiltonian. For the process $`\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })_{0+}+\gamma `$, the maximum cross section is found to be $`1.6\mu \mathrm{b}`$ at $`p_\mathrm{\Omega }0.5`$ GeV, while that for the process $`\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })_{0+}+\eta `$, which requires a threshold $`p_\mathrm{\Omega }>0.88`$ GeV due to finite $`\eta `$ mass, reaches a maximum value of $`5.5\mu \mathrm{b}`$ at $`p_\mathrm{\Omega }1.1`$ GeV. Because of the mismatch between the large initial and small final relative momenta between the two omegas in the reaction $`\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })_{0+}+\eta `$, which involves only strong interactions, its cross section is only slightly larger than that for the reaction $`\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })_{0+}+\gamma `$, which involves a much weaker electromagnetic interaction.

To include the diomega production reactions in the AMPT transport model, we parametrize their cross sections as functions of center of mass energy $`\sqrt{s}`$ (in GeV), i.e.,

$$\sigma _{\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })+\gamma }=87\left(1\frac{\sqrt{s}_0}{\sqrt{s}}\right)^{0.91}\left(\frac{\sqrt{s}_0}{\sqrt{s}}\right)^{24.6}\mu \mathrm{b},$$
(1)
where $`\sqrt{s}_0=2m_\mathrm{\Omega }=3.345`$ GeV, and

$$\sigma _{\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })+\eta }=1.5\times 10^4\left(1\frac{\sqrt{s}_0}{\sqrt{s}}\right)^2\left(\frac{\sqrt{s}_0}{\sqrt{s}}\right)^{34.9}\mu \mathrm{b},$$
(2)
where $`\sqrt{s}_0=m_\eta +2m_\mathrm{\Omega }B_{(\mathrm{\Omega }\mathrm{\Omega })}=3.777`$ GeV.

Due to small number of multistrange baryons in heavy ion collisions, the productions of $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$ as well as of the diomega are treated perturbatively in the AMPT model , i.e., the collision dynamics is assumed to be unaffected by the production of these particles. For $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ production, since its abundance is extremely small, we assume that the production/annihilation rate of $`\eta (\gamma )`$ is unaffected. Taking a total scattering cross section between two $`\mathrm{\Omega }`$s to be $`\sigma _{\mathrm{\Omega }\mathrm{\Omega }}=20`$ mb, the two $`\mathrm{\Omega }`$s will collide if their distance of closest approach is $`\sqrt{\sigma _{\mathrm{\Omega }\mathrm{\Omega }}/\pi }`$. If $`P_\mathrm{\Omega }`$ and $`P_\mathrm{\Omega }^{}`$ are the production probabilities of the colliding $`\mathrm{\Omega }`$s, the probability of the resulting diomega is then $`P_{(\mathrm{\Omega }\mathrm{\Omega })_{0+}}=P_\mathrm{\Omega }P_\mathrm{\Omega }^{}\mathrm{\Gamma }_{\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })_{0+}}`$, where $`\mathrm{\Gamma }_{\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })_{0+}}=(\sigma _{\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })+\gamma }+\sigma _{\mathrm{\Omega }+\mathrm{\Omega }(\mathrm{\Omega }\mathrm{\Omega })+\eta })/\sigma _{\mathrm{\Omega }\mathrm{\Omega }}`$. The probabilities of each of the interacting omegas are reduced by $`P_{(\mathrm{\Omega }\mathrm{\Omega })_{0+}}`$. The annihilation cross section of a diomega via the reaction $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}+\eta \mathrm{\Omega }+\mathrm{\Omega }`$ is obtained using detailed balance of Eq. (2). If a diomega collides with a eta, the diomega is totally annihilated and two omegas are produced with isotropic momentum distribution and each having a probability corresponding to the diomega production probability of $`P_{(\mathrm{\Omega }\mathrm{\Omega })_{0+}}`$. Since diomega production is treated perturbatively, the final results do not depend on the value of $`\sigma _{\mathrm{\Omega }\mathrm{\Omega }}`$ used in the AMPT model.

In Fig. 1, we show the time evolution of the abundances of midrapidity $`\mathrm{\Omega }^{}`$ and $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ in central Au+Au collisions at the RHIC energy of $`\sqrt{s}=130A`$ GeV. As is evident from the figure, most of the omegas and $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ are produced within 10 fm/c after hadronization when the energy density is high. Beyond this time the diomega yield saturates to a value of $`2.8\times 10^7`$ per event.

From the study of time evolution of the production and absorption rates of particles, we found  that most strange hadrons, including the $`\mathrm{\Omega }`$s, reach chemical equilibrium in collisions at the RHIC energy of $`\sqrt{s}=130A`$ GeV. In Fig. 1 we depict the time evolution of the production (solid lines) and absorption rates (dashed lines) of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$. It is observed that the absorption rate of a diomega by colliding with a eta meson is always smaller compared to its production rate at all times. This stems from additional cross section in the production channel due to the electromagnetic process which has no threshold. We thus infer that the small number of diomega $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ that may be produced in the relativistic collisions at RHIC do not reach chemical equilibrium.

The rapidity distributions of the produced hadrons peak at midrapidity  because of increased multiple collisions and rescattering at central rapidities compared to that at large rapidities. In Fig. 2 we show the rapidity distribution of the multistrange particles $`\mathrm{\Xi }`$ and $`\mathrm{\Omega }`$, along with that of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$. The rapidity distribution of particles with increasing strangeness content gradually become narrower as these particles are successively produced from their parent that collide more frequently at small rapidities where the baryon density is high.

We have also studied the energy dependence of diomega production in relativistic heavy ion collisions. In Fig. 3, we show the yield of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ at midrapidity in Au+Au collisions at various energies. Also shown in the figure are the yields of $`\mathrm{\Omega }^{}`$ and $`\mathrm{\Xi }^{}`$ as functions of energies. It is seen that while the $`\mathrm{\Omega }^{}`$ number increases by about a factor of three from the SPS energy ($`\sqrt{s}=17A`$ GeV) to the highest RHIC energy ($`\sqrt{s}=200A`$ GeV), the $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ number increases by about a factor of 4, indicating that the dibaryon $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ number reveals a much faster rate of increase with $`\sqrt{s}`$ as compared to $`\mathrm{\Omega }^{}`$.

The AMPT model gives a lower bound on the diomega production probability since there might be other production channels, such as $`\mathrm{\Omega }^{}+\mathrm{\Xi }^{}(\mathrm{\Omega }\mathrm{\Omega })_{0+}+K^0`$, which have been neglected in the present calculation because of their unknown cross sections. If the diomega production cross section is larger, then diomegas may reach chemical equilibrium with omegas in heavy ion collisions. For most other particles including multistrange baryons, equilibrium thermal models have been successfully employed  to explain experimental data for their yields and ratios over a wide range of energies from AGS to RHIC. Adopting the statistical model that is based on the grand canonical ensemble with complete thermal, chemical, and strangeness equilibrium, we have found that the results from the AMPT model for the $`K^+/\pi ^+`$, $`\overline{p}/p`$, and $`K^{}/K^+`$ ratios, which are about 0.18, 0.65, and 0.89, respectively, at midrapidity in Au+Au collisions at energy $`\sqrt{s}=130A`$ GeV, can be approximately described with a temperature $`T170`$ MeV, baryon chemical potential $`\mu _B37`$ MeV, and strange chemical potential $`\mu _S10`$ MeV. If we assume that diomegas are also in chemical equilibrium with omegas, then the ratio $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}/\mathrm{\Omega }^{}`$ is $`7.4\times 10^5`$. With the omega number of about 0.41 at midrapidity, this leads to a diomega production probability of $`3.0\times 10^5`$ per event, which is two orders of magnitude higher than that obtained in our transport model.

Diomega production in heavy ion collisions can also be studied using the coalescence model  based on the omega phase space distribution at freeze out as obtained in the AMPT model. Assuming a harmonic oscillator wave function , the Wigner density for the diomega is $`\rho _D(𝐫,𝐪)=8\mathrm{exp}(r^2/d^2q^2d^2)`$, where $`𝐫=𝐫_1𝐫_2`$ and $`𝐪=(𝐩_1𝐩_2)/2`$ are given in the c.m. system of $`\mathrm{\Omega }`$-$`\mathrm{\Omega }`$, and the value of $`d=0.69`$ fm corresponds to a rms radius of 0.84 fm for the diomega. For central Au+Au collisions at $`\sqrt{s}=130A`$ GeV, the yield of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ at midrapidity is found to be $`2.6\times 10^5`$ and is comparable to that from the thermal model. The similar results obtained from both the thermal and coalescence models are not surprising as it was shown in Ref. , that the two models are equivalent when matter is in thermal and chemical equilibrium and the binding energy of the composite particle is much smaller than the temperature. Since the AMPT model predicts a hadronic matter at freeze out that is close to thermal and chemical equilibrium, the diomega yield from the coalescence model thus should be similar to that given by the thermal model.

Our estimate for the production probability of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ are well within the limits of the present detectors used at RHIC energies. Therefore, this exotic object can, in principle, be detected in present and future experiments. The fact that $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ has a large strangeness content with a binding energy of $`116`$ MeV, it is stable against strong hadronic decays and possess weak decays: $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}\pi ^{}+\mathrm{\Xi }^0+\mathrm{\Omega }^{}`$ and $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}\pi ^0+\mathrm{\Xi }^{}+\mathrm{\Omega }^{}`$. Because three-body decay are involved, the final state phase-space would be suppressed. In the sudden approximation, the mean lifetime of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ was found  to be about four times longer than the free $`\mathrm{\Omega }`$ lifetime of $`0.822\times 10^{10}`$ sec. Apart from these conventional decay modes, the nonmesonic decay $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}\mathrm{\Xi }^{}+\mathrm{\Omega }^{}`$ is also possible; the estimated  lifetime of $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ for this process is twice the free $`\mathrm{\Omega }`$ lifetime. Thus, instead of direct observation, the $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ may also be detected in the $`\mathrm{\Xi }^{}\mathrm{\Omega }^{}`$ invariant mass distribution. The observation of which could provide useful information of the unknown $`\mathrm{\Omega }\mathrm{\Omega }`$ interaction strength. Our study thus opens up the intriguing possibility of detecting the new dibaryon $`(\mathrm{\Omega }\mathrm{\Omega })_{0+}`$ in heavy ion collisions at the RHIC energies.

This work of SP and CMK is supported by the National Science Foundation under Grant No. PHY-9870038, the Welch Foundation under Grant No. A-1358, and the Texas Advanced Research Program under Grant No. FY99-010366-0081, while that of ZYZ is supported by the National Natural Science Foundation of China and the Chinese Academy of Sciences.


format: text

text:
# 1 Introduction

## 1 Introduction

Boolean networks were introduced by Kauffman  as simplified models of the complex interaction in the regulatory networks of living cells. The binary variable $`\sigma _i`$ encodes the activity of the effective “gene” $`i`$; expressed or not expressed. Depending upon the initial state, the system evolves to one of the several limit cycles. In the biological picture, the different limit cycles are interpreted as different cell types. One of Kauffman’s motivations for investigating these networks was the idea that the structure of genetic networks present in nature is not only determined by selection. Rather, a good fraction of the network functionality is inherent in the ensemble of regulatory networks as such. In fact, he found an ensemble of critical Boolean networks “on the edge of chaos” that captures some features observed in nature. These Boolean networks show a remarkable stability; in most cases small perturbations in the state of the network do not change the trajectory to a different limit cycle. This is desirable in the biological interpretation since stability of genetic regulatory networks against small fluctuations is a crucial property. Another striking observation is that the number of limit cycles for the critical Boolean networks grows as a square-root of the system size . This is an analogy to multicellular organisms, where it is found empirically that the number of cell-types also grows approximately as the square-root of the genome-size.

The model also exhibits analogies  with infinite range spin glasses . In the framework of an annealed approximation , some of the previous numerical observations concerning a phase transition between a frozen and a chaotic phase in the model could be understood. The average limit Hamming distance $`d_h`$, the number of bit-wise differences between two random configurations, was used as an order parameter. In the frozen phase one has $`d_h=0`$ for infinite systems whereas in the chaotic phase one has $`d_h0`$. The parameter driving the transition is the probability $`p`$ that two different inputs to a Boolean variable $`\sigma _i`$ give raise to different values. In  the annealed approximation was extended to provide distributions for the number and the length of the limit cycles. Also, good agreement between the results from the annealed approximation and the numerical calculations was demonstrated in the chaotic phase. An alternative order parameter $`s`$, the fraction of variables that are stable, i.e. evolve to the same fixed state independently of the initial state, was introduced in . In the infinite size limit one has $`s=1`$ in the frozen phase, whereas $`s1`$ in the chaotic phase. In  the concept of relevant variables was introduced. A variable $`\sigma _i`$ is not relevant, if it is stable and/or no variable state depends on $`\sigma _i`$. The relevant variables are of interest since they contain all information about the asymptotic dynamics of the network, i.e. the number of limit cycles and their cycle lengths.

In this work we focus on the stability of the Kauffman model and how this property is related to the stable core of the network. The probability that inversion of a single variable will make the system end up in a different limit cycle is known to be small and approaches zero for large networks. However, we find that if the network is reduced to its relevant core, this probability in drastically raised and increases slightly with the system size.

To facilitate this study we introduce a method that removes variables that cannot be relevant by inspection of transition functions and network connectivity. The resulting reduced network contains all relevant variables and possibly some irrelevant ones. Since all relevant variables are included it will have exactly the same asymptotic dynamics as the original network even though the total number of variables is drastically reduced. We find that resulting number of variables is close to the true number of relevant variables. This indicates that properties of the stable core can mostly be understood by the comparatively trivial interactions detected in the decimation procedure.

The decimation procedure can also be used to reduce the bias in the estimate of some observables like for example the number $`n_c`$ of limit cycles. Different from earlier works we do not observe a $`\sqrt{N}`$ scaling, but rather a linear growth of $`n_c`$ with the system size.

## 2 The Kauffman Model

A random Boolean network essentially is a cellular automaton with $`N`$ binary state variables $`\sigma _i`$. These evolve synchronously according to the transition functions $`f_i\left(\left\{\sigma \right\}\right)`$, which are chosen randomly at time $`t=0`$ and are then kept fixed. In the Kauffman model $`f_i`$ are constrained to depend on at most $`K`$ different randomly chosen input variables:

$$\sigma _i\left(t\right)=f_i[\sigma _{v_i^1}\left(t1\right),\mathrm{},\sigma _{v_i^K}\left(t1\right)],$$
(1)
for every variable $`\sigma _i`$. The integers $`\{v_i^1,\mathrm{},v_i^K\}`$ define the input connections to variable $`i`$.

The transition function $`f_i`$ maps each possible combination of input signals to Boolean output values. These output values are independently set to $`true`$ or $`false`$ with probabilities $`p`$ and $`1p`$ respectively. This makes some functions independent of some or all of its $`K`$ input variables. Furthermore, depending on $`K`$ and $`p`$, a finite fraction of the state variables $`\sigma _i`$ are not used by any of the transition functions.

The random Boolean network is a deterministic system. Given the state variables at some time, the future trajectory of the $`\sigma _i`$ is known. The volume of the state space is finite, therefore all trajectories must posses a limit cycle. Besides the stability of the system the number of limit cycles, the length distribution of the cycles and transient trajectories are well established observables for this model. In numerical simulations it is in general not possible to probe the models whole state space, except for very small systems. The volume of the $`\left\{\sigma \right\}`$ state space grows exponentially and the number of graphs $`\left\{f_i\right\}`$ even super-exponentially. The commonly used strategy for exploring this model therefore contains two approximations.

1. A small fraction of all possible networks is used as a representative ensemble.
2. For each network only a subset of the state space is probed.

Point 2 introduces a systematic bias to the number of limit cycles since not all of them will be found. In the results section we will re-analyze the number of cycles after decimation of irrelevant nodes. This allows full enumeration of state space for up to $`N=32`$. In this way we get an improved estimate for the scaling of the respective observables with the system size.

## 3 The Decimation Procedure

It is well known that some variables in a Kauffman network evolve to the same steady state independently of the initial configuration. These stable variables are clearly irrelevant for the asymptotic behavior of the network. The same holds for those variables that do not regulate any other variable, i.e. no transition functions is dependent on them. In fact, as pointed out in , for a variable to be relevant it has to be unstable and regulate some unstable variables that in turn regulate others and so on. In other words, a variable is relevant if and only if it is unstable and regulates other relevant variables.

By definition, in the frozen phase the fraction of stable variables goes to unity as $`N`$ goes to infinity. Therefore, a large fraction of the variables are likely to be stable even for finite $`N`$. Since the irrelevant variables includes all stable variables, a considerable part of a network does not affect the asymptotic dynamics at all. The process of identifying the irrelevant variables can be divided into two separate steps. Firstly, the stable variables are identified. Secondly, the variables that do not regulate any unstable variables are identified.

Identifying stable variables is in principal easy, but computationally demanding. In  this was done by performing simulations of the dynamics of the system and monitoring which variables were in the same state in all probed limit cycles. However, finding all limit cycles essentially means that all $`2^N`$ possible states have to be probed, which is possible only for very small networks. Since a variable that is stable within the probed limit cycles may change state within some of the unprobed limit cycles, searching a fraction of state space will in some cases overestimate number of stable variables.

Here we introduce an alternative method, which by pure inspection of the connectivity and the transition functions of a network identifies variables that must be stable. The basis for our approach is that transition functions dependent on no input variables give a constant output, i.e. the corresponding variable is stable.

As stated above, some transition functions are independent of all their input variables, i.e. they are constants. This means that the corresponding variables will be stable (after the initial time step) and a transition function that is dependent on such variable will receive a constant signal. By replacing the stable input variable with the corresponding constant value, the number of input variables is reduced. For each replaced variable the input state space is reduced by a factor $`1/2`$ and within this subspace the rule may be independent of yet other input variables. If in the end even this rule become a constant, the corresponding variable is stable (after a transient time), and can be replaced by a constant. Therefore, we have to repeat this procedure until no more stable variables are found. We summerize the method as follows:

1. For every updating function, $`f_i`$, remove all inputs it does not depend upon.
2. For those $`f_i`$ with no inputs, clamp the variable $`\sigma _i`$ to the corresponding constant value.
3. For every $`f_i`$, replace the clamped inputs with the corresponding constant.
4. If any variable has been clamped, repeat from step 1.

For a pseudo-code description of this method see Appendix (A).

It is clear that our method sometimes does not find all stable variables. We see an example of such a situation in Fig. 1. Here the inputs to a function are coupled logically and hereby confined to a subspace of possibilities. Within this subspace the otherwise unstable variable is stable. The figure illustrates just one of the possible couplings between inputs.

Once the stable variables are identified and removed from the network the non-regulating variables can be removed iteratively. Since our method keeps all relevant variables the resulting network will have exactly the same asymptotic dynamics as the original network.

## 4 Results

Let us start by analyzing the size of the stable core as a function of the system size $`N`$. In Fig. 2 the the size of the stable core $`N^{}`$, identified by the decimation procedure described above, is shown. Each data point is averaged over $`10^4`$ instances of networks. For comparison the size of the stable core $`N^+`$ as estimated by the method used in  is also plotted. The latter procedure is based on observations of the dynamics of the full network and identification of nodes acquiring the same constant value independently of the start-configuration. Since only a small part of the state space can be probed in practice, the number $`N^+`$ is biased to overestimate the true size $`\eta `$ of the stable core. On the other hand, our decimation procedure underestimates $`\eta `$ because some configurations, like the one depicted in Fig. 1, which may lead to stable variables are not identified. Therefore, we have $`N^+\eta N^{}`$.

It is somewhat surprising to observe $`N^+N^{}`$, which indicates that properties of the stable core can mostly be understood by the comparatively trivial interactions detected in the decimation procedure. The probability $`s`$ for a node to belong to the stable core can be estimated by using Eq. (2) in

$$s\left(t+1\right)=\underset{k=0}{\overset{K}{}}s\left(t\right)^{Kk}\left(1s\left(t\right)\right)^k\left(\genfrac{}{}{0pt}{}{K}{k}\right)p_k,$$
(2)
where $`p_k`$ is the probability that a transition function for given values of $`Kk`$ of its input variables is independent of its other $`k`$ inputs. This equation describes the growth of the stable core with the time. At $`t=0`$ only nodes which happen to have a constant transition function are stable. At later times non-constant transition functions, which receive inputs from stable nodes, can acquire a constant value. In  Eq. (2) was used as a self-consistency equation for infinite systems, i.e. letting $`t\mathrm{}`$. In a finite system, the iteration has to stop at some time $`T`$, which reflects a characteristic length in the network, the maximal distance a signal can flow before it reaches all nodes. The length scale is set by the average distance (in number of links) a signal can travel. The signal pathway in a sparse directed random graph with only a few loops is approximately a branched polymer, where it is known (see e.g. ) that the average distance grows polynomial, i.e. $`TcN^\gamma `$. We have fitted the constants $`c`$ and $`\gamma `$ numerically to our data and find $`\gamma =0.32\left(3\right)`$.

After removing the $`s\left(T\right)`$ stable nodes, the decimation procedure also eliminates the leaves in the interaction, i.e. those non-constant vertices with out-degree $`q\left(t=0\right)=0`$. This changes the out-degree of the remaining nodes. Therefore, this procedure is repeated until no more nodes with $`q\left(t\right)=0`$ are found. The fraction $`P_l`$ of (indirect) leaves can be estimated by the self consistent equation

$$P_l=\underset{q=1}{\overset{\mathrm{}}{}}P\left(q|\stackrel{~}{N},\stackrel{~}{K}\right)P_l^q,$$
(3)
where $`\stackrel{~}{N}=N\left(1s\right)`$ is the number of nodes after removing the $`\eta `$ constant nodes, $`\stackrel{~}{K}`$ is the average in-degree, an $`P\left(q|\stackrel{~}{N},\stackrel{~}{K}\right)`$ the distribution of the out-degree $`q`$ given in Appendix B. We solved Eq. (2) and (3) numerically, the resulting graph is also shown in Fig. 2, which is in very good agreement with our numerical results.

One of the important features of Kauffman’s model is the intrinsic stability of critical Boolean networks. How does decimation affect this behavior? While the network reduction does not change the number and the length of limit cycles, the size of the basins of attraction has to be reduced because the state space is shrunken by orders of magnitude. To get a quantitative picture, we analyze the network stability with respect to the inversion of one randomly chosen variable, after the state trajectory has reached a limit cycle. The probability to end up in a different limit cycle compared to the undisturbed system is shown in Fig. 3. If a limit cycle has not been found within $`10^5`$ steps the network is discarded. For the full network we observe the well known stability, the probability to end up in a a different limit cycle asymptotically approaches zero for large lattices. By contrast, for the decimated network the error-probability grows slowly with the network size and the stability is essentially lost. This means, the tolerance against perturbations observed for Kauffman networks is mostly generated by the stable core: in most cases the perturbed signal is “lost” in the stable core and the full network remains unaffected. It has recently been argued  that the in-homogeneous, for example scale-free, geometry of real world networks is underlying the stability of these systems. Here we find the opposite: stability is primarily generated dynamically by the propagation (Eq. 2) of the stable core in the network logic. The homogeneous geometry plays only a secondary role: if just a geometric reduction of the network is performed, i.e. the leaf-variables are removed (see the discussion of of Eq. (3)), the error-sensitivity is almost unchanged compared to the full network.

The decimation of constant variables from the network enables us to probe a much larger fraction of the state space for a given network. Therefore one may expect to get an better estimate for the number of limit cycles $`n_c`$, which with the commonly used method tends to be underestimated, because some limit cycles may have been missed due to the huge state space. By decimating the networks we can fully enumerate the state space for $`N32`$ and hereby get an unbiased estimate. For larger systems we use the standard method with $`1000`$ restarts on each of the reduced networks. Not unexpectedly we observe a small discontinuity in the curve at the point were the simulation scheme is changed. In Fig. 4 we plot $`n_c`$ as a function of the network size $`N`$. We do not find the quite often cited $`\sqrt{N}`$ behavior for this observable. Rather, we find a linear growth with $`N`$. A possible explanation for the different results obtained in some earlier works may be the bias introduced by the standard method in combination with lower computational power.

## 5 Summary

The source of the remarkable error tolerance of critical Kauffman model is identified as the “dynamics of the stable core”. While this seems to be a contradiction in terms, it quite nicely describes the percolation-like process, which underlies the propagation of the “stability” signal. Starting from the relatively few nodes with transition functions which do not at all depend on their inputs, the islands of frozen states grow in time by the interaction with the already stable nodes. This process is only limited by the finite size of the system. A small fluctuation in the state of the system will most probably not propagate through the stable core and therefore in most cases has no effect. We demonstrate this by reducing given random networks to the dynamical core, where most of the stable, irrelevant variables have been removed. The stability against small fluctuations for these networks is reduced by orders of magnitude and will probably go to zero for infinite networks. It is interesting to observe that these effects are mostly driven by the network logic and not by the network geometry.

For the identification of the dynamical core we have developed a decimation procedure, which is based on inspection of the networks connectivity and logic. The relatively simple procedure works surprisingly well. The results for the size of the stable core are in very good agreement with the values obtained by observing the dynamics of state-space trajectories in the full network .

As a by-product we use the reduced networks to get an improved estimate for the number of limit cycles as a function of the network size. We find that the number of limit cycles grows linearly with $`N`$, which is in sharp contrast to the square-root behavior reported by other groups. Even though this $`\sqrt{N}`$ behavior was an interesting analogy with multi-cellular organisms (with approximately $`\sqrt{N}`$ different cell types for genomes with genome size $`N`$), our result does in no way reduce the importance of Kauffman networks as an example of self-organized order.

Acknowledgments: We have benefitted from discussions with C. Peterson. This work was in part supported by the Swedish Foundation for Strategic Research and the Knut and Alice Wallenberg Foundation through the SWEGENE consortium.

## Appendix A Pseudo Code of the Decimation Procedure

We remind ourselves of the definition of the Kauffman model. $`\sigma _i\left(t\right)`$ denotes the state of variable $`i`$ at time $`t`$ which is determined by $`\sigma _i\left(t\right)=f_i[\sigma _{v_i^1}\left(t1\right),\mathrm{},\sigma _{v_i^K}\left(t1\right)]`$ with $`\left\{v\right\}_i=\{v_i^1,\mathrm{},v_i^K\}`$ defining the connectivity of the network.

We represent the network with a number of nodes, each containing an updating function $`f_i`$ and the list of nodes sending it signals $`\left\{v\right\}_i`$. The function decimate reduces the network $`Net`$ to its relevant core by first removing the stable variables and then those variables that do not regulate any variable.

$`node_i=[f_i,\left\{v\right\}_i]`$  $`Net=\{node_1,\mathrm{},node_N\}`$  decimate($`\mathrm{𝑁𝑒𝑡}`$)  $`\mathrm{𝐝𝐨}`$    $`\mathrm{𝑜𝑙𝑑𝑁𝑒𝑡}`$ := $`\mathrm{𝑁𝑒𝑡}`$   $`\mathrm{𝐟𝐨𝐫}`$ each $`n`$ $``$ $`\mathrm{𝑁𝑒𝑡}`$ $`\mathrm{𝐝𝐨}`$    remove\_unused\_inputs($`n`$, $`\mathrm{𝑁𝑒𝑡}`$)    remove\_if\_no\_inputs($`n`$, $`\mathrm{𝑁𝑒𝑡}`$)  $`\mathrm{𝐰𝐡𝐢𝐥𝐞}`$ $`\mathrm{𝑜𝑙𝑑𝑁𝑒𝑡}`$ $``$ $`\mathrm{𝑁𝑒𝑡}`$  $`\mathrm{𝐝𝐨}`$    $`\mathrm{𝑜𝑙𝑑𝑁𝑒𝑡}`$ := $`\mathrm{𝑁𝑒𝑡}`$   $`\mathrm{𝐟𝐨𝐫}`$ each $`n`$ $``$ $`\mathrm{𝑁𝑒𝑡}`$ $`\mathrm{𝐝𝐨}`$    remove\_if\_no\_outputs($`n`$, $`\mathrm{𝑁𝑒𝑡}`$)  $`\mathrm{𝐰𝐡𝐢𝐥𝐞}`$ $`\mathrm{𝑜𝑙𝑑𝑁𝑒𝑡}`$ $``$ $`\mathrm{𝑁𝑒𝑡}`$  remove\_unused\_inputs($`n`$, $`\mathrm{𝑁𝑒𝑡}`$)  $`\mathrm{𝐟𝐨𝐫}`$ each $`i`$ $``$ $`\{v\}_n`$ $`\mathrm{𝐝𝐨}`$   $`\mathrm{𝑓𝑙𝑎𝑔}`$ := $`\mathrm{𝗍𝗋𝗎𝖾}`$   $`\mathrm{𝐟𝐨𝐫}`$ each configuration of inputs $`\mathrm{𝐝𝐨}`$    $`\mathrm{𝐢𝐟}`$ flipping $`\sigma _i`$ changes outcome of $`f_n`$ $`\mathrm{𝐭𝐡𝐞𝐧}`$     $`\mathrm{𝑓𝑙𝑎𝑔}`$ := $`\mathrm{𝖿𝖺𝗅𝗌𝖾}`$   $`\mathrm{𝐢𝐟}`$ $`\mathrm{𝑓𝑙𝑎𝑔}`$=$`\mathrm{𝗍𝗋𝗎𝖾}`$ $`\mathrm{𝐭𝐡𝐞𝐧}`$    replace $`\sigma _i`$ with constant value $`\mathrm{𝗍𝗋𝗎𝖾}`$ in $`f_n`$    remove $`i`$ from $`\{v\}_n`$  remove\_if\_no\_inputs($`n`$, $`\mathrm{𝑁𝑒𝑡}`$)  $`\mathrm{𝐢𝐟}`$ $`\{v\}_n`$=$`\mathrm{}`$ $`\mathrm{𝐭𝐡𝐞𝐧}`$   $`\mathrm{𝐟𝐨𝐫}`$ each $`m`$ with $`n`$ $``$ $`\{v\}_m`$ $`\mathrm{𝐝𝐨}`$    replace $`\sigma _n`$ with constant value $`f_n\left(\mathrm{}\right)`$ in $`f_m`$    remove $`n`$ from $`\{v\}_m`$  remove\_if\_no\_outputs($`n`$, $`\mathrm{𝑁𝑒𝑡}`$)  $`\mathrm{𝐢𝐟}`$ $`n`$ $``$ $`\{v\}_m`$ for every $`m`$ $`\mathrm{𝐭𝐡𝐞𝐧}`$   remove $`n`$ from $`\mathrm{𝑁𝑒𝑡}`$

## Appendix B In- and Out-degree distribution

The reduced number $`\stackrel{~}{K}`$ of inputs after the decimation described in Eq. (2) is the expectation value of the in-degree for the number of inputs from a non-stable variable:

$$\stackrel{~}{K}=\frac{1\left(1s\left(T\right)\right)+2\left(1s\left(T\right)\right)^2}{1s\left(T\right)}=1+2\left(1s\left(T\right)\right).$$
(B.4)
The out-degree distribution for a node in the random network can be understood by enumerating the number of ways the $`NK`$ links can be distributed over this node and the $`N1`$ remaining nodes, weighted by the corresponding probabilities to choose the nodes:

$$P\left(q|N,K\right)=\left(\frac{1}{N}\right)^q\left(\frac{N1}{N}\right)^{NKq}\left(\genfrac{}{}{0pt}{}{NK}{q}\right).$$
(B.5)


text:
# On the local unitary equivalence of states of multi-partite systems

## I Introduction

Quantum entanglement is the key ingredient in quantum information processing . Bipartite entanglement is quite “well” understood by now, but the investigation of the multi-partite case is just beginning.

Several quantum information protocols such as quantum teleportation and quantum dense coding are based on the equivalence of some classes of states of bi-partite systems under the action of local (one-particle) unitary operations . So the question of local unitary equivalence of multi-partite states of quantum systems is of importance. Recently several authors  have studied the so-called polynomial invariants of local unitaries. In the present paper we discuss the problem of local unitary equivalence. We present a set of necessary conditions which the two multi-partite states have to fulfill in order to declare them locally unitary equivalent.

The state of a quantum system is usually represented by a density matrix. We know that the elements of the density matrix depend on the choice of the basis in Hilbert space, $``$, of the system. Physically the basis can be thought of as a complete set of observables. For example, in the case of spin-1/2 particle, the choice of the basis in Hilbert space corresponds to the direction of the Stern-Gerlach experimental setup .

A local unitary transformation corresponds to a change of a basis in each of the subsystems. Such unitary transformations are known as passive transformations, because there is no physical action performed on the quantum system itself. These transformations simply reflect the choice of our point of view rather than any specific manipulation of the physical system. It reflects the fact that two locally-unitary-equivalent states have the same matrix form, only the choice basis of subsystems is different. An active unitary transformation corresponds to the most general dynamical map of the states for isolated systems. It means that under the dynamics described by local unitaries each of the subsystems evolves independently. Physically it corresponds to the absence of interactions between the subsystems. This property allows us to prepare from a given state of the multi-partite system some other state only by using local unitaries (here we assume that the subsystems can be separated by space-like intervals). From a mathematical point of view there is no difference between the active and passive unitary transformations.

In this paper we address the question: “Under what conditions are the two states states, $`\varrho `$ and $`\sigma `$, of a multi-partite system locally unitary equivalent?” We present a set of conditions which have to be satisfied in order that the two states are locally unitary equivalent. In addition we study whether it is possible to prepare a state of a multi-qudit system which is divided into two parts $`A`$ and $`B`$ by unitaries acting only on the systems $`A`$ and $`B`$, separately. In Section II we present examples of some invariants together with some applications for pure states. A criterion about the local unitary equivalence for mixed states is given in Section III. In Section IV we investigate the connection of bipartite and multi-partite entanglement via local unitaries for system of $`N`$ qudits.

## II Local unitary equivalence

Let us consider two states $`\varrho `$ and $`\sigma `$ of an $`N`$party system described by the Hilbert space $`_{A_1\mathrm{}A_N}`$. The question of local unitary equivalence can be mathematically formulated in the following way: Does there exist a unitary transformation $`U=U_1\mathrm{}U_N`$, such that $`\varrho =U\sigma U^{}`$?

In order to study this problem of local unitary equivalence we firstly, simplify our task. Specifically, we divide the total Hilbert space in the following way

$`_{A_1\mathrm{}A_N}=_{A_1}\mathrm{}_{A_N}=_X_Y.`$ (1)
where $`X`$ and $`Y`$ are two sets of subsystems $`A_k`$. If we would have some criteria for bi-partite systems being locally unitary equivalent, then we can decide whether the following is true

$`\varrho =U_XU_Y\sigma U_X^{}U_Y^{}.`$ (2)
If the answer is “no”, then the states are not locally unitary equivalent. If the answer is “yes”, then we continue our investigation by checking whether the pair of states

$`\varrho _X=\mathrm{Tr}_Y\varrho `$ $`\mathrm{and}`$ $`\sigma _X=\mathrm{Tr}_Y\sigma `$ (3)
$`\varrho _Y=\mathrm{Tr}_X\varrho `$ $`\mathrm{and}`$ $`\sigma _Y=\mathrm{Tr}_X\sigma `$ (4)
are locally unitary equivalent. That is, we again, divide each of the subsystems $`X,Y`$ into two smaller sets and we use the criteria for bi-partite system. Using this method recursively we finally partition the multi-partite system into its components $`A_k`$ and give an answer to the question of locally unitary equivalence. From here it follows that it is enough to have a criterion for the bi-partite composite systems with the help of which multi-partite systems can be studied. In what follows , we present some necessary conditions for states to be locally unitary equivalent.

### A Invariants of local unitaries

1. Eigenvalues. It is well known that the necessary and sufficient condition for two states to be unitarily equivalent is that they have the same set of eigenvalues. Let us denote by $`eig(\varrho )`$ the set of eigenvalues of the density matrix $`\varrho `$. Then

$`\varrho ^{}=U\varrho U^{},\mathrm{iff}eig(\varrho ^{})eig(\varrho ).`$ (5)
2. Eigenvalues of the reduced states. The reduced density matrices $`\varrho _X`$ and $`\varrho _Y`$ change after the action of the local unitary transformation $`U_XU_Y`$ as follows: for all $`U_X`$ and $`U_Y`$

$`\varrho _X^{}`$ $`=`$ $`\mathrm{Tr}_Y\left[U_XU_Y\varrho U_X^{}U_Y^{}\right]=U_X\varrho _XU_X^{},`$ (6)
$`\varrho _Y^{}`$ $`=`$ $`\mathrm{Tr}_X\left[U_XU_Y\varrho U_X^{}U_Y^{}\right]=U_Y\varrho _YU_Y^{}.`$ (7)
Now we can define the first invariant: local unitary transformations do not change the eigenvalues of the reduced density matrices, which describe states of the subsystems.

For any state $`\varrho ^{}`$ such that $`\varrho ^{}=U_XU_Y\varrho U_X^{}U_Y^{}`$ it is valid that

$`eig(\varrho _X)`$ $``$ $`eig(\varrho _X^{}),`$ (8)
$`eig(\varrho _Y)`$ $``$ $`eig(\varrho _Y^{}).`$ (9)
We note that this relation is valid for any operator $`\varrho `$, not only for density operators.
3. Eigenvalues of the reduced eigenprojectors. Quantum states can be expressed in a a spectral form. Specifically, let $`U=U_XU_Y`$ and $`\varrho =_k\lambda _kP_k`$, then $`\varrho ^{}=U\varrho U^{}=_k\lambda _kP_k^{}`$, where $`P_k`$ and $`P_k^{}=UP_kU^{}`$ are eigenprojectors belonging to the eigenvalue $`\lambda _k`$. Since the unitary transformation under consideration is local, i.e. $`U=U_XU_Y`$, we obtain the implication

$`\varrho ^{}=U\varrho U^{}`$ $`\mathrm{then}`$ $`dim(P_k)=dim(P_k^{})`$ (12)
$`eig(\mathrm{Tr}_YP_k)eig(\mathrm{Tr}_YP_k^{})`$
$`eig(\mathrm{Tr}_XP_k)eig(\mathrm{Tr}_XP_k^{}),`$
where $`dim(P_k)`$ denotes the dimension of the projection $`P_k`$ defined as the number of its nonzero eigenvalues.
4. Partial transposition. We use the following identity for any operators $`K_X,L_Y,M_X,N_Y`$

$`(K_XL_Y\varrho M_XN_Y)^{T_Y}=K_XN_Y^T\varrho ^{T_Y}M_XL_Y^T`$ (13)
where by $`T_Y`$ we denote the partial transposition of the subsystem $`Y`$. A similar identity is fulfilled also for the partial transposition, $`T_X`$, of the system $`X`$. Note, that the partial transposition depends on the choice of the local basis. Introducing the local unitary transformation $`U_XU_Y`$ we obtain

$`(\varrho ^{})^{T_Y}`$ $`=`$ $`(U_XU_Y\varrho U_X^{}U_Y^{})^{T_Y}`$ (14)
$`=`$ $`U_X(U_Y^{})^T\varrho ^{T_Y}U_X^{}U_Y^T`$ (15)
$`=`$ $`W\varrho ^{T_Y}W^{}`$ (16)
where $`W=U_X(U_Y^T)^{}`$ is a unitary operator. Its unitarity follows from the fact that the transposed unitary operator is a unitary operator. The unitary transformation leaves the eigenvalues of the operator $`\varrho ^{T_Y}`$ unchanged. Therefore, we obtain the second invariant associated with density matrices which are invariant under local unitary transformations. Namely, if $`\varrho ^{}=U_XU_Y\varrho U_X^{}U_Y^{}`$ then

$`eig(\varrho ^{T_Y})`$ $``$ $`eig(\varrho ^{T_Y})`$ (17)
$`eig(\varrho ^{T_X})`$ $``$ $`eig(\varrho ^{T_X}).`$ (18)

### B Local unitary equivalence for pure states

Let us consider two pure states $`|\psi `$ and $`|\varphi `$ of a bi-partite system written in the Schmidt-decomposed form

$`|\psi `$ $`=`$ $`{\displaystyle \underset{k}{}}\sqrt{\lambda _k}|k_X|k_Y`$ (19)
$`|\varphi `$ $`=`$ $`{\displaystyle \underset{k}{}}\sqrt{\mu _k}|k^{}_X|k^{}_Y.`$ (20)
From the equality $`\lambda _k=\mu _k`$ (for all $`k`$) the validity of the second invariant property follows. That is the reduced states have the same set of eigenvalues, namely $`\{\lambda _k\}`$. The first invariant property is fulfilled trivially for pure states. Moreover, let us define local unitaries by the equations $`U_X|k_X=|k^{}_X`$ and $`U_Y|k_Y=|k^{}_Y`$. Then the states (20) are locally unitary equivalent and $`|\varphi =U_XU_Y|\psi `$.

Theorem 1
Let us consider two pure states $`|\psi `$ and $`|\varphi `$ of a bi-partite system. These state are locally unitary equivalent if and only if their coefficients in the Schmidt decomposition are equal, i.e. the sets of eigenvalues of reduced states coincide.

Example 1.
This simple theorem can be rather useful. In particular, we can use it to prove that any GHZ state is locally equivalent to one EPR pair and some ancilla system in an arbitrary state. It means we are able to prepare a GHZ state from an EPR pair a vice versa . Let us define two states of $`N`$ qubits $`|GHZ_{1\mathrm{}N}=(|0^{(N)}+|1^{(N)})/\sqrt{2}`$ and $`|\mathrm{\Psi }_{1\mathrm{}N}=|EPR_{12}|\varphi _{3\mathrm{}N}`$, where $`|EPR=(|00+|11)/\sqrt{2}`$. Since $`\mathrm{Tr}_{2\mathrm{}N}|GHZGHZ|=\mathrm{Tr}_{2\mathrm{}N}|\mathrm{\Psi }\mathrm{\Psi }|=\frac{1}{2}𝟙`$ we know that the eigenvalues of the reduced states coincide. Specifically, let us consider a splitting of the whole system into a single qubit denoted as “1” and the rest of the system. It means that these states are locally (with respect to the given splitting) unitary equivalent.

Example 2.
Let us consider a state (for more details see )

$`|X={\displaystyle \frac{1}{\sqrt{2^n}}}{\displaystyle \underset{k=0}{\overset{2^n1}{}}}|k_X|\alpha _k_Y`$ (21)
of $`N=2n+1`$ qubits located on a ring. This system is divided into two subsystem of $`n`$ neighboring qubits (subsystem $`X`$) and $`(n+1)`$ neighboring qubits (subsystem $`Y`$). In this notation $`|k_X`$ represents the state of $`n`$ qubits, if $`k`$ is written in a binary form. The state $`|\alpha _k_Y=_l(1)^{d_k+d_l+d_{k,l}}|l_Y`$ where $`|l_Y`$ is the state of the rest of the $`n+1`$ qubits, if $`l`$ is written in a binary form. The number $`k`$ and $`l`$ define the state of the ring in the sense of its binary expression. From here we can determine the state of each individual qubit. The numbers $`k`$ give us the qubits in the state $`|1`$ and we can write $`\{i,\mathrm{},j\}`$ to represents the position of the qubits in this state. Similarly, $`l`$ determines the set $`\{r,\mathrm{},s\}`$. We define $`d_k=d_{\{i,\mathrm{},j\}}:=_{i<j}d(i,j)`$, $`d_{k,l}=_{a\{i,\mathrm{}j\},b\{r,\mathrm{},s\}}d(a,b)`$, where $`d(i,j)`$ is the shortest distance between two qubits (in the state $`|1`$) on a ring (i.e. number of steps needed to get from qubit in position $`i`$ into the position $`j`$ around the ring, or the smallest number of qubits lying between the $`i`$-th and $`j`$-th position). The state (21) is one of the eigenstates of the Ising model (for more details see ).

For the purpose of our example it is important that Eq. (21) represents the Schmidt form of the state $`|X`$ with all Schmidt coefficients equal to $`1/\sqrt{2^n}`$ and the Schmidt number (i.e. number of nonzero Schmidt coefficients) is equal to $`2^n`$. Let us define a state

$`|\mathrm{\Psi }=|\varphi _{Y_{n+1}}\left({\displaystyle \underset{j=1}{\overset{n}{}}}|EPR_{X_jY_j}\right),`$ (22)
where $`|EPR_{X_jY_j}=(|0_{X_j}|0_{Y_j}+|1_{X_j}|1_{Y_j})/\sqrt{2}`$ is the EPR pair shared between the subsystems $`X`$ and $`Y`$, such that $`j`$ represents the position of the qubit in this subsystems. We can express the state (22) in the Schmidt form

$`|\mathrm{\Psi }={\displaystyle \frac{1}{\sqrt{2^n}}}{\displaystyle \underset{k=0}{\overset{2^n1}{}}}|k_X|\beta _k_Y,`$ (23)
where $`|\beta _k_Y=|k_Y^{}|\varphi _{n+1}`$ and $`Y^{}`$ denotes the subsystem $`Y`$ without one of its qubit at the $`(n+1)`$-th position. Since the set of Schmidt coefficients for these states coincides, we conclude that the states (21) and (22) are locally unitarily equivalent. We stress that this equivalence is with respect to the splitting of the ring into two parts $`X`$ and $`Y`$. Therefore, from the point of view of entanglement between systems $`X`$ and $`Y`$ these states are equivalent, and we are able to prepare $`n`$ EPR pairs between these systems by a local transformation $`U=𝟙_𝕏𝕌_𝕐`$, where $`U_Y`$ is determined by the equality $`U_Y|\alpha _k_Y=|\beta _k_Y`$. Finally we note, that this transformation is not uniquely defined, because of the degeneracy of Schmidt coefficients.

### C Counterexample for mixed states

Unfortunately, Theorem 1 is valid only for pure states and cannot be generalized for arbitrary impure states. In what follows we will show a simple counterexample. The idea behind this example is that we cannot create entanglement by a local unitary transformation and that the spectral decomposition of a density matrix is unique. It means the separable pure state in the spectral decomposition cannot evolve into entangled one.

Example 3.
Let us consider two states $`\varrho _{XY}`$ and $`\sigma _{XY}`$ of a bi-partite system $`_𝒳_𝒴`$, where $`dim_Y=2`$ and $`dim_Y=8`$. This means the system of four qubits is divided into a one-qubit and a three-qubit subsystems. Let the set of eigenvalues be the same

$`eig(\varrho _{XY})=eig(\sigma _{XY})=\{{\displaystyle \frac{1}{4}},{\displaystyle \frac{3}{8}},{\displaystyle \frac{5}{16}},{\displaystyle \frac{1}{16}}\}`$ (24)
and the corresponding eigenvectors be

$`\left(\begin{array}{c}\frac{1}{\sqrt{4}}\left(|0,1+|0,2+|0,4+|1,0\right)\\ |1,7\\ |0,5\\ \frac{1}{\sqrt{3}}(|1,1+|1,2+|1,4)\end{array}\right)`$ (29)
for $`\varrho _{XY}`$ and

$`\left(\begin{array}{c}|1,7\\ \frac{1}{\sqrt{3}}(|0,7+|1,3+|0,5)\\ |0,0\\ |0,1\end{array}\right)`$ (34)
for $`\sigma _{XY}`$. It is easy to check that these two states also have the same sets of eigenvalues of their reduced density matrices, i.e.

$`eig(\varrho _X)`$ $`=`$ $`eig(\sigma _X)=\{{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}}\}`$ (35)
$`eig(\varrho _Y)`$ $`=`$ $`eig(\sigma _Y)=\{{\displaystyle \frac{3}{8}},{\displaystyle \frac{5}{16}},{\displaystyle \frac{1}{4}},{\displaystyle \frac{1}{16}}\}`$ (36)
So for the state under considerations the conditions of Theorem 1 are satisfied. On the other hand, in contradiction to the general properties of local unitaries, we see the violation of the property the entanglement creation prohibited by local unitary transformation. In other words, the condition of the same set of reduced eigenvalues of eigenprojectors (invariant A3) is not fulfilled.

## III Local unitary equivalence for the mixed states

We have shown that Theorem 1 is not valid for all states. Certainly, it can be trivially generalized to a class of the so called factorizable states, given by density operators $`\varrho _{XY}=\varrho _X\varrho _Y`$. On the other hand, it is not clear whether the invariants presented in Section II provide us with a definite answer for the problem of the local unitary equivalence for mixed states. To clarify the situation we formulate the following theorem

Theorem 2
Suppose two non-degenerate states $`\varrho `$ and $`\sigma `$ of the composite system $`X+Y`$ are given, i.e.

$`\varrho ={\displaystyle \underset{k}{}}\lambda _k|\psi _k\psi _k|,\sigma ={\displaystyle \underset{k}{}}\mu _k|\varphi _k\varphi _k|.`$ (37)
and $`\lambda _k\lambda _l,\mu _k\mu _l`$ for $`kl`$. Let us express the vectors $`|\psi _1`$ and $`|\varphi _1`$ in their Schmidt bases and fix these two basis on $`_X_Y`$. If for each $`k`$, $`\lambda _k=\mu _k`$ and the coefficients of the corresponding eigenvectors $`|\psi _k=_{mn}\alpha _{mn}^k|mn`$ and $`|\varphi _k=_{mn}\beta _{mn}^k|m^{}n^{}`$ (where $`|mn`$ and $`|m^{}n^{}`$ are the previously mentioned fixed bases) coincide, i.e. $`\alpha _{mn}^k=\beta _{mn}^k`$, then the states $`\varrho `$, $`\sigma `$ are locally unitary equivalent. The local transformation is given by $`U_X|m=|m^{}`$ and $`U_Y|n=|n^{}`$.

From the construction it is clear that the theorem is valid. Note that the choice of the fixed basis is given by the Schmidt expression of one of the eigenvectors, but the result doesn’t depend on it. In the degenerate case we cannot use the argument with the Schmidt decomposition, because then the eigenprojectors are no longer one dimensional.

We conclude this section by saying that in general by using invariants of local unitaries we are able only to determine whether the states under consideration are not locally unitary equivalent.

## IV Local connection of multi-partite and bipartite entanglement

In this section we generalize Example 2 of Section III. Specifically, let us consider a pure state $`|\psi _{A_1,\mathrm{},A_N}`$ of a composite system consisting of $`N`$ qudits. Let us split this system into two multi-partite subsystems $`X`$ and $`Y`$. By calculating the entropy of one of the subsystems we can determine the degree of entanglement between these two parties. An interesting question is whether with the help of a local transformation with respect to the splitting we are able to prepare a state where the entanglement is shared only between the specific pairs of qudits such that one particle of each pair belongs to the system $`X`$ while the second belongs to the system $`Y`$. In other words, the qudits belonging to the subsystem $`X`$ ($`Y`$) are mutually separable. It means that the multi-partite entanglement shared between the subsystems $`X`$ and $`Y`$ can be arranged by bi-partite entanglement shared between the pairs of qudits just by performing an appropriate local unitary transformation $`U=U_XU_Y`$ (see Fig. 1).

Let us denote by $`A_1,\mathrm{},A_N`$ the qudits of the whole system and by $`X_1,\mathrm{},X_n`$ ($`Y_1,\mathrm{},Y_m`$) the qudits in a subsystem $`X`$ ($`Y`$). Note that $`n+m=N`$. We also introduce a function $`E_{XY}(|\psi _{XY})`$ which describes the degree entanglement shared between the systems $`X`$ and $`Y`$. Suppose $`nm`$. The most general pure state of the whole system, for which the entanglement between the subsystems $`X`$ and $`Y`$ is equal to the entanglement shared by pairs of qudits is given by the expression

$`|\mathrm{\Phi }_{XY}=\left({\displaystyle \underset{j=1}{\overset{n}{}}}|\chi _j\right)|rest_{Y_{n+1},\mathrm{}Y_m}`$ (38)
where $`|\chi _j=_k\sqrt{a_{k_j}}|\chi _k_{X_j}|\chi _k_{Y_j}`$ is state written in Schmidt basis. The entanglement between the systems $`X`$ and $`Y`$ which are in the state (38) is given can be represented by the expression

$`E_{XY}(|\mathrm{\Phi }_{XY})={\displaystyle \underset{jk}{}}a_{k_j}\mathrm{log}a_{k_j}.`$ (39)
The equality $`E_{XY}(|\psi _{XY})=E_{XY}(|\mathrm{\Phi }_{XY})`$ is not sufficient to determine whether the states $`|\psi _{XY}`$ and $`|\mathrm{\Phi }_{XY}`$ are locally unitary equivalent. Let us rewrite both states in their Schmidt bases

$`|\psi _{XY}`$ $`=`$ $`{\displaystyle \underset{r=0}{\overset{d^n1}{}}}\sqrt{q_r}|\alpha _r_X|\beta _r_Y,`$ (40)
$`|\mathrm{\Phi }_{XY}`$ $`=`$ $`{\displaystyle \underset{r=0}{\overset{d^n1}{}}}\sqrt{p_r}|r_X\left(|r_{Y_1,\mathrm{}Y_n}|rest_{Y_{n+1},\mathrm{},Y_m}\right),`$ (41)
where $`d=dim(_{A_j})`$, $`r`$ is a base d number which represents the state of $`n`$ qudits, e.g. $`r=4`$ in a binary form represents a state of 5 qubits $`|0,0,1,0,0`$, or in the base 3 form it represents a state of 5 qutrits $`|0,0,0,1,2`$.

The eigenvalues of the reduced states are $`\{q_r\}`$ for the state $`|\psi _{XY}`$ and $`\{p_r=a_{k_1}\mathrm{}a_{k_n}\}`$ for the state $`|\mathrm{\Phi }_{XY}`$, where $`r=(k_1,\mathrm{},k_n)`$ in the number base-$`d`$ form, i.e. $`k_j=0,1,\mathrm{},d1`$. Theorem 1 implies that these sets of eigenvalues have to be identical in order to ensure local unitary equivalence. Our aim is to find some condition under which the equations

$`q_r=p_r=p_{(k_1,\mathrm{},k_n)}=a_{k_1}\mathrm{}a_{k_n}.`$ (42)
are satisfied. In what follows we have $`d^n`$ equations with $`nd`$ undetermined variables. These equation, in general, do not have a solution. In that case the states are not locally unitary equivalent and the entanglement shared by systems $`X`$ and $`Y`$ cannot be prepared by local unitaries from bi-partite entanglement shared between the subsystems. If the system of equations (42) has a solution then the local unitary transformations are given by relations

$`U_X|\alpha _r_X`$ $`:=`$ $`|\alpha _r_X`$ (43)
$`U_Y|\beta _r_Y`$ $`:=`$ $`|r_{Y_1,\mathrm{},Y_n}|rest_{Y_{n+1},\mathrm{},Y_m}.`$ (44)
We note, that from above, the local unitary equivalence between the $`|X`$ state and $`n`$ EPR pairs shared between systems $`X`$ and $`Y`$ follows (see Section II).

## V Conclusion

In this paper we have analyzed the problem of finding under which conditions two states states $`\varrho `$ and $`\sigma `$ of a multi-partite system are locally unitary equivalent. We have presented a set of conditions which have to be satisfied in order that the two states are locally unitary equivalent. We have presented a set of parameters which play the role of invariants under the action of local unitary transformations. In addition, we have shown how these invariants can be used to characterize the local unitary equivalence of states.

We have shown that the problem of local unitary equivalence for multi-partite systems can be reduced to the problem of analysis of bi-partite systems. Specifically, it is enough to investigate the case of bi-partite local unitary equivalence and to show sufficient conditions of this property for non-degenerate density matrices.

We left two questions opened: the local unitary equivalence for degenerate density matrices and the characterization of the connection between multi-partite and bipartite entanglement, since, as it is well known that there exists an intrinsic multi-partite entanglement .

Finally, we have analyzed multi-qudit systems. We have shown, that when analyzing local unitary equivalence of two states it is useful to consider the splitting of the whole system into two subsystems. $`X`$ and $`Y`$. In this case a system of equations of the form (42) results. We have shown that these equations might not have a solution from which it directly follows that the state under consideration cannot be prepared by local unitary transformations from any number of partially entangled pairs shared between two systems $`X`$ and $`Y`$ (i.e. these states cannot be prepared just using a bi-partite entanglement as a resource and local transformations on subsystems $`X`$ and $`Y`$).

###### Acknowledgements.

We thank Mark Hillery for discussions and for reading the manuscript. V.B. thanks Chris Fuchs for a stimulating discussion. This work was supported by the IST project EQUIP under the contract IST-1999-11053 and by the grant GACR 201/01/0413.


text:
# Mesoscopic thermal transport through a weak link

## I introduction

Recently there has been considerable interest in the physics of phonons in mesoscopic and nanoscale systems. In a beautiful experiment, Schwab et al.  observed a quantization of the thermal conductance in freely suspended one-dimensional dielectric wires, corresponding to a conductance given by $`\pi k_\mathrm{B}^2T/6\mathrm{}`$ per transmitted vibrational mode, analogous to the well-known electrical conductance quantization in units of $`e^2/2\pi \mathrm{}`$ per spin-resolved channel in one-dimensional mesoscopic conductors . The thermal conductance quantization can be understood  by using a thermal analog of the Landauer-Büttiker formula .

In this paper we consider mesoscopic thermal transport through a wire or weak link in the limit where the phonon transmission probability through the wire is small (the so-called nonadiabatic regime) and the thermal conductance is much less than $`\pi k_\mathrm{B}^2T/6\mathrm{}`$, the value corresponding to one fully propagating channel. The geometry we consider is shown schematically in Fig. 1.

The mechanical bridge may consist of one or more chemical bonds, or by a narrow wire of dielectric material, both of which we model by a harmonic spring of stiffness $`K`$ . We obtain a general expression for the thermal current that can be regarded as a thermal analog of the well-known formula, derived by Schrieffer et al., for the electrical current through a tunneling barrier. Our result can also be interpreted as an application of the thermal Landauer formula in the weak-tunneling limit, with the energy-dependent phonon transmission probability calculated microscopically.

## II thermal current in the weak-tunneling limit

We consider two macroscopic dielectrics, L and R, held at fixed temperatures $`T_\mathrm{L}`$ and $`T_\mathrm{R}`$. The Hamiltonian of the isolated solids is $`H_0=H_\mathrm{L}+H_\mathrm{R},`$ where

$$H_I\underset{n}{}\omega _{In}a_{In}^{}a_{In},I=\mathrm{L},\mathrm{R}.$$
(1)
The $`a_{nI}^{}`$ and $`a_{nI}`$ are phonon creation and annihilation operators for the left and right sides. The vibrational modes of the isolated bodies are labeled by $`n`$ and have energies $`\omega _{In}`$.

The two bodies are connected by a harmonic spring with stiffness $`K`$,

$$\delta H=\frac{1}{2}K\left(u_\mathrm{L}^zu_\mathrm{R}^z\right)^2.$$
(2)
Here $`u_I^z`$ is the normal component of the displacement field $`𝐮(𝐫)`$ at the surface of body $`I`$ at the point of connection to the weak link, with the local surface normal taken to be in the $`z`$ direction. Our mesoscopic weak-link model is illustrated in Fig. 2.

The system is described by the Hamiltonian $`H=H_0+\delta H`$. The thermal current is calculated by defining a thermal current operator $`\widehat{I}_{\mathrm{th}}`$ according to

$$\widehat{I}_{\mathrm{th}}_tH_\mathrm{R}=i[H,H_\mathrm{R}].$$
(3)
The expectation value of $`\widehat{I}_{\mathrm{th}}`$ is the energy per unit time flowing from the left to the right body. To leading order in $`\delta H`$ (the weak-coupling limit) we obtain

$$I_{\mathrm{th}}=\frac{2\pi K^2}{\mathrm{}}_0^{\mathrm{}}𝑑ϵϵN_\mathrm{L}(ϵ)N_\mathrm{R}(ϵ)\left[n_\mathrm{L}(ϵ)n_\mathrm{R}(ϵ)\right],$$
(4)
an expression analogous to the formula derived by Schrieffer et al.  for the electrical current through a tunneling barrier. Here $`n_\mathrm{L}(ϵ)`$ and $`n_\mathrm{R}(ϵ)`$ are Bose distribution functions with temperatures $`T_\mathrm{L}`$ and $`T_\mathrm{R}`$, and

$$N_I(\omega )\frac{1}{\pi }\mathrm{Im}D_I(\omega )$$
(5)
is the local spectral density, defined in terms of the Fourier transform of the retarded correlation functions

$$D_I(t)i\theta (t)[u_I^z(t),u_I^z(0)]_0$$
(6)
for the isolated macroscopic bodies L and R. The spectral density defined in (5) is different than the ordinary thermodynamic DOS, even in a homogeneous system. The quantity (5), however, is the one relevant here.

Our result (4) shows that the (linear) thermal conductance of a weak mesoscopic link, defined by

$$G_{\mathrm{th}}\underset{T_\mathrm{L}T_\mathrm{R}}{lim}\frac{I_{\mathrm{th}}}{T_\mathrm{L}T_\mathrm{R}},$$
(7)
will vary at low temperatures as a power-law in $`T`$. If we let $`\alpha `$ be the exponent characterizing the power-law phonon spectral density of the reservoirs at low energies, then

$$G_{\mathrm{th}}T^{2\alpha +1}.$$
(8)
For example, $`\alpha =1`$ at the planar surface of a semi-infinite isotropic elastic continuum.

A simple application of our theory to a nanostructure consisting of a cylindrical neck of Si material connecting two semi-infinite Si crystals has been given in Ref. .

## III discussion

The classical theory of thermal conduction, based on the heat equation and on the concept of a local thermal conductivity, is entirely inapplicable to mesoscopic dielectrics. In a mesoscopic system, thermal resistance is caused by elastic scattering of phonons, whereas in an infinite, disorder-free crystal it is caused by anharmonicity. In the non-mesoscopic regime (for example, at higher temperatures), the thermal conductivity of a narrow bridge is determined by the bridge material’s bulk thermal conductivity $`\kappa `$, but in the mesoscopic regime it is determined by the mechanical properties of the bridge material, through the stiffness $`K`$.

## ACKNOWLEDGMENTS

This work was supported by NSF CAREER Grant No. DMR-0093217, and by a Research Innovation Award and Cottrell Scholars Award from the Research Corporation. It is a pleasure to thank Andrew Cleland, Steve Lewis, and Aleksandar Milosevic for useful discussions.


text:
# Rossby-Haurwitz waves of a slowly and differentially rotating fluid shell

## 1 Introduction

The investigation of $`r`$ modes in rotating stars has seen a renewed interest in recent years since they were shown to be unstable to the emission of gravitational waves (Andersson 1998, Friedman and Morsink 1998). In the last three years a number of papers have outlined the basic properties of the $`r`$-mode oscillations, their role in triggering and feeding the instability, and the importance of this instability in the emission of gravitational waves as well as in the slowing down of rapidly rotating, hot neutron stars (see Andersson and Kokkotas 2001 and Friedman and Lockitch 2001 for some recent reviews). Recent studies, on the other hand, have focussed on more intricate but also more realistic aspects of the instability. Such aspects include the interaction of the instability with magnetic fields (Spruit 1999, Rezzolla et al. 2000), nonlinear and secular effects (Rezzolla et al. 2001a, 2001b; Stergioulas and Font 2001; Lindblom et al. 2001), the presence of differential rotation in Newtonian stars (Karino et al. 2001, hereafter KYE), of a continuous frequency spectrum in the slow-rotation approximation of uniformly rotating relativistic stars (Kojima and Hosonuma 1999) and the lack of discrete, physically plausible mode solutions within the continuous part of the spectrum (Ruoff and Kokkotas 2001). We here pay attention to one of these more “subtle” aspects and focus on the role played by differential rotation in Newtonian stars. In particular, we address the question of whether a large degree of differential rotation in slowly rotating stars turns the eigenvalue problem for $`r`$ modes into a singular eigenvalue problem, thus preventing the existence of $`r`$ modes. This is an important question because of the possible observation of a singular problem in $`r`$-mode oscillations that has emerged from recent calculations of differentially rotating axi-symmetric Newtonian stars (KYE). More intriguingly, similar evidence is converging also from different approaches such as those looking at purely axial $`r`$ modes in rotating non-isentropic relativistic stars in the slow-rotation approximation (Yoshida 2001, Yoshida and Futamase 2001; Ruoff and Kokkotas, 2001a, 2001b) although a regularization may be possible in some cases (Lockitch and Andersson 2001).

Part of the subtlety introduced by differential rotation in the eigenvalue problem is due to the fact that, above a certain degree of differential rotation, the eigenvalue equations become very stiff as a result of large radial and polar gradients of the rotational angular frequency. Under these circumstances, the numerical solution requires increasingly high accuracy and computational costs. To circumvent this problem, we have here resorted to a simpler model based on a differentially rotating, thin shell. A thin-shell model has been adopted also by Levin and Ushomirsky (2001) to investigate the effect of electromagnetic radiation-reaction for $`r`$ modes on a uniformly rotating shell. This approach replaces the set of partial differential equations of a multidimensional star model with an ordinary differential equation which can be solved to much higher accuracy and with modest computational costs. Although much simpler to solve, the $`r`$-mode eigenvalue problem for a thin shell incorporates many of the mathematical properties of the corresponding eigenvalue problems for multidimensional Newtonian stars or for slowly-rotating relativistic stars. As will become clear in the following Sections, adopting a thin-shell model has been very valuable for gaining insight into the behaviour of $`r`$ modes in differentially rotating fluids.

## 2 Basic Equations

The thin-shell model adopted here is based on the following simplifying assumptions: (i) the background star is slowly rotating and we will omit terms in the hydrodynamical equations of $`𝒪(\mathrm{\Omega }^2)`$; (ii) the fluid is incompressible and inviscid; (iii) the shell is spherical with radius $`R`$ for every rate of differential rotation; (iv) the fluid motion is constrained on the shell, that is, we have no radial component of velocity.

In the inertial (nonrotating) frame the fluid velocity $`\stackrel{}{v}`$ is the composition of the unperturbed velocity of the differentially rotating shell $`\stackrel{}{v}__0`$, where, in spherical polar coordinates, $`v__0^i=(0,0,\mathrm{\Omega }(\theta )R\mathrm{sin}\theta )`$, with the velocity perturbation $`\delta \stackrel{}{u}`$ which has only tangential components, i.e. $`\delta u^i=(0,u^\theta ,u^\varphi )`$. As a result, the components of the fluid velocity can be written as

$$v^i=(0,u^\theta ,\mathrm{\Omega }R\mathrm{sin}\theta +u^\varphi ),$$
(1)
In the simple model considered here, the fluid motion in the shell is fully described by the Euler and continuity equations. The latter, in particular, is trivially satisfied by the background motion and assumes the following form in its perturbed part

$$_\theta (u^\theta \mathrm{sin}\theta )+_\varphi u^\varphi =0.$$
(2)
As customary with this type of problem, it is here convenient to introduce the fluid vorticity $`\stackrel{}{\xi }\frac{1}{2}\times \stackrel{}{v}`$, whose radial component has the form

$$\xi ^r=\mathrm{\Omega }\mathrm{cos}\theta +\frac{1}{2}\frac{d\mathrm{\Omega }}{d\theta }\mathrm{sin}\theta +\frac{1}{2r\mathrm{sin}\theta }\left[_\theta (u^\varphi \mathrm{sin}\theta )_\varphi u^\theta \right].$$
(3)
It is a simple exercise to show that the Euler equation can then be cast into a conservation equation for the vorticity $`_t\stackrel{}{\xi }+\times (\stackrel{}{\xi }\times \stackrel{}{v})=0`$, whose radial component is

$$\left(_t+\frac{v^\theta }{r}_\theta +\frac{v^\varphi }{r\mathrm{sin}\theta }_\varphi \right)\xi ^r=0.$$
(4)
Because we are interested in the harmonic modes, we assume that the tangential components of the perturbation velocity can be written as $`u^\theta ,u^\varphi \mathrm{exp}(i\sigma t+im\varphi )`$, so that the linearized version of equation (4) is

$$i(\sigma m\mathrm{\Omega })\left[\frac{d}{d\theta }(u^\varphi \mathrm{sin}\theta )imu^\theta \right]+2u^\theta \mathrm{sin}\theta \frac{d}{d\theta }\left(\mathrm{\Omega }\mathrm{cos}\theta +\frac{1}{2}\frac{d\mathrm{\Omega }}{d\theta }\mathrm{sin}\theta \right)=0.$$
(5)
We can now use the continuity equation (2) to rewrite the radial component of the vorticity conservation equation as

$$\frac{d}{d\mu }\left[(1\mu ^2)\frac{d\chi }{d\mu }\right]\frac{m^2}{1\mu ^2}\chi \frac{2m\mathrm{\Omega }}{\sigma m\mathrm{\Omega }}\left[1+\left(\frac{2\mu }{\mathrm{\Omega }}\right)\frac{d\mathrm{\Omega }}{d\mu }\left(\frac{1\mu ^2}{2\mathrm{\Omega }}\right)\frac{d^2\mathrm{\Omega }}{d\mu ^2}\right]\chi =0,$$
(6)
where we have introduced the new coordinate $`\mu \mathrm{cos}\theta `$ and where $`\chi u^\theta \mathrm{sin}\theta `$. Together with regular boundary conditions at the northern and southern poles ($`\mu =\pm 1`$), equation (6) accounts for the eigenvalue problem of the fluid normal modes on the shell. In the case of a uniformly rotating shell, $`d\mathrm{\Omega }/d\mu =0=d^2\mathrm{\Omega }/d\mu ^2`$, and equation (6), first derived by Haurwitz in 1940 (see also Stewartson and Rickard, 1969), reduces to Legendre’s equation whose regular solution is $`\chi =P_^l^_m(\mu )`$, with $`P_^l^_m(\mu )`$ being the associated Legendre functions. Note that in the case of uniform rotation the eigenfrequencies obey the well known dispersion relation for $`r`$ modes

$$\sigma =m\mathrm{\Omega }\frac{2m\mathrm{\Omega }}{l(l+1)},$$
(7)
and this justifies calling these modes Rossby-Haurwitz waves. For simplicity, however, hereafter we will refer to them as $`r`$ modes.

The solution of equation (6) in the case of a differentially rotating shell depends on the law of differential rotation chosen and on its first and second polar derivatives. The choice of a law of differential rotation is, in this sense, somewhat arbitrary but we have here followed previous work on the subject (Eriguchi and Müller 1985, KYE) and modeled the differential rotation through either a “$`j=\mathrm{const}.`$” law

$$\mathrm{\Omega }=\frac{A_j^2+1}{A_j^2+1\mu ^2}\mathrm{\Omega }__\mathrm{E},$$
(8)
or a “$`v=\mathrm{const}.`$” law111These differential laws take their names from the fact that through the variation of the parameter $`A_{j,v}`$ they represent families of rotation laws spanning the range between uniform rotation (for $`A_{j,v}\mathrm{}`$) and differential rotation with constant specific angular momentum or constant linear velocity (for $`A_{j,v}0`$), respectively.

$$\mathrm{\Omega }=\frac{A_v+1}{A_v+\sqrt{1\mu ^2}}\mathrm{\Omega }__\mathrm{E}.$$
(9)
In both cases $`\mathrm{\Omega }__\mathrm{E}`$ is the angular velocity at the equator ($`\mu =0`$) and the parameter $`A_{j,v}>0`$ accounts for the degree of differential rotation so that uniform rotation (i.e. $`\mathrm{\Omega }/\mathrm{\Omega }__\mathrm{E}=\mathrm{const}.`$) is reached for $`A_{j,v}\mathrm{}`$. Note that the use of a law of differential rotation is physically plausible as long as such a law does not violate Rayleigh’s stability criterion for rotating inviscid fluids: $`d\left(\varpi ^2\mathrm{\Omega }\right)^2/d\varpi >0`$, where $`\varpi `$ is the cylindrical radial coordinate. It is straightforward to check that both expressions (8) and (9) satisfy Rayleigh’s criterion for all values of $`A_{j,v}>0`$.

The solution of equation (6) is considerably simpler than the solution of the corresponding set of partial differential equations describing normal modes of a differentially rotating axi-symmetric star and discussed by KYE. Nevertheless, most of the mathematical properties that are found in the solution of the system of partial differential equations can already be found in (6). A particularly important property of equation (6), which has been encountered also when dealing with $`r`$ modes of differentially and rapidly rotating stars (KYE), is that it may become a singular eigenvalue problem at the angular position $`\mu _s`$ for which

$$\sigma m\mathrm{\Omega }(\mu _s)=0=\omega _{\mathrm{ph}}\mathrm{\Omega }(\mu _s),$$
(10)
where $`\omega _{\mathrm{ph}}\sigma /m`$ is the phase velocity of the mode. The condition (10), which can be interpreted as the appearance of a “corotation point” (i.e. a point on the shell at which the perturbation pattern rotates at the same angular velocity as the background shell) has been interpreted by KYE as the cause impeding the numerical solution of the eigenvalue problem for a sufficiently large degree of differential rotation222KYE also noted that the difficulty in finding a convergent numerical solution was dependent on the rate of rotation of the background models, with rapidly rotating models providing solutions for comparably smaller values of $`A_{j,v}`$.. Interestingly, equation (6) offers close analogies also with the corresponding equation obtained in full General Relativity for a slowly and uniformly rotating relativistic star (Yoshida 2001, Yoshida and Futamase 2001; Ruoff and Kokkotas, 2001a, 2001b). The analogy is brought about by the fact that the corrections due to the relativistic dragging of inertial frames are mathematically similar to the corrections due to differential rotation and introduce a similar coefficient which could become singular for certain rates of rotation333Note that although equation (6) can in principle be singular, a global solution can still be found if series expansion techniques, such as the Frobenius method, are employed in the vicinity of the singular point $`\mu _s`$..

## 3 Strategy of the Numerical Solution

The numerical procedure adopted in the solution of equation (6) is straightforward and needs particular care only if a singular point should appear during the solution. In general, we can numerically integrate equation (6) as a two point boundary value problem with boundary conditions at both poles, $`\mu =\pm 1`$, being given by the requirement that the solution $`\chi `$ is regular there. The two solutions found are then matched at an arbitrary point $`\mu __\mathrm{M}`$ in the domain $`(1,1)`$ following the standard procedure in the solution of an eigenvalue problem. In particular, for a trial value of the eigenfrequencies we look for a zero of the Wronskian evaluated at $`\mu __\mathrm{M}`$, with different zeros representing sequences of different mode numbers. Once a zero is found, the eigenfrequency is used in equation (6) to calculate the corresponding eigenfunction. This procedure is repeated for different degrees of differential rotation so that a sequence is built describing the eigenfrequencies of the rotating shell from uniform rotation up to extreme differential rotation. The results of these calculations are discussed in detail in the following Section.

## 4 Numerical Results

This Section briefly presents the results of the numerical solution of equation (6) for the two laws of differential rotation discussed above.

### 4.1 $`j=\mathrm{const}.`$ differential rotation

We show in figure 1 the values of the phase velocity $`\omega _{\mathrm{ph}}`$ normalized to the value of the angular velocity at different latitudes on the shell. (This normalization is particularly convenient as it allows one to detect immediately whether corotation, corresponding to $`\omega _{\mathrm{ph}}/\mathrm{\Omega }(\mu )`$ 1, takes place or not.) As expected, in the limit of uniform rotation, the curves converge to $`2/3`$, the value given by the dispersion relation (7) for an $`m=l=2`$, $`r`$ mode. For progressively smaller values of $`A_j`$, on the other hand, the curves split in response to the different angular velocities at different angular positions on the shell. Surprisingly, corotation is reached only at the equator (i.e. $`\mu =0`$) and then only in the limit of $`A0`$, which is not physically interesting.

This can be clearly seen in figure 2, where we show the deviation away from unity of the ratio $`\omega _{\mathrm{ph}}/\mathrm{\Omega }(\mu =0)`$. For both curves, which correspond in the limit of uniform rotation to an $`m=l=2`$ and an $`m=2`$, $`l=3`$ mode respectively, the behaviour for small values of the parameter $`A`$ is well fitted by a power law $`\omega _{\mathrm{ph}}/\mathrm{\Omega }(\mu =0)KA^n`$, where $`n3.76`$ for $`m=2`$, ’$`l=2`$’, and $`n53.4`$ for $`m=2`$, ’$`l=3`$’ ($`K`$ is here a positive constant). As a result, the deviation away from unity tends to zero only in the limit of $`A0`$ and we therefore conclude that no corotation appears in these modes in the range of validity of a $`j=\mathrm{const}.`$ law of differential rotation.

The small inset in figure 1 shows a comparison between the ratio $`\omega _{\mathrm{ph}}/\mathrm{\Omega }(\mu =0)`$ obtained with the present shell approach and the corresponding quantity (indicated with filled squares) obtained by KYE with their numerical code for a slowly rotating stellar model with axis ratio 0.95. There is a rather good agreement, at least for the range of differential rotation rates in which the mode calculation was possible. Such an agreement brings confidence about the relevance of the results obtained with the shell model also for multidimensional star models.

When looking at the results of KYE, it becomes apparent how the behaviour of the eigenfrequencies might have suggested the appearance of a corotation point. As mentioned in the Introduction, when the degree of differential rotation increases past a certain threshold, very large radial and polar gradients appear in the equations for the eigenvalue problem. The solution of very stiff equations using finite difference techniques is a very difficult task and it is therefore not surprising that KYE were not able to find convergent solutions to the eigenvalue problem for $`A0.3`$. By making use only of an ordinary differential equation, the shell approach bypasses this difficulty and provides an accurate solution for any value of $`A`$.

### 4.2 $`v=\mathrm{const}.`$ differential rotation

In the case of the $`v=\mathrm{const}.`$ law, the absence of a corotation point is even more evident. In figure 3 the eigenfrequency is plotted as in figure 1. There, we clearly see that the ratio $`\omega _{\mathrm{ph}}/\mathrm{\Omega }(\mu )`$ is below unity everywhere on the shell. The comparison with the results by KYE, shown in the small inset of figure 3, is less good than the one seen for a $`j=\mathrm{const}.`$ differential rotation law but the overall behaviour is rather similar.

In addition to the differential rotation laws (8) and (9) we have also investigated a quadratic differential rotation law of the type $`\mathrm{\Omega }/\mathrm{\Omega }__\mathrm{P}=1+(\mu ^21)/B`$, where $`B>0`$ and uniform rotation is reached for $`B\mathrm{}`$. Also in this case, no evidence for corotation was found for the values of the parameter $`B`$ satisfying the Rayleigh stability criterion (i.e $`B>2`$).

## 5 Conclusion

We have used a thin-shell model to investigate the behaviour of Rossby-Haurwitz waves (the analogues of $`r`$ waves on thin shells) in the presence of differential rotation. Our simplified approach replaces the set of partial differential equations for a multidimensional stellar model with a single ordinary differential equation that can be solved to arbitrary accuracy. The numerical solution of the eigenvalue problem has shown that there is no evidence of corotation, even for asymptotically large values of differential rotation. As a result, the eigenvalue problem for Rossby-Haurwitz waves never becomes singular, as had instead been suggested by recent calculations of differentially rotating axi-symmetric Newtonian stars. We have found that, for this simplified model, $`r`$ modes can in principle coexist with arbitrarily large rates of differential rotation, the eigenfrequencies and eigenfunctions being suitably modified in response to the degree of differential rotation.

The relevance of the results found here for more realistic stellar models is not easy to assess, although the equations solved here show many of the mathematical features of the corresponding equations for Newtonian differentially rotating stars, or for slowly rotating relativistic stars.

It is pleasure to acknowledge fruitful discussions with M. Abramowicz, N. Andersson, S. Bonazzola, Y. Eriguchi, K. Kokkotas, J. Miller, J. Ruoff and N. Stergioulas. Financial support for this research has been provided by the MURST and by the EU Programme “Improving the Human Research Potential and the Socio-Economic Knowledge Base” (Research Training Network Contract HPRN-CT-2000-00137).

## References