text:
# Elasticity of Stiff Polymer Networks

## Abstract

We study the elasticity of a two-dimensional random network of rigid rods (“Mikado model”). The essential features incorporated into the model are the anisotropic elasticity of the rods and the random geometry of the network. We show that there are three distinct scaling regimes, characterized by two distinct length scales on the elastic backbone. In addition to a critical rigidiy percolation region and a homogeneously elastic regime we find a novel intermediate scaling regime, where the elasticity is dominated by bending deformations.

The elasticity of cells is governed by the cytoskeleton, a partially crosslinked network of relatively stiff filaments forming a several $`100`$ nm thick shell called the actin cortex alberts . While the statistical properties of single cytoskeletal filaments are by now relatively well understood wilhelm-frey:96 ; legoff-etal:02 , theoretical concepts for the elasticity of stiff polymer networks are still evolving. One major open question is to understand how stresses and strains are transmitted in such networks. In synthetic gels that are formed by rather flexible chain molecules the response to macroscopic external forces is – on the level of single filaments – isotropic and entropic in origin. It is generally believed that macroscopic stresses are transmitted in such a way that the local deformations within the network stay affine, i.e. that the end-to-end distance of individual filaments follows the macroscopic shear deformation doi-edwards:86 . In contrast, the building blocks of the actin cortex are semiflexible polymers, whose hallmark is an extremly long persistence length $`\mathrm{}_p`$, which is comparable to the total contour length $`\mathrm{}`$. As a consequence, the response of such stiff polymers to external forces shows a pronounced anisotropy kroy-frey:96 . Consider a semiflexible polymer with one end clamped at a fixed orientation. When forces are applied at the other end transverse to the tangent vector at the clamped end, the response may be characterized by a transverse spring coefficient $`k_{}(\mathrm{})=3\kappa /\mathrm{}^3`$ proportional to the bending modulus $`\kappa `$. Whereas this response is of purely mechanical origin, the linear response for longitudinal forces is due to the presence of thermal undulations, which tilt parts of the polymer contour with respect to the force direction. The corresponding effective spring coefficient $`k_{}(\mathrm{})=6\kappa ^2/(k_BT\mathrm{}^4)`$ is proportional to $`\kappa ^2/T`$ indicating the breakdown of linear response for very stiff filaments. In a typical network one expects the distance between crosslinks $`\mathrm{}_c`$ to be much smaller than the persistence length and filament length. Hence we have $`k_{}(\mathrm{}_c)/k_{}(\mathrm{}_c)=2\mathrm{}_p/\mathrm{}_c1`$, i.e. the elastic response of the filaments is indeed highly anisotropic.

These anisotropic elastic properties of individual filaments suggests that the macroscopic elasticty of networks will not only depend on the number of crosslinks and the density of filaments, but also on the geometry and architecture of the network. For very regular networks such as a triangular lattice the longitudinal spring coefficient $`k_{}`$ will dominate the macroscopic moduli mackintosh-kaes-janmey:95  since the network can not be strained without a change of the end-to-end distance of individual polymers. In other regular network architectures, the softer bending modes would be dominant satcher-dewey:96 . Naturally, this will lead to a very different prediction for the elastic modulus of the network. It is not at all obvious what type of network geometry (elongation dominated versus bending dominated) is relevant in less ideal structures with a significant amount of disorder such as in typical cytoskeletal networks.

As a first step towards understanding the elasticity of stiff polymer networks we consider a two-dimensional model defined as follows (see Fig. 1).

We generate the random network by placing $`N`$ line-like objects of equal length $`\mathrm{}`$ on a plane with area $`A=L^2`$ such that both position and orientation of the filaments are uniformly randomly distributed. Periodic boundary conditions in both directions are used. Upon increasing the line density $`\rho =N\mathrm{}/A`$ there is a critical threshold $`\rho _c`$ for geometric percolation pike-seager:74a . Numerical simulations leroyer-pommiers:94  show that the correlation length $`\xi (\rho \rho _c)^\nu `$ of the incipient infinite percolation cluster scales with a critical exponent $`\nu =4/3`$, identical to the value obtained for random site percolation on a lattice stauffer-aharony:94 . Transport of scalar quantities like the conductivity is also in the same universality class as lattice models balberg-binenbaum-wagner:84 . In order to study the transport of non-scalar quantities like shear stress we need to specify how forces are transmitted between the building blocks of the network. In our “Mikado model” the building blocks are homogeneous elastic rods characterized by a Young modulus $`E`$ and a circular cross-section of radius $`r`$. Wherever two sticks intersect they are connected by a crosslink with zero extensibility. In the cytoskeleton one finds a variety of linker proteins with a range of mechanical properties limozin-sackmann:02 . Here we restrict ourselves to crosslinks that either fix the relative orientation of the rods (“stiff crosslinks”) or allow free rotation (“free hinges”). Similar to thermally fluctuating semiflexible polymers, the elastic response of a stick segment between two neighboring crosslinks is characterized by length dependent force constants for compression or elongation, $`k_{\text{comp}}(\mathrm{}_c)=\pi r^2E/\mathrm{}_c`$, and bending, $`k_{\text{bend}}(\mathrm{}_c)=k_{}(\mathrm{}_c)=(3/4)\pi r^4E/\mathrm{}_c^3`$. The distance between two crosslinks $`\mathrm{}_c`$ shows a Poissonian distribution, where the average distance of crosslinks along a filament scales as the inverse of the line density, $`\overline{\mathrm{}_c}=\pi /\rho `$ pike-seager:74a . While this is a purely mechanical model, it still captures the essential feature that for typical densities of the network the compressional stiffness is much larger than the bending stiffness, $`k_{\text{comp}}(\mathrm{}_c)/k_{\text{bend}}(\mathrm{}_c)=(4/3)\mathrm{}_c^2/r^21`$. It leaves out steric effects due to thermal fluctuations of the filaments, which give rise to the plateau modulus in solutions hinner-etal:98 .

Consider the energy of the network as a function of the deviations of the positions of all intersection points and the rod orientations at the intersection points from their initial values. For small deformations of the network, this function can be approximated by a quadratic form that vanishes for vanishing deviations, as—by construction—the undeformed network is not prestressed. To analyse the elastic properties of the model network, a shear deformation respecting the periodic boundary conditions is enforced by demanding that corresponding points on the left and right boundary of the simulation cell undergo equal displacements while the displacements of corresponding points on the upper and lower boundary of the cell must agree vertically but differ horizontally by a distance $`\mathrm{\Delta }=\gamma L`$, where $`\gamma `$ is the shear strain. The orientation of the rods at corresponding points on the boundary are required to be equal. The remaining degrees of freedom are then allowed to relax, i.e. the harmonic approximation to the energy of the network is minimized in the presence of the constraints. The derivative of the resulting energy of the deformed state with respect to the strain $`\gamma `$ is proportional to the shear modulus. In pinciple, this reduces the determination of the shear modulus of a given network to the solution of a linear equation. However, for interesting parameters (thin rods), the problem is numerically highly unstable as we are searching for the lowest point of a complicated high-dimensional valley with extremely steep slopes but hardly varying base altitude. This problem is best left to one of the commercially available finite element solvers which have seen many years of careful optimization and testing. The results presented below were obtained using the program Nastran by MSC Software.

In the following discussion we take the rod length $`\mathrm{}`$ as unit length and $`\kappa /\mathrm{}^3`$ as unit for the elastic modulus. Then the independent parameters are the densitiy $`\rho `$, the system size $`L`$ and the aspect ratio $`\alpha =r/\mathrm{}`$ of the rods. Note that the latter is a measure of the relative magnitude of compressional to bending stiffness.

We start with an analysis of the elasticity close to the percolation threshold. For stiff crosslinks we find that the percolation threshold is the same for rigidity as for connectivity percolation, $`\rho _c=5.71`$. For free hinges a higher line density of filaments is needed, $`\rho _c=6.7`$, for the network to become rigid. This agrees well with recent results, $`\rho _c=6.68`$, for stiff fiber networks latva-kokko-timonen:01 , where the crosslinks are fixed in space but the angles between the fibers can vary. In both cases, we find that the shear modulus $`G`$ vanishes as the line density of sticks approaches the critical value $`\rho _c`$, according to a power law $`G(\rho \rho _c)^\mu `$. For our numerical analysis on a finite lattice we expect the shear modulus to obey the following finite size scaling law

$`G=L^{\mu /\nu }h(L/\xi ),`$ (1)
where the scaling function behaves as $`h(x)x^{\mu /\nu }`$ and $`h(x)1`$ for large and small values of the scaling variable $`x=L/\xi `$, respectively. Fig. 2 shows that the data collapse works very well for densities ranging from values close to the percolation threshold $`\rho _c`$ up to $`\rho 20`$. For the data shown, $`L`$ ranges from $`2`$ to $`30`$.

For larger densities, systematic deviations are clearly visible. This will turn out to be a very interesting observation, as we will discuss in detail below. We get the best data collapse in the critical region if we choose the values $`2.4\pm 0.2`$ and $`2.3\pm 0.2`$ for the critical exponent $`\mu /\nu `$ in the case of stiff crosslinks and free hinges, respectively. Since the difference between the exponents is within the statistical error, we can make no definite conclusion whether networks with free hinges and stiff crosslinks belong to different universality classes for elasticity percolation.

The rigidity exponent $`\mu 3.15\pm 0.2`$ is significantly lower than in other classes of continuum percolation models, like the “Swiss-cheese model”, where $`\mu 5`$ feng-halperin-sen:87 ; benguigui:86 . It is also lower than the value $`\mu 4`$ for lattice models with bond-bending forces arbabi-sahimi:93b ; stauffer-aharony:94 . Hence it seems likely that the “Mikado model” constitutes a new universality class for rigidity percolation. Similar results have independently been found in Ref. head-levine-mackintosh:03 .

Now we come back to the above mentioned systematic deviations from the scaling law, Eq. 1, at densities above $`\rho 20`$. To understand these better, let us have a closer look at the shear modulus as a function of the elastic moduli of individual filaments for densities not too close to the percolation threshold. In this regime the shear modulus becomes independent of system size for moderately large systems; for the following results we have chosen systems satisfying $`L/\xi 200`$. Fig. 3 shows the shear modulus as a function of $`\alpha `$ for a series of densities; we have communicated a preliminary version of these data in Ref. frey-kroy-wilhelm-sackmann:97 . Note that $`k_{\text{bend}}(\mathrm{})`$ is effectively kept constant since we are measuring all elastic constants in units of $`\kappa /\mathrm{}^3`$.

There are two strikingly different regimes. For high densities and/or thick rods ($`\alpha 0.1`$), where compressional stiffness is lower or comparable to the bending stiffness (lower right part of Fig. 3), the shear modulus scales linearly with the filament compressional modulus and the number of filaments per unit area, $`G(\rho \rho _c)\alpha ^2`$. Such a linear regime has also been found in a series of studies on random fiber networks raisanen-etal:97 . It is by now well established that the elastic modulus can be described quantitatively in terms of effective medium models astrom-etal:00 . Hence, in the high line density regime the network behaves as a homogeneously elastic medium, dominated by the compressional modulus of the individual filaments. As a consequence, local deformations follow a macroscopic shear in an affine way. This has to be contrasted with the elastic behavior for slender rods with low aspect ratios ($`\alpha 10^5`$ for the higher densities), where bending becomes the softer mode. Then, one finds an extended plateau region, which broadens significantly with lowering the line density, where the shear modulus becomes completely independent of $`k_{\text{comp}}(\mathrm{})\alpha ^2k_{\text{bend}}(\mathrm{})`$ frey-kroy-wilhelm-sackmann:97 . This strongly suggests that in this regime the macroscopic elasticity of the network is dominated by bending stiffness of individual filament. This conclusion is corroborated by the observation that almost all of the energy stored in the deformed network is accounted for by transverse deformation of the rods (compare Fig. 1). Another remarkable feature of this plateau regime is the strong dependence of the shear modulus on line density. We find $`G(\rho \rho _c)^\mu ^{}`$ with a rather large exponent $`\mu ^{}6.7`$. From the above analysis it may seem as if the anomalous elasticity in the plateau regime and the homogeneous elasticity in the affine regime are two separate phenomena, and one might wonder how one emerges from the other. To analyze this relation we try a crossover scaling ansatz,

$$G=(\rho \rho _c)^\mu ^{}g[\alpha (\rho \rho _c)^\nu ^{}]=\xi ^{\mu ^{}/\nu ^{}}\stackrel{~}{g}(\alpha /\xi ^{}),$$
(2)
where we have defined a new length scale $`\xi ^{}(\rho \rho _c)^\nu ^{}`$. For this ansatz to reduce to the modulus expected in the affine region, the scaling function $`g(x)`$ needs to scale as $`g(x)x^2`$ for $`x1`$ and the exponents need to obey the scaling relation $`\mu ^{}=2\nu ^{}+1`$. In the plateau regime, $`g(x)`$ is expected to be constant. As shown in Fig. 4, we obtain an excellent scaling collapse for over almost eight orders of magnitude in the scaling variable $`x=\alpha /\xi ^{}`$ using $`\nu ^{}=2.83`$ or equivalently $`\mu ^{}=6.67`$ and the critical line density $`\rho _c5.71`$, associated with connectivity percolation. Additionally, the scaling function $`g(x)`$ displays the expected behavior. Meeting both of these requirements is highly nontrivial, and gives strong evidence for the anomalous scaling law in Eq. 2.

The existence of such a broad scaling regime far from the percolation threshold is a surprising and intriguing feature of stiff polymer networks. At the moment we are lacking a complete understanding of its physical origin. In particular, the geometrical significance of the new length scale $`\xi ^{}`$ remains unclear. One may speculate that the anomalous scaling behavior is a subtle consequence of the interplay between quenched randomness of the network structure and long-range correlation effects induced by the stiffness of the filaments. An immediate consequence of the scaling form, Eq. 2, is the existence of a crossover line density $`\rho _{\text{cross}}`$ scaling as $`\mathrm{}\rho _{\text{cross}}\alpha ^{1/\nu ^{}}`$, where we have re-introduced units of length $`\mathrm{}`$. This implies that increasing filament length at constant line density drives the system towards the affine regime, in accord with Ref. head-levine-mackintosh:03 .

While these results for an idealized two-dimensional model are certainly not straightforwardly applicable to three-dimensional cytoskeletal networks, one may still try to get an idea for the scales involved. We expect that network densities can be compared roughly by using the average distance $`\mathrm{}_c`$ between intersections as a measure: A cytoskeletal network might have $`\mathrm{}_c0.1\mu \text{m}`$ with typical filament lengths of $`2\mu \text{m}`$. These values correspond to a two-dimensional line density of $`\rho 20`$ and an aspect ratio of $`\alpha 0.002`$, which would place a typical actin network in the bending dominated intermediate regime.

Understanding the full complexity of cytoskeletal networks certainly merits further theoretical and experimental work. Building on the knowledge gained from our idealized model, future investigations may among many other questions want to address three-dimensional systems, polydispersity, thermal fluctuations or even the kinetics of the crosslinking molecules.

We acknowledge M. Alava and K. Kroy for useful discussions and comments, and P. Benetatos for a critical reading of the manuscript.


format: text

text:
# Preferential attachment in the protein network evolution

## Abstract

The Saccharomyces cerevisiae protein-protein interaction map, as well as many natural and man-made networks, shares the scale-free topology. The preferential attachment model was suggested as a generic network evolution model that yields this universal topology. However, it is not clear that the model assumptions hold for the protein interaction network. Using a cross genome comparison we show that (a) the older a protein, the better connected it is, and (b) The number of interactions a protein gains during its evolution is proportional to its connectivity. Therefore, preferential attachment governs the protein network evolution. Evolutionary mechanisms leading to such preference and some implications are discussed.

The analysis of networks has attracted great interest in recent years. Many man-made networks, including the World Wide Web1 , scientific2  and movie actor3  collaborations, and linguistic4  networks, have been shown to be scale free, with different nodes having widely different connectivities5 ; 6 ; 7 . Networks of biological origin, such as metabolic interaction8  and protein-protein interaction networks9 , also share this property. The emergence of the scale-free topology in such diverse examples calls for a universal explanation, based on generic principles, applicable to all the different networks studied. This was achieved by the growing network model, suggested by Barabási and Albert10 , which assumes the continuous creation of new nodes and their preferential attachment to previously well-connected nodes. An exact solution for the dynamics of the model demonstrates the emergence of the scale-free topology from these generic assumptions, given an asymptotically linear attachment kernel11 ; 12 . The model assumptions seem self-evident for social networks. A direct test for some of these networks have validated the preferential attachment principle, and shown an approximate linear kernel13 ; 14 . However, it is less clear how this model can be justified for natural networks, such as the biological networks. While the dynamic growth of the network can be understood on an evolutionary time scale10 , the preferential attachment assumption is far from obvious, as the interactions are not formed based on a conscious choice.

In this work, we focus on the Saccharomyces cerevisiae (bakers’ yeast) protein-protein interaction network, which is often used as a model for a biological interaction network. A cross-genome comparison is employed to obtain a classification of the yeast proteins into different age groups. We observe a correlation between a protein’s age and its network connectivity, in accordance with the growing network picture. Furthermore, this classification enables us to directly observe the preferential attachment phenomenon. Signs of this phenomenon have been previously observed through analysis of divergent pairs of duplicated genes wagner . We thus conclude that the Barabási-Albert model is indeed relevant for describing the evolution of the yeast protein-protein interaction map. We further discuss implications of this phenomenon to the governing rules of protein evolution.

We start by classifying the whole database of 6294 bakers’ yeast proteins15  into four age groups. For this purpose, we pick three other model organisms for which a fully sequenced genome and a comprehensive list of proteins are available, and are of varying evolutionary distance from the baker’s yeast. The evolutionary distance between two organisms can be extracted from the phylogenetic tree (the ”tree of life”) describing the evolutionary branching process16  (see figure 1): Escherichia coli17  belongs to the Bacteria branch (estimated time of diversion 4 Giga-years ago, Gya), Arabidopsis thaliana18  belongs the Plants branch (estimated diversion 1.6Gya), while Schizosaccharomyces pombe19  (fission yeast) and the bakers’ yeast belong to different sub-phyla on the Fungi branch (estimated diversion 1.1Gya). A cross-genome comparison between these organisms is employed in order to estimate the age of each bakers’ yeast protein. We assume that a protein created at a certain time in a certain ancestor organism will have descendants in all organisms that diverged from this ancestor. For example, proteins that are older than the first (Bacteria) diversion should have descendants in all four organisms, while those created after the fission-yeast diversion are expected to have descendants in the bakers’ yeast alone. While the descendant proteins continue to evolve and diverge, they still show higher sequence similarity than a random pair of proteins.

For each of the bakers’ yeast proteins, we search for similar proteins in the other three organisms (see details below), and use the results to classify it into one of four age groups. Proteins with no fission-yeast similarities are expected to be relatively new (group 1, 872 proteins); those with similarities only in fission-yeast are expected to have an ancestor prior to the diversion and are therefore older (group 2, 665 proteins); those with fission-yeast and Arabidopsis similarities are even older (group 3, 2079 proteins); and those with analogues in all three organisms form the oldest group of proteins (group 4, 2678 proteins), with ancestors that predate the first diversion. Only a small fraction (less than 10%) of the similarities were not consistent with the evolutionary timeline. Note that our age-group classification is not sensitive to duplication events duplications , and thus new proteins generated by duplication are here classified as old.

Here are some brief technical details on the similarity search done. We use the standard definitions for the similarity distance between sequences, and employ the standard Protein-BLAST program20 . The program is given a query sequence (in our case: the yeast protein) and a reference database (the set of all proteins of the other organism), and compares the query sequence to each of the database sequences, in search for shared patterns. Each found match gets a score (termed “E-score”), which is the expected number of same or higher quality matches given a randomized database. The probability to get a match of same or higher quality for a random sequence is

$$P(E|\mathrm{random}\mathrm{pair})=1\mathrm{exp}(E),$$

where $`E`$ is the E-score. The lower this probability, the higher the confidence that the sequences similarity (or the match) is indeed due to a common ancestor for both sequences. We considered two proteins to be similar if the E-score of their match was lower than the cutoff value $`E_c=0.7`$, corresponding to $`P(E_c|\mathrm{random}\mathrm{pair})0.5`$.

In the following, we use the obtained age-group classification of the yeast proteins to analyze the structure of the protein-protein interaction network. We use a published database of yeast protein-protein interactions21 , and first look at the average connectivity. Figure 2 shows a clear dependency of the connectivity on the protein age, with older proteins having significantly more interactions. While group 1 proteins (newest) have only 0.5 links per protein, group 4 proteins (oldest) have 6.2 links per protein. This supports the picture of the growing network model, where the older a node the higher its probability to gather interactions with other late-coming proteins.

A direct test of the second assumption of the growing network model, namely, the preferential attachment principle, requires detailed information on the network development, which is beyond our reach. However, the above classification provides us with snapshots of the growing network at three points in its evolution, enabling an insight into the evolution of protein interactions. We study the sub-network defined by group 4 proteins and the links connecting them, recording the connectivity of each old protein on this sub-network. This sub-network was used as a model for the interaction map at an early stage of the evolution process (the time of divergence of the Bactria branch). The number of links of each old protein to the newer proteins (groups 1,2,3) is the number of links acquired since that time. We then looked at the number of new links a node gathered as a function of its connectivity in the old network. A similar analysis is done for the sub-networks defined by groups 3 and 4 combined (proteins with an Arabidopsis analogue), and for groups 2,3 and 4 combined (proteins with fission-yeast analogue). As Figure 3 shows, the number of new links tends to increase with the number of links in the old network, which is a signature of preferential attachment. The number of new links appears to be approximately linear in the connectivity, suggesting a linear preferential attachment kernel, and consistent with the scale free topology11 .

The growing network paradigm suggests a dynamic model for preferential attachment: that is, all nodes are created equal and the attachment probability is related to the actual current connectivity (”rich get richer” model) as defined by the network dynamics. An alternative model22  suggests a static explanation in which each node has a different intrinsic fitness that determines its ability to interact and doesn’t change as the network grows. In this model both the actual connectivity and the attachment probability of a protein depend on its intrinsic fitness. Given an appropriate distribution of the fitness parameter, this model can explain the results of figure 3 (”good get richer” model), but it is not consistent with the age-dependence shown in figure 2. While the growing network model predicts that older nodes will be better connected, connectivity in the static model is related solely to the node fitness, and age and connectivity shouldn’t be correlated. Thus, our results (figure 2) support the first option as a model for the protein interaction evolution. Gene duplication was also suggested as an explanation for the scale-free topology of the protein interaction network23 ; 24 ; 25 ; 26 . However, since duplication events are not detected by our age-group classification, our results show that the proteins network structure cannot be attributed solely to evolution by duplication.

The question of the evolutionary mechanism leading to the dynamic preferential attachment remains: how does becoming better connected make a protein more attractive for future interactions, and why is the preference linear in the number of connections? We suggest two possible mechanisms that partially answer these questions:

(i) The more connections a node acquires, the stronger is the selective pressure to make it more connectable. On the molecular level this can be understood as a tendency to increase the number of protein attachment domainsdomains  (such as the WWWW  or proline-richproline  domains), or to improve the existing domains such that they bind to more target proteins. In this mechanism the preferential attachment is related to the physico-chemical properties of the highly-connected protein. In order to test this possibility, one can look at the distribution of domains and other reoccuring patterns in the set of highly-connected proteins, and check whether connectability can be traced to sequence motifs. However, the lack of well-studied interaction network for other organisms and the partial understanding of attachmnet properties of protein domains limits our ability to perform such study.

(ii) Many protein interactions are actually physical interactions that change or regulate the functionality of the interacting parties, such as phosphorylation and complex formation. The number of potential distinct operation modes of a protein increases exponentially with the number of its regulating proteins, and similarly the number of potential variants of a given complex increases exponentially with the number of its building-block proteins. Therefore, the more connected a protein, the stronger the selection towards creating a protein to interact with it. Here, the phenomena relates to the biological functionality of the protein. This mechanism can be validated by the following experiment: current technology enables us to dig out proteins that form a complex together with a given target proteintap . One can look at the different complexes generated under varying conditions and study the different combinations obtained, that is, how many distinct complexes were formed using the target protein. Then, it is possible to study how many new structures have been made available by each complex member. We predict that the contribution of each new member will be multiplicative, i.e., the number of new structures will be, on average, proportional to the total number of structures.

The preferential attachment phenomenon demonstrates an important principle in the process of evolution. It dynamically leads to the formation of big protein complexes and pathways, which introduce high complexity regulation and functionality. New systems are not generated as self-interacting modules of new proteins; rather, new proteins tend to connect to the old well-connected hubs of the network and modify existing functional units. Indeed, 267 of the 872 group 1 proteins (31%, versus 12% of group 4) have no interactions documented in the database, indicating a very low number of actual interactions. Thus, we get information on protein’s centrality based on its sequence alone. This information is helpful in analyzing the protein interaction network given the partial information available.

It was shown that the higher the connectivity of a node, the higher its probability to be essential, i.e. to have a lethal knockout phenotype9 . As mentioned above, highly connected nodes tend to be older. We find that essential proteinss also tend to be older: only 8% of the newest proteins are essential, in contrast to 20% of the oldest proteins ($`\chi ^2`$-test p-value $`310^{20}`$).

In conclusion, we show that the protein networks evolve by creating new, unconnected links, which attach to the existing network according to the linear preferential attachment principle. This explains the scale-free topology shared by the network, and has implications for understanding the evolutionary mechanisms. The correlation of the protein’s age to its centrality opens new possibilities for deriving information on the interaction network topology based on sequence data.

###### Acknowledgements.

We thank A.-L Barabási, Y. Kliger, A. Lipshtat and S. Redner for valuable discussions and comments.


text:
# SL(2,Z) tensionless string backgrounds in IIB string theory

## 1 Introduction

In string theory the analogue of a massless particle is a tensionless string . Just like particles often can be treated as approximately massless in a high energy description, tensionless strings are expected to be of interest for the high energy behavior of string theory . Various action formulations with or without supersymmetry have been constructed , and their quantum properties have been studied . One interesting result is that if the (spacetime) conformal symmetry of the classical formulation is to survive quantization, the spectrum is topological . (This brings to mind the comments by Witten on the relation of the results in  to topological theories ).

It should also be noted that in the AdS/CFT context, an interesting relation between the spectrum of tensionless strings and higher spin theories has been advocated .111cf.  and references therein for an overview on highes spin theories. More abstractly, tensionless strings have recently been seen to arise in the framework of intersecting M5-branes  and earlier in the context of D-branes , although for these tensionless strings an action is lacking.

An important question is the role of these strings in the solution space of supergravity. In this brief note we show that they arise as sources to exact solutions therein. It proceeds as follows: First, we describe Schwarz’ classical solution to $`D=10`$ IIB supergravity theory  in terms of a set of $`SL(2,)`$ background metrics, which is the starting point for our calculation. Then, in section 3, we perform a Lorentz transformation orthogonally to the string. We derive the limit as the boost achieves the speed of light. It turns out that the solution is Minkowski space up to a plane-fronted shock wave at the position of the (boosted) string. We show that this gives a solution to the string equations to all orders and that the limit is compatible with Einstein’s equations. We derive the energy momentum tensor and show that is has the expected behavior. We end with a few comments and open questions for future study.

## 2 The original solution

We consider $`D=10`$ IIB supergravity theory with the bosonic action

$`S_{10}^{\mathrm{IIB}}={\displaystyle \frac{1}{2\kappa ^2}}{\displaystyle \mathrm{d}^{10}x\sqrt{g}\left[R+\frac{1}{4}\text{tr}(^1)\frac{1}{12}^T\right]},`$ (1)
where $`F_5`$ is put to zero.222The charges of the $`F_5`$ field strength are carried by 3-branes, which is not of interest here. $``$ is a $`SL(2,)`$ matrix of the scalar fields and $``$ is a vector of the 3-forms:

$`=\mathrm{e}^\varphi \left(\begin{array}{cc}\left|\lambda \right|^2& \chi \\ \chi & 1\end{array}\right),=\mathrm{d}=(H^{(1)},H^{(2)}),`$ (4)
with $`\lambda =\chi +\mathrm{ie}^\varphi `$. Here, $`\varphi `$ is the dilaton (NS-NS) field and $`\chi `$ is the R-R scalar field. The action (1) is invariant under the following action of $`SL(2,)`$:

$`\mathrm{\Lambda }\mathrm{\Lambda }^T,(\mathrm{\Lambda }^T)^1.`$ (5)
Schwarz generalized the results of Dabholkar et al.  and found a solution given by the $`SL(2,)`$ set of background metrics :

$`\mathrm{d}s^2=A_q^{3/4}\left(\mathrm{d}t^{\mathrm{\hspace{0.17em}2}}+(\mathrm{d}x^1)^2\right)+A_q^{1/4}\mathrm{d}x\mathrm{d}x.`$ (6)
where $`\mathrm{d}x=(\mathrm{d}x^2,\mathrm{},\mathrm{d}x^9)`$,

$`A_q=1+{\displaystyle \frac{\mathrm{\Delta }_q^{1/2}Q}{3r^6}},\mathrm{\Delta }_q^{1/2}=𝐪^T^1𝐪,𝐪=(q_1,q_2).`$ (7)
Herein, $`r^2=xx`$ and $`Q`$ is the fundamental $`_{\mu \nu }`$ electric charge. The restriction to $`SL(2,)`$ follows from the Dirac quantization condition and that $`q_1`$ and $`q_2`$ are relative prime. The corresponding fields read

$`_{01}=^1𝐪\mathrm{\Delta }_q^{1/2}A_q^1,`$ (8)
$`\lambda ={\displaystyle \frac{q_1\chi _0q_2|\lambda _0|^2+\mathrm{i}q_1\mathrm{e}^{\varphi _0}A_q^{1/2}}{q_1q_2\chi _0+\mathrm{i}q_2\mathrm{e}^{\varphi _0}A_q^{1/2}}},`$ (9)
where we used the fact that the only non-zero components of $``$ can be taken to be $`_{01}=_{10}`$ . The metric has a singularity at $`r=0`$, which is interpreted as an infinitely long source string with action

$`S={\displaystyle \frac{T_q}{2}}{\displaystyle \mathrm{d}^2\xi \left[\mathrm{\Delta }_q^{1/2}^aX^\mu _aX^\nu G_{\mu \nu }+ϵ^{ab}_aX^\mu _bX^\nu _{\mu \nu }^T𝐪\right]}.`$ (10)
Here, $`T_q=\mathrm{\Delta }_q^{1/2}Q`$ is the string tension. The background fields in (10) are actually string condensates. They arise as string loop effects . The above string solutions were interpreted by Witten as bound states of F(undamental) and D-strings .

## 3 Performing the boost

We study a boost of the metric (6) in a direction orthogonal to the string. Without loss of generality, we perform a Lorentz transformation in the $`z=x^9`$ direction. The form of the transformation is the usual

$`t^{}=\gamma (t+vz),z^{}=\gamma (x+vt),\gamma =\left(1v^2\right)^{1/2}.`$ (11)
In the transformed coordinates, the metric reads

$$\begin{array}{c}\mathrm{d}s^{}=\frac{A_q^{3/2}\mathrm{\Delta }_q^{1/2}Q\gamma ^2}{3r^6}\left(\mathrm{d}t^{}v\mathrm{d}z^{}\right)^2+A_q^{3/4}(\mathrm{d}x^1)^2\hfill \\ \hfill +A_q^{1/4}\left((\mathrm{d}t^{})^2+\mathrm{d}\stackrel{~}{x}\mathrm{d}\stackrel{~}{x}+(\mathrm{d}z^{})^2\right),\end{array}$$
(12)
where we introduced $`\stackrel{~}{x}=(x^2,\mathrm{}x^8)`$, $`r^6=\left(\gamma ^2(z^{}vt^{})^2+\stackrel{~}{x}\stackrel{~}{x}\right)^3`$ and

$`A_q^{}=1+{\displaystyle \frac{\mathrm{\Delta }_q^{1/2}Q}{3r^6}}.`$ (13)
If we boost (12) to the speed of light, i. e. let $`v1`$, we immediately find that the energy diverges. To avoid this problem, the energy has to be rescaled as $`ϵ=\gamma ^1ϵ_0`$ or, equivalently, $`Q=\gamma ^1Q_0`$, following the lines of Aichelburg and Sexl . We notice that $`\lambda `$, $``$ and hence, $`\mathrm{\Delta }_q`$ tend to their vacuum expectation values $`\lambda _0`$, $`_0`$ and $`\mathrm{\Delta }_{q,0}`$, respectively. This implies, in particular, that the dilaton field $`\varphi `$ and the R-R scalar field $`\chi `$ tend to their (constant) vacuum expectation values. The string tension becomes $`T_q=\mathrm{\Delta }_q^{1/2}Q\gamma ^1\mathrm{\Delta }_{q,0}^{1/2}Q_0^{}`$ and vanishes as $`v1`$. Thus, the result may be interpreted as a background metric with a tensionless string as its source. Finally, we find $`A_q^{}1`$ in the boost limit. Therefore, the only term we need to consider more carefully is the first term in (12), i. e. the limit of

$`\gamma r^6={\displaystyle \frac{\gamma ^5}{\left(\left(z^{}vt^{}\right)^2+\gamma ^2\rho ^2\right)^3}},`$ (14)
where $`\rho ^2=\stackrel{~}{x}\stackrel{~}{x}`$. Integrating both sides, one finds that the right hand side of (14) tends to

$`\underset{v1}{lim}\gamma r^6={\displaystyle \frac{3\pi }{8}}\left(\stackrel{~}{x}\stackrel{~}{x}\right)^{5/2}\delta \left(z^{}t^{}\right),`$ (15)
Thus, in the boost limit, the metric becomes

$`\mathrm{d}s^{}={\displaystyle \frac{\pi \mathrm{\Delta }_{q,0}^{1/2}Q_0^{}}{8\rho ^5}}\delta \left(z^{}t^{}\right)\left(\mathrm{d}t^{}\mathrm{d}z^{}\right)^2+(\mathrm{d}x^1)^2(\mathrm{d}t^{})^2+(\mathrm{d}z^{})^2+\mathrm{d}\stackrel{~}{x}\mathrm{d}\stackrel{~}{x},`$ (16)
or, equivalently, in coordinates $`u=z^{}t^{},v=z^{}+t^{}`$,

$`(\mathrm{d}s^{})^2=\mathrm{d}u\mathrm{d}v+(\mathrm{d}x^1)^2+\mathrm{d}\rho ^2+\rho ^2\mathrm{d}\mathrm{\Omega }+{\displaystyle \frac{\pi \mathrm{\Delta }_{q,0}^{1/2}Q_0^{}}{8\rho ^5}}\delta \left(u\right)\mathrm{d}u^2.`$ (17)
The solution is still an $`SL(2,)`$ set of background metrics. Moreover, it has the typical pp-wave metric structure

$`\mathrm{d}s^2=\mathrm{d}u\mathrm{d}v+K(x^2,\mathrm{},x^8,u)\mathrm{d}u^2+{\displaystyle \underset{i=1}{\overset{8}{}}}(\mathrm{d}x^i)^2.`$ (18)
For boosts in other directions we find an expected factor of $`\mathrm{sin}^2\varphi `$ in front of the coefficient function $`K`$, where $`\varphi `$ is the angle between the boost direction and the direction of the string, i. e. $`x^1`$. Especially, the only direction significantly different from the one described above ($`\varphi =\pi /2`$) is the direction parallel to the string ($`\varphi =0`$) where we recover flat Minkowski space and in which case the $`SL(2,)`$ degeneracy vanishes.

### 3.1 The B field

According to (8), we only consider the transformation of $`_{01}`$. The Lorentz transformation (11) generates four non-zero components in the transformed $`^{}`$ field, all of which diverge to infinity in the limit $`v1`$. This apparent problem is resolved by using the gauge freedom to perform the transformation

$`_{01}\stackrel{~}{}_{01}=_{01}^1𝐪\mathrm{\Delta }_q^{1/2}=^1𝐪\mathrm{\Delta }_q^{1/2}\left(A_q^11\right)\stackrel{v1}{}0.`$ (19)
The gauge fixed $``$ field vanishes in the limit of the boost, and hence does $`=\mathrm{d}`$.

### 3.2 Energy momentum tensor

The boosted metric $`G_{\mu \nu }`$ given by (17) is ”mostly flat” and contains a plane fronted shock wave at $`u=0`$, i. e. along the position of the string, encoded in the coefficient function $`K`$. This implies that $`G_{\mu \nu }`$ solves Einstein’s equations in empty space except for $`\rho =0`$ which contains a coordinate singularity. Moreover, since $`1/\rho ^5`$ is the Green’s function for the Laplacian in seven dimensions, i. e.

$`\mathrm{\Delta }{\displaystyle \frac{1}{\rho ^5}}={\displaystyle \frac{16}{3}}\pi ^3\delta (\rho ),`$ (20)
we obtain

$`R_{uu}={\displaystyle \frac{1}{2}}\mathrm{\Delta }K={\displaystyle \frac{1}{3}}\pi ^4\mathrm{\Delta }_{q,0}^{1/2}Q_0^{}\delta (\rho )\delta (u).`$ (21)
All other components of the Ricci tensor vanish. Thus, since the Ricci scalar vanishes, the energy momentum tensor of the system is read off from Einstein’s equations:

$`T_{uu}={\displaystyle \frac{1}{24}}\pi ^3\mathrm{\Delta }_{q,0}^{1/2}Q_0^{}\delta (\rho )\delta (u).`$ (22)
This solution has the expected form for an energy momentum tensor of an infinitely long (classic) string at rest. However, the result has to be compared to the solution directly derived from the action. The energy momentum tensor for non-vanishing tension is proportional to $`T_q`$. In the case of vanishing tension, it can be obtained from the action

$`S={\displaystyle \mathrm{d}^2\sigma \mathrm{\Lambda }\left(_aX^\mu ^aX^\nu G_{\mu \nu }\right)^2}.`$ (23)
This is the original source term (10) written with the help of a Lagrange multiplier and having performed the limit, i. e. $`T_q,_{\mu \nu }0`$. A better approach is to rewrite the action with the help of a world-sheet vector density $`V^a`$ of weight $`\frac{1}{2}`$

$`S={\displaystyle \mathrm{d}^2\sigma V^aV^b_aX^\mu _bX^\nu G_{\mu \nu }}.`$ (24)
The tensor quantity $`V^aV^b`$ has determinant zero and $`V^a`$ satisfies the following equations of motion:

$`V^b_aX^\mu _bX^\nu G_{\mu \nu }=0,_a\left(V^aV^b_bX^\mu \right)=0.`$ (25)
For simplicity, we restrict ourselves to the transverse gauge on the world-sheet by using the reparametrization symmetry, i. e. fixing

$`V^a=\left(\begin{array}{c}v\\ 0\end{array}\right).`$ (28)
Moreover, by using the diffeomorphisms on the world-sheet, we fix the constant $`v`$ to $`v=1`$. This is, in fact, the so called diffeomorphism gauge, yielding

$`T_{\mu \nu }(x^\sigma )=_\tau X_\mu _\tau X_\nu \delta ^{(8)}(X^\sigma x^\sigma ),`$ (29)
where the eight-dimensional delta function covers the space transverse to the string, cf. . Here, we already integrated out the world-sheet directions. This implies that $`X^u`$ and $`X^1`$ are fixed to the values of $`x^u`$ and $`x^1`$. Since $`X^2,\mathrm{},X^8,X^u=0`$, the only non-vanishing contribution arises from

$`X_u=G_{uv}X^v=X^v=\alpha \tau .`$ (30)
From this, we obtain

$`T_{uu}=_\tau X_u_\tau X_u\delta (\rho )\delta (u)=\alpha ^2\delta (\rho )\delta (u),`$ (31)
and thus, by comparing to (22),

$`\alpha ^2={\displaystyle \frac{1}{24}}\pi ^3\mathrm{\Delta }_{q,0}^{1/2}Q_0^{}.`$ (32)
Therefore, the momentum density of the (boosted) string is fixed by (32). The more important outcome is that (17) is ”compatible” with taking the limit directly inside the action.

## 4 Conclusions

We found an $`SL(2,)`$ family of string backgrounds describing tensionless strings by performing an infinite boost to the solution of Schwarz . The outcome is of pp-wave form incorporating a delta-like singularity along the position of the string. We showed that the anti-symmetric tensor field $``$ vanishes in this limit and hence do the 3-forms. The dilaton as well as the R-R scalar field become constant. We calculated the energy momentum tensor and found that is has the expected properties.

Hayashi and Samura performed a Penrose limit  on solutions of the type we found . Considerations similar to our were done by Horne et. al. . However, their solution did not contain a delta-like singularity and this was believed to be true for any boosted string.

As there is not much known about tensionless strings themselves, it is an interesting open question for further studies to look at the spectrum that arises from this new solution. On the other hand, tensionless strings arise naturally in the context of intersecting branes. Consequently, one might think of relating the presented results to tensionless strings arising in M-theory. One could also envisage connecting this to the context of string scattering from D-branes, cf. Garousi and Myers .

###### Acknowledgments.

UL acknowledges initial discussions with Bo Sundborg on the topic of this article. The authors further gratefully acknowledge discussions with James Gregory and Chris Hull. The research of UL is supported in part by EU contract HPNR-CT-2000-0122 and by VR grant 650-1998368.


text:
# Detection of Pulsed X-ray Emission from XMM-Newton Observations of PSR J0538+2817

## 1 INTRODUCTION

Of the $`>1000`$ radio pulsars detected so far, 51 have also been found in the X-rays Becker & Aschenbach (2002). These X-ray emitting pulsars represent a wide range of ages ($`10^37\times 10^9`$ yrs), magnetic field strengths ($`10^810^{13}`$ G), periods (1.6 – 530 ms) and spectral properties. In particular, while non-thermal (power law) emission is supposed to dominate in the youngest objects (like the Crab), the soft X-ray radiation of rotation-powered pulsars with ages as large as $`10^510^6`$ yrs should be dominated by thermal emission from the neutron star (NS) surface (see e.g. Becker & Pavlov 2002 for a recent review). These “middle-aged pulsars” are old enough for their magnetospheric emission to have decreased and become fainter than the thermal one. At the same time, they are young enough that the NS surface is not too cool for its thermal radiation to still be detectable in the X-rays.

There are three middle-aged pulsars for which a thermal component, interpreted as thermal emission from the surface, has certainly been observed: Geminga, PSR B0656+14 and B1055-52 (dubbed the three Musketeers by Becker & Trümper 1997). Due to its age, a dominant thermal component was also expected from PSR J0538+2817.

PSR J0538+2817 is a 143 ms pulsar located in the direction of the Galactic anti-center Anderson et al. (1994). The source was detected during an untargeted pulsar survey at 430 MHz using the 305 m Arecibo radio telescope Foster et al. (1995). PSR J0538+2817 is located at a distance of 1.2 kpc (using the model of Cordes & Lazio 2002), with a spin down age of $`6\times 10^5`$ yrs. The source is believed to be associated with the supernova remnant S147 Anderson et al. (1996) as the distance to the pulsar and the remnant are consistent, as are their ages. Sun et al. (1996) detected X-ray emission from PSR J0538+2817 using data from the ROSAT All Sky Survey (RASS). However, the photon statistics from the RASS were not good enough to perform timing or spectral studies. Marginal evidence was found for pulsed gamma-ray emission at the radio period from photon arrival time analysis of EGRET data Sun et al. (1996).

PSR J0538+2817 is a factor 5–10 more distant than Geminga, and only the large collective area of XMM-Newton has now made it possible to obtain good spectral information. In this paper we report on XMM-Newton observations of PSR J0538+2817 and present timing and spectral analyses of the source.

## 2 OBSERVATIONS

PSR J0538+2817 was observed with XMM-Newton for $``$18 ks on 2002 March 8 as a part of the Guaranteed Time Program. The following analysis uses data from the three European Photon Imaging Camera (EPIC) instruments: two EPIC MOS detectors Turner et al. (2001) and the EPIC PN detector Struder et al. (2001). The thin optical blocking filter was used on the PN. To provide a temporal resolution of 6 ms, which is needed for timing studies of millisecond pulsars, the PN was operated in small window mode. The MOS1 was operated in full window (imaging) mode with a time resolution of 1.4 s. In order to obtain better temporal resolution the MOS2 was operated in timing mode; in this mode data from the central CCD are collapsed into a one-dimensional row to achieve a 1.5 ms time resolution. The medium filter was used for both MOS observations.

Observations of PSR J0538+2817 were taken with the Resolution Grating Spectrometer (RGS), however we did not detect enough photons for a meaningful analysis. We also do not detect the pulsar in the images taken with the OM. We are unable to determine a realistic upper limit for the optical detection of PSR J0538+2817 as the position of the pulsar is coincident with the stray light effect in the OM image.

We reduced the EPIC data with the XMM-Newton Science Analysis System (SAS v5.3.0). To maximize the signal-to-noise we filtered the data to include only single, double, triple and quadruple photon events for MOS, and only single and double photon events for PN. We find no evidence for any extended emission in a spatial analysis of the MOS and PN data, indicating that PSR J0538+2817 is consistent with being a point source from XMM-Newton observations.

We extracted data for the source within a radius of 30$`\mathrm{}`$ for the PN and MOS1 instruments. While this radius encircles 100% of the flux for PSR J0538+2817 in the MOS1, only 89% of the flux is encompassed for the PN. However, due to the position of the source on the PN detector we are restricted in our choice of radii. Hence, the flux measured from the PN is corrected for in the following analyses. Selection of the photons from the MOS2 event file was achieved by extracting 100% of the flux within a rectangular region centered on the source.

In our temporal analysis we consider data within the 0.2–15 keV energy range for the EPIC MOS2 and PN instruments. The timing studies were performed on barycentrically corrected MOS2 and PN event files.

We used MOS1 and PN data to produce spectra of PSR J0538+2817. We subtracted a background which was extracted from an annulus around the source fiducial region. We also created the corresponding photon redistribution matrix (RMF) and ancillary region file (ARF).

## 3 TIMING ANALYSIS

We determine a predicted pulse period for PSR J0538+2817 from the radio measurements Anderson et al. (1996) at the epoch of our XMM-Newton observations (MJD 2452341.5), assuming a linear spin-down rate. It is possible however that glitches could alter the period from the predicted value. Therefore we searched for a pulsed signal in the PN ($``$6500 photons) and MOS2 ($``$3300 photons) data over a wider frequency range centered on the derived radio period, P = 143.158657 ms ($`f=6.9852569`$ Hz).

We have employed two methods in our search for pulsed emission from PSR J0538+2817. In the first method we implement the $`Z_n^2`$ test Buccheri et al. (1983). In the second method we calculate the Rayleigh Statistic (see de Jager 1991; Mardia 1972) and then calculate the Maximum Likelihood Periodogram (MLP) using the $`\mathrm{\Delta }C`$-statistic Cash (1979) to determine significant periodicities in the datasets (see Zane et al. 2002). The advantage of this method is that the $`\mathrm{\Delta }C`$-statistic makes no assumptions about the distribution of the data, and is itself distributed as $`\chi ^2`$ which allows one to determine errors for the periods detected directly from the periodogram.

We searched the PN data (0.2–15 keV) for a periodic signal using the $`Z_n^2`$ test, with the number of harmonics $`n`$ being varied from 1 to 10. We find that in order to optimize the signal-to-noise ($`S`$/$`N`$) ratio we must use the value $`n=1`$. The resulting $`Z_1^2`$ periodogram is shown in Fig. 1 (top, left panel). We find a peak at $`f=6.9852548_{0.0000124}^{+0.0000104}`$ Hz, which, within the quoted 90% uncertainty range, is consistent with the predicted radio frequency given above, $`f=6.9852569`$ Hz. The $`Z_1^2`$-statistic for this peak is 65.67, which has a probability of chance occurrence of $`5.5\times 10^{15}`$. Frequencies of this order are detected at the values expected, given radio measurements, in similar objects (J.A. Kennea 2003, private communication).

The corresponding $`Z_1^2`$ periodogram for the MOS2 data (0.2–15 keV) is shown in Fig. 1 (top, right panel). We note that the noise level is much greater in the MOS2 periodogram compared to that for the PN data. The largest peak near to the predicted radio frequency, and the second largest peak in the periodogram, is at $`f=6.9852939_{0.0000148}^{+0.0000117}`$ Hz. This value, within the 90% confidence limit quoted, is not consistent with the frequency found in the PN data, or the derived radio frequency. If we consider the 68% uncertainty range, $`f=6.9852939_{0.0000262}^{+0.0000243}`$ Hz, the frequency is still not consistent with the predicted and PN frequencies. The probability of chance occurrence for this peak is $`1.5\times 10^3`$. The largest peak in the MOS2 periodogram occurs at $`f=6.9939901_{0.0000133}^{+0.0000103}`$ Hz. However, it is difficult to reconcile this value with the predicted radio frequency. Both peaks have a low significance when compared with the noise level.

Assuming that the pulsed signal we have detected in the PN data is real we have investigated the energy range over which we can detect a periodic signal for the PN data. We separated the PN data into 4 energy bands, 0.2–1, 1–2, 2–3 and 3–10 keV, and generated $`Z_1^2`$ periodograms for each dataset. Our results show that we can detect a pulsed signal in the 0.2–1 and 1–2 keV ranges, but at energies greater than 2 keV we do not find a dominant peak in the periodograms near to the predicted radio frequency. We find a peak at $`f=6.9852548_{0.0000139}^{+0.0000095}`$ Hz for the 0.2–1 keV data (Fig. 1, bottom, left panel), and at $`f=6.9852499_{0.0000126}^{+0.0000127}`$ Hz for the 1–2 keV data (Fig. 1, bottom, right panel). Both of these frequencies are consistent with the frequency found previously for the PN data (0.2–15 keV), and the derived radio frequency, within the 90% uncertainty range quoted. The probability of chance occurrence for the 0.2–1 and 1–2 keV peaks are $`2.2\times 10^9`$ and $`3.0\times 10^8`$, respectively. The frequencies being consistent within different energy ranges is expected. The lower amplitude and significance of detection of pulsed emission in the separate energy ranges is simply due to the fact that fewer events are detected.

We also searched for a periodic signal using the MLP method. Initially we constructed the MLPs using the PN and MOS2 data in the 0.2–15 keV energy range to compare with our $`Z_1^2`$ results. We find a peak for the PN data at $`f=6.9852540_{0.0000039}^{+0.0000049}`$ Hz, which is detected at the 8-sigma level (Fig. 2, left panel). The frequency of the peak is consistent, within the 90% confidence limit, with the $`Z_1^2`$ value for the PN data, and with the predicted radio frequency.

The largest peak in the MOS2 MLP occurs at $`f=6.9852931_{0.0000185}^{+0.0000112}`$ Hz (Fig. 2, right panel). While this is not consistent with the PN and predicted radio frequencies, it is in agreement with the values found from the $`Z_1^2`$ test for the MOS2 data. While the peak could represent a true periodicity, its significance above the noise is low.

We used the MLP method on the PN data which had been separated into different energy ranges, given above. The results are in excellent agreement with the $`Z_1^2`$ test, again with the significance and amplitude of the detection decreasing with number of events.

While our results indicate that we have detected pulsed emission from PSR J0538+2817 in the PN data, the absence of a periodic signal at a similar frequency, within the 90% confidence limit, in the MOS2 data is worrying. However, we note that the MOS2 receives $`40`$% of the flux that is received by the PN, and that the quantum efficiency of the MOS instruments is poorer than that for the PN, especially in the softer energies where we are detecting the pulsed signal. Despite these factors, the most important consideration for this analysis is the $`S`$/$`N`$ ratio for the PN and MOS2 data. Our PN data has $`S`$/$`N6`$, while for the MOS2 we have $`S`$/$`N2.5`$. This suggests that the lack of signal in the MOS2 data is hampering the period search; our results for the MOS2 temporal analysis being consistent with noise is perhaps not surprising.

The frequency found from the $`Z_1^2`$ test for the PN data in the 0.2–15 and the 0.2–1 keV band, $`f=6.9852548`$ Hz, deviates from the derived radio frequency by 0.0000021 Hz, or 0.000043 ms. We have folded the PN light curves in the 0.2–15, 0.2–1, 1–2 and 2–10 keV bands using our detected frequency, and the radio ephemeris (MJD 2449564.5).

The 0.2–15 keV light curve (Fig. 3, top panel) has a pulse fraction of $`18\pm 3`$%. The pulse profile is clearly asymmetric with a slower rise to and faster fall from maximum. Indicative of this asymmetry, a sinusoidal fit to the light curve can be excluded with $`\chi _\nu ^2=4.0`$. The arrival time of the radio pulse is consistent with the maximum of the X-ray peak. The light curves in the 0.2–1 and 1–2 keV bands (Fig. 3, second and third panels) have pulse fractions of $`30\pm 3`$% and $`27\pm 3`$%, respectively.

In order to investigate the asymmetry in the 0.2–15 keV light curve we determined whether or not there is a phase shift between the 0.2–1 and 1–2 keV light curves. Although the PN light curve in the 0.2–15 keV energy range was not well described by a sinusoid, as a first estimate to the maxima of the 0.2–1 and 1–2 keV light curves we have employed the method of fitting with a sinusoid, as these two light curves are more symmetrical. While the fit for the 1–2 keV light curve is improved by a factor of $``$ two the fit is still not optimal, and the value of $`\chi _\nu ^2=3.6`$ for the 0.2–1 keV light curve still indicates an asymmetry. However, for estimation purposes we have marked the position of the maxima of the cosine curves for each band in Fig. 3. The maximum in the 0.2–1 keV band occurs at $`\varphi =0.94\pm 0.05`$, and at $`\varphi =0.86\pm 0.05`$ in the 1–2 keV band. While this indicates that the maxima in the two bands are shifted by 0.08 in phase, corresponding to a pulse shift of 29, the uncertainty on the peak positions does not rule out the null hypothesis that there is zero phase shift. We note that this shift in phase is small compared to the pulse width. If the shift is real, this could be the origin of the asymmetry in the 0.2–15 keV light curve.

In an attempt to improve the sparse information we have from the sinusoidal fit to the PN light curve in the 0.2–15 keV band, we have also modeled the brightness variation with a symmetrical polar cap model Cropper et al. (2001), which includes gravitational bending effects. The input parameters for the model are the angular size of the caps, the angle between the magnetic and rotation axis and that between the viewing direction and the rotation axis. The best fit to the data occurs for 20, 11 and 31 for the three angles respectively, with a $`\chi _\nu ^2`$ of 1.71. The compactness parameter $`M/2R`$ was set to a typical value of 3. The $`\chi _\nu ^2`$ for the polar cap model is significantly better than for the sinusoidal fit with the same number of free parameters, but falls short of that required because of the symmetry of the model prediction.

The absence of any sharp features in the light curve, and the absence of any asymmetry in our model results in degeneracies between the fit parameters and render the formal uncertainties on their values to be large. By exploring the parameter space, we find that a range of viewing directions between 10 and 60 is permitted, with the angle between magnetic and rotation axis constrained to between 20 and 35 over this range. The change in viewing direction is in fact accommodated by a change in cap size to produce the observed pulsation, with the semi-angle subtended increasing from a negligible size to 30 as the viewing angle increases from 10 to 60. This analysis is therefore limited given the systematic differences between the model and the data. Nevertheless, the general point that we note is that the method allows us to determine an upper limit for the emitting region. The pole cap size has to be constrained to a full angle of $``$65, otherwise the observed pulsation amplitude cannot be reproduced.

Another factor which could be responsible for the asymmetry is aberration Cordes (1978). This effect has been seen in radio observations of pulsars Gangadhara & Gupta (2001) and is due to corotation which causes the radiation beams to bend toward the azimuthal direction (the direction of rotation of the pulsar) so that emission is received earlier than if there were no rotation. The components emitted at higher altitudes in the light cone are shifted to an earlier phase due to this aberration, resulting in an asymmetry in the location of the leading and trailing components with respect to the center of the pulse profile. For the same cone, the lower frequency radiation comes from a higher altitude, and hence is observed first, than the higher frequency radiation. The effects of aberration increase with altitude. If we assume that the emission we observe from PSR J0538+2817 originates from close to the surface of the neutron star i.e. a heated polar cap, and not at a height which is near to that of the light cylinder, the effects of aberration should be negligible.

## 4 SPECTROSCOPY

Data from MOS1 and PN were used to create spectra of PSR J0538+2817. In order to preserve any features in the spectra while providing sufficient signal-to-noise, the spectra were re-binned to ensure a minimum of 15 counts in each energy bin. We fit the combined MOS1 and PN spectra with the X-ray spectral analysis package XSPEC (v11.1.0) using various models, each with interstellar absorption ($`N_H`$) accounted for by XSPEC’s photoelectric absorption model. The spectra were fit between 0.2–3 keV.

We used three different methods to determine the value of $`N_H`$ to be used in the model fits: (1) $`N_H`$ was allowed to run free, (2) $`N_H`$ was fixed at the estimated fraction of galactic absorption in the direction and at the distance of the pulsar with a value of 1.3$`\times 10^{21}`$ cm-2 ***This was achieved by determining the galactic absorption in the direction of the pulsar Dickey & Lockman (1990), and calculating the fraction of this column depth at the distance of the pulsar., (3) $`N_H`$ was fixed at the value of 2.5$`\times 10^{21}`$ cm-2 as measured from a star near to the pulsar in position and distance.

PSR J0538+2817 is classed as a “middle-aged” pulsar and is thought to perhaps have similar spectral properties to the three Musketeers Becker & Trümper (1997); Becker & Pavlov (2002). The spectra of PSR B0656+14 Becker & Pavlov (2002) and B1055-52 Pavlov et al. (2002) are similar and can be well-fit with a model which comprises of two blackbodies and a power law (PL), known as the TS+TH+PL model. The current interpretation is that the PL component is due to some non-thermal magnetospheric emission, while the two blackbodies represent a soft thermal component (TS) emerging from most of the NS surface and a hard thermal component (TH) which originates from heated polar caps. Geminga is not bright enough to allow one to discriminate between the different components of the temperature distribution at the surface: its thermal spectrum is well-fit by a single blackbody plus power law.

In order to study the spectra of PSR J0538+2817, and to investigate the pulsar’s similarity to the three Musketeers we fit the PN and MOS1 spectra with a variety of models. These were as follows, a single blackbody, a single power law, a blackbody plus power law, a two blackbody, and the TS+TH+PL model.

We find that the multi-temperature models yield the same value for the temperature for both components indicating that only one blackbody is required. The large uncertainties on the power law index for the single power law, blackbody plus power law and TS+TH+PL fit indicate that the PL component is not required.

Our results indicate that the spectra can be well-fit with a single blackbody (Fig. 4, top panel). Table 1 shows the results for the fits with different $`N_H`$ values as described above. Comparing the $`\chi _\nu ^2`$ values for the three fits, and the values of $`N_H`$, suggests that the estimate of the fraction of galactic absorption in the direction of PSR J0538+2817 is too low (fit (2)).

From our model fit we find $`N_H`$ = ($`2.63_{0.02}^{+0.02}`$)$`\times 10^{21}`$ cm-2 and $`T^{\mathrm{}}=(2.09_{0.07}^{+0.06})\times 10^6`$ K for fit (1) in which $`N_H`$ is allowed to run free. The fit results in $`\chi _\nu ^2`$ = 1.08 for 320 d.o.f. For fit (3) in which $`N_H`$ is fixed at the measured value for a star at a similar distance and position as the pulsar of $`N_H`$ = 2.51$`\times 10^{21}`$ cm-2 gives $`T^{\mathrm{}}=(2.12_{0.03}^{+0.04})\times 10^6`$ K, with $`\chi _\nu ^2`$ = 1.08 for 321 d.o.f. Employing the F-test we find a probability of 0.4 which indicates that freeing $`N_H`$ is not warranted. In Fig. 5 we show the error ellipses in the temperature-absorbing column plane which confirm this result.

We find unabsorbed luminosities of 4.0$`\times 10^{32}`$ ergs s-1 and 2.7$`\times 10^{32}`$ ergs s-1 from XSPEC in the ROSAT 0.1–2.4 and 0.5–10 keV bands, and a blackbody luminosity of 4.6$`\times 10^{32}`$ ergs s-1 for fit (3). The emitting radius resulting from the single blackbody fit is $`R^{\mathrm{}}=1.68\pm 0.05`$ km. This indicates that the emission originates from a region which is smaller than the surface of the NS.

However, the radiation radius is representative of the stellar radius only when the spectral properties of the emergent radiation are consistently accounted for. To investigate whether one can increase the emitting radius so that it represents the whole of the neutron star we have fit the data with a pure-H, non-magnetized atmospheric model Pavlov et al. (1992); Zavlin et al. (1996) which is publicly available. We note that the observed modulation in the light curve argues against the emission being generated by the full surface of the NS as implied by an atmospheric model, however the size of the emitting region determined from the model is dependent on the value of the normalization constant and can therefore represent a smaller region. Typically, these models are found to require lower color temperatures and larger radii with respect to a simple blackbody. The results of the fits are given in Table 2. In the fits we let $`N_H`$ run free, but we fix the mass, $`M_{NS}`$, and radius, $`R_{NS}`$, of the NS. For $`M_{NS}`$ = 1.4 M$`\mathrm{}`$ and $`R_{NS}`$ = 10 km we find $`N_H`$ = ($`4.18_{0.06}^{+0.01}`$)$`\times 10^{21}`$ cm-2 and $`T_{eff}`$ = ($`0.68_{0.01}^{+0.10}`$)$`\times 10^6`$ K, with $`\chi _\nu ^2`$ = 1.09 for 320 d.o.f (Fig. 4, bottom panel). For $`M_{NS}`$ = 1.4 M$`\mathrm{}`$ and $`R_{NS}`$ = 18 km we find $`N_H`$ = ($`4.17_{0.05}^{+0.01}`$)$`\times 10^{21}`$ cm-2 and $`T_{eff}`$ = ($`0.57_{0.01}^{+0.06}`$)$`\times 10^6`$ K, with $`\chi _\nu ^2`$ = 1.10 for 320 d.o.f.

The atmospheric model gives lower temperatures than the blackbody model which are more in agreement with the predicted values from standard cooling curves (see e.g. Tsuruta 1998; Page 1998; Slane, Helfand & Murray 2002; see also Section 5). The standard cooling model is based on neutrino emission via the modified URCA process, n-n and p-p neutrino bremsstrahlung, crust neutrino bremsstrahlung, and plasmon neutrino processes. Assumptions about the form of the nucleon-nucleon force and the many-body technique used to calculate the equation of state must be made in the cooling models. The standard curve considered here (see Fig. 6, solid line) uses an equation of state of moderate stiffness. However, the distance to the pulsar derived from the atmospheric model fits are too low. We find that the fit with $`R_{NS}`$ = 10 km gives a pulsar distance of 0.26 kpc, and 0.38 kpc with $`R_{NS}`$ = 18 km (see Section 1). While this suggests that the atmospheric model can be ruled out, we note that the model we have used is non-magnetized, and with $`B7\times 10^{11}`$ G for PSR J0538+2817 the applicability of this model is questionable. Further studies with a magnetic model must be performed to address this.

We note that our spectral fits have been obtained from data which are averaged over the spin pulse. However, the pulsation can in principle be caused by a local variation in the emission properties of the NS surface, or by some absorption in a frame rotating with the NS. These effects may have implications for any spectral fitting to a mean spectrum.

## 5 DISCUSSION

### 5.1 Results from the Spectral Analysis

We have analyzed spectra of PSR J0538+2817 from the MOS1 and PN instruments on XMM-Newton. We find that the spectra are well described by a single blackbody with $`T^{\mathrm{}}=2.12\times 10^6`$ K and $`N_H`$ = 2.5$`\times 10^{21}`$ cm-2. The temperature inferred from the model fit results in an emitting radius of $`R^{\mathrm{}}=1.68`$ km. It is difficult to reconcile this radius with that of the whole neutron star. The higher temperature implied by the blackbody fit could indicate that there is a heating mechanism occurring at the poles i.e. magnetospheric processes, inverse Compton scattered photons above the polar cap producing electron-positron pairs resulting in surface heating by the returning positrons (see Greiveldinger et al. 1996; Harding & Muslimov 2002).

An inferred temperature from a blackbody fit which results in a radius which is smaller than the surface of the neutron star has been noticed in several other sources Becker & Pavlov (2002). In general, fits with un-magnetized, pure-H atmospheric cooling models give a temperature lower by a factor 1.5–3, but often they also give radii that are too large (or, equivalently, distances which are too small). Therefore, caution must be taken in whether they are directly applicable. The same is true in the assumption that the thermal emission is due to accretion from the interstellar gas, since pure-H cooling and accretion spectra are almost indistinguishable in the X-rays Treves et al. (2002). If the atmosphere is dominated by heavy-elements, the fit values from the model spectra are close to those from a blackbody, but in this case one also expects absorption lines and edges Becker & Pavlov (2002).

In the case of PSR J0538+2817, we have found that a pure-H, un-magnetized atmospheric model also provides an excellent fit to the spectral data. The resulting temperature is a factor $`3`$ lower than that derived from the blackbody fit. The value is therefore high, but certainly in better agreement with that predicted by the cooling models (see Fig. 6). However, these fits imply either very small distances to the pulsar ($`200400`$ pc) or unacceptably large NS radii, so the applicability of these spectral models to PSR J0538+2817 is questionable. A similar conclusion has been reached by Pavlov et al. (2002) in the case of PSR 0656+14, based on Chandra data. We note that in the case of PSR J0538+2817 fitting with a magnetized must be performed to be able to rule out conclusively an atmospheric model.

PSR J0538+2817 has often been claimed to be potentially similar to other middle-aged sources, PSR B0656+14, B1055-52 and Geminga (the three Musketeers, see Table 3). As noted in Section 4 the power law component is detected in the spectra of the three Musketeers. Results from ROSAT indicate that all of the three objects exhibit a pulse shift of $`100^{}`$ at $`E0.40.6`$ keV, and an increase in the pulse fraction (confirming the different nature of the thermal and non-thermal emission). However, recent results for B1055-52 deduced from XMM-Newton observations indicate that the pulse phase energy dependence is not as strong as suggested by ROSAT. The surface temperatures of the three Musketeers (as taken by the blackbody fit) are in the range (0.3–1.2)$`\times 10^6`$ K which is typically too high to be reconciled with that predicted by standard cooling curves Becker & Pavlov (2002). We note that the temperature determined for PSR J0538+2817 from the blackbody fit is even higher. The inferred radii for the three Musketeers are compatible with the equation of state of a neutron star, or smaller. However, as discussed above this issue is extremely delicate since the values of peak temperature, star radius (and distance) are strongly dependent on the assumption about the spectral model.

In the light of this situation, PSR J0538+2817 is an extremely interesting object. Its spectrum only shows a thermal component: the non-magnetospheric power law has not been detected by XMM-Newton. This result must be confirmed by further deep pointings: for instance, we note that in the case of B1055-52 the ROSAT spectrum was well-fit by either two blackbodies Greiveldinger et al. (1996) or a flatter power law Wang et al. (1998). Only very recently have Chandra observations of B1055-52 Pavlov et al. (2002) revealed the spectral energy distribution in detail; the Chandra spectrum is well-fit by the TS+TH+PL model. We note that B1055-52 is at a distance of $`0.7`$ kpc and faint, PSR J0538+2817 is also faint and at a larger distance, hence both sources are difficult to observe.

On the other hand, if confirmed, the absence, or very low level, of a hard tail in the spectrum of PSR J0538+2817 is peculiar. While PSR J0538+2817 and B1055-52 have similar ages, PSR J0538+2817 may represent a transition object between those discussed above (PSR B0656+14, B1055-52 and Geminga) and older ($`>10^6`$ yrs) sources (e.g. B0950+08, B1929+10, J0437-4715), based on its spectral properties. The latter have a NS surface (excluding the cap) which is much cooler and therefore not detectable in the X-ray range. However, the polar caps may be heated by accretion or bombardment of particles thus producing the observed X-rays. In this case, the soft X-ray spectrum may be dominated by the TH component discussed above. This interpretation may also help to reconcile the temperature inferred from the spectral fit of PSR J0538+2817 which is higher than the predicted surface temperature of a NS at the age of PSR J0538+2817 (and makes the comparison with the cooling curves unfounded).

If the temperature from our blackbody fit to the spectrum of PSR J0538+2817 is taken as the TH component, we find an upper limit, from the 90% uncertainty level, of 0.56$`\times 10^6`$ K for the surface temperature (TS) of the neutron star. We find that the 2$`\sigma `$ upper limit for the unabsorbed flux in the 3.0–10.0 keV band for a power law with a photon index of 2.0 is 5.9$`\times 10^{14}`$ erg cm-2 s-1.

### 5.2 The Pulse Profile

The XMM-Newton observation of PSR J0538+2817 provides the first evidence for pulsed X-ray emission at a period which is consistent, within errors, to the predicted radio pulse period. The pulse profile is broad and slightly asymmetric, with a slower rise and more rapid decline.

Pulsed X-rays from PSR J0538+2817 could be produced via one of two mechanisms. The first is due to relativistic particles produced either in regions above the polar caps of the neutron star or in outer gaps. These particles are accelerated by the magnetic field which is aligned with the electric field and/or the drop in potential generated by the Deutsch wave. Non-thermal X-ray emission results from the effects of curvature, synchrotron and inverse Compton radiation on the accelerated particles. The other mechanism is thermal emission which may be produced at the stellar surface due to the initial cooling of the neutron star, or reheating of the polar caps by the bombardment of back-flowing particles and internal friction.

The particular X-ray emission mechanism can sometimes be determined from the pulse shape. High amplitude, sharp, narrow pulses can only be produced by a highly beamed, and therefore relativistic population of electrons, while low amplitude, quasi-sinusoidal pulses can be produced by either method. Since the latter characteristics are those seen in the PN light curve of PSR J0538+2817, in principle we cannot discriminate between the two mechanisms. However, the absence of a power law tail in the spectrum, argues toward a thermal origin.

The observed pulsation may also be due to absorption in the rotating magnetosphere, instead of being an intrinsic characteristic of the pulsar emission (see Cropper et al. 2002 for a discussion of a similar effect in the case of RX J0720–3125). However, the thermal emission of PSR J0538+2817 is hotter than that of RX J0720–3125, and this implies a higher ionization parameter (by a factor of $`20`$). Moreover PSR J0538+2817 has a faster rotation period (by a factor of $`>`$ 50) which limits the maximum extent of the light cylinder and therefore the volume of absorbing material. The high densities of largely neutral material required, together with the rather weak magnetic field in this system to confine it, argue against absorption as the origin of the pulsations.

The simplest explanation for the flux variation is therefore that it is due to hot, rotating polar caps. The asymmetry in the pulse shape in Fig. 3 (top panel) cannot be produced by any symmetrical structure on the NS, nor by any beaming effects, as long as these operate symmetrically with respect to the magnetic or rotation axis. On the other hand, if the emission originates in a hot spot (or cap), the observed asymmetry can be qualitatively understood by requiring that the trailing edge of the cap is hotter and brighter than the leading edge. The DC component of the pulse profile implies that there is never a time when the emission is “off” (see Fig. 3). If the polar cap model is correct this suggests that part of the cap is always visible.

## 6 CONCLUSIONS

We find that the spectra of PSR J0538+2817 can be well-fit with a blackbody. The value for the emitting radius indicates that the emission is from a region which is smaller than the surface of the NS, i.e. a hot spot. We investigated whether a pure-H atmospheric model would provide model temperatures to be in line with that expected from the cooling curve. We find that the model fits the spectra well, but the distance to the pulsar determined from the fit is too small, unless we use unrealistically large radii for the NS. We do not find any evidence for a non-thermal component in the spectra of PSR J0538+2817. However, we note that the power law component was detected in B10522-52 only after deeper observations. This may also be true for PSR J0538+2817.

We have detected for the first time pulsed X-ray emission from PSR J0538+2817, the period of which is consistent, within errors, with the predicted radio period. Whether the mechanism for X-ray emission is via non-thermal magnetospheric synchrotron or thermal emission from the hot surface is not clear, as the broad asymmetric pulse shape found for PSR J0538+2817 does not favor either model. However, the absence of a power law component may indicate that the emission is thermal in origin.

We find an asymmetric pulse profile which may be due to hot, rotating polar caps, with a temperature structure which is hotter at the trailing edge. If the polar cap model is correct, the small radius we determine for the emission region is expected.

## 7 ACKNOWLEDGMENTS

This work is based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. The authors acknowledge support from the Institute of Geophysics and Planetary Physics (IGPP) program at LANL and NASA grants S-13776-G and NAG5-7714. K. McGowan acknowledges S. Trudolyubov, and S. Zane acknowledges R. Turolla, for helpful discussions.


text:
# A modified SPH approach for fluids with large density differences

## 1 Motivation

SPH is a Lagrangian particle method for solving the equations of hydrodynamics that was invented by Lucy  and Gingold and Monaghan . Instead of discretising space with a grid, the matter is discretised into so-called particles which move with the fluid and do not exchange mass. This method is especially suited for compressible flows with irregular boundaries. SPH has been used for many astrophysical problems with great success, and has also been applied to other fields of physics, such as e.g. the simulation of liquids  and solids .

Due to its Lagrangian nature, simulating several non-mixing fluids is a straightforward extension to SPH. Each particle gets initially marked with the phase it belongs to, and these marks do not change with time.

We are interested in simulating compressible flows with large density differences in contact discontinuities, e.g. due to an interface between a jet of liquid near the speed of sound in a surrounding gas. Unfortunately, standard SPH as e.g. presented in  breaks down in this case, as explained in this text. Previous approaches to handle large density differences in SPH (e.g. ) rely on ad-hoc countermeasures.

For subsonic flows,  presents a mesh-free numerical method for gas-liquid phase interfaces which is based on the MPS method . It is restricted to incompressible flows, whereas we are interested in compressible flows. We are also interested in similar spatial resolutions in both phases of the flow. This requirement is different from e.g. astrophysical collapse simulations, where one usually wants a higher resolution in the denser regions.

Another point of interest for us is having differing spatial resolutions, i.e. particles with different masses, in regions where the mass density does not change. This has been used on several occasions (e.g. in ). However, this can lead to rather large errors with standard SPH.

We will in the following present a modification of SPH that interprets certain numerical quantities in a different manner , leading to a stable, robust, and accurate evolution even in these cases.

## 2 Describing multi-phase flows

One commonly used way of introducing the SPH discretisation (see also ) starts out by considering an arbitrary field $`f(𝐱)`$. This field is first *smoothed* by folding it with a kernel $`W(𝐱)`$, which leads to the smoothed field $`f(𝐱)`$

$`f(𝐱)`$ $`:=`$ $`{\displaystyle d^3x^{}f(𝐱^{})W(𝐱𝐱^{})}`$ (1)
where the kernel $`W(𝐱)`$ must be normalised according to $`d^3xW(𝐱)=1`$. One usually chooses kernels that have approximately the shape of a Gaussian, and that have compact support for reasons of efficiency. The size of the domain of support is called the *smoothing length* and is usually denoted with the letter $`h`$.

In the next step, the smoothing integral is discretised at $`N`$ particle positions $`𝐱_i`$, which can in principle be chosen freely, but should of course be “reasonably” distributed. This leads to the *SPH approximation* $`\stackrel{~}{f}(𝐱)`$

$`\stackrel{~}{f}(𝐱)`$ $`:=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}V_jf_jW(𝐱𝐱_j)`$ (2)
of the field. The volumes $`V_j`$ are the discrete counterparts of the volume element $`d^3x^{}`$ in the integral above, and have to be chosen so as to be consistent with the spatial distribution of the particles. It is customary to assign a certain mass $`m_i`$ to each particle, and then replace the volumes through

$`V_i`$ $`=`$ $`{\displaystyle \frac{m_i}{\rho _i}}`$ (3)
where $`\rho _i`$ is the discretised mass density assigned to the particle. This is motivated by the fact that the particles do not exchange mass, which leads to the time evolution equation

$`{\displaystyle \frac{d}{dt}}m_i`$ $`=`$ $`0`$ (4)
making $`m_i`$ a natural choice for one of the primary variables. In order to make SPH Lagrangian, the particles have to move with the flow, leading to

$`{\displaystyle \frac{d}{dt}}𝐱_i`$ $`=`$ $`𝐯_i`$ (5)
as the time evolution equation for the particle positions $`𝐱_i`$. Here $`𝐯_i`$ is the discretised fluid velocity field.

Assuming that the desired spatial resolution is about the same in all fluid phases, one would choose similar particle spacings there, leading to similar particle volumes $`V_i`$. If the mass densities in the different fluids are about the same, no further problems arise, and SPH as usual can be used to describe them. However, if the difference in mass density is large (say, about one order of magnitude or more), then the particle masses will differ by the same factor, leading to problems at the phase interfaces. These problems are mostly caused by inaccuracies in the mass density. These inaccuracies are substantial for largely different particle masses. Below we describe how they come about, and how they can be avoided.

### 2.1 Standard SPH

There exist in principle two different methods for obtaining the discretised mass density $`\rho _i`$ in the “standard” SPH formalism. Both start out by considering the approximate mass density $`\stackrel{~}{\rho }`$ at the particle position $`𝐱_i`$. This quantity is obtained from eqn. (2) by choosing $`f(𝐱)=\rho (𝐱_i)`$ and applying eqn. (3), leading to

$`\stackrel{~}{\rho }(𝐱_i)`$ $`=`$ $`{\displaystyle \underset{j}{}}m_jW_{ij}`$ (6)
where the abbreviation $`W_{ij}:=W(𝐱_i𝐱_j)`$ has been introduced.

The first and conceptually simpler method to obtain $`\rho _i`$ from this is by setting $`\rho _i:=\stackrel{~}{\rho }(𝐱_i)`$, leading to

$`\rho _i`$ $`=`$ $`{\displaystyle \underset{j}{}}m_jW_{ij}\text{.}`$ (7)
This method is often used for astrophysical problems when there are free boundaries, i.e. when the matter distribution extends into vacuum. However, it is not suited for a phase interface with a large density discontinuity. The smoothing inherent in eqn. (7) will smooth out the density jump over a region of the size of the smoothing length $`h`$ in either direction of the discontinuity. Particles in this region will then “see” a density that is much less or much larger from the real density in their phase. When this density is used to calculate the pressure through the two phases’ equations of state, the pressure will be very wrong (as shown in figure 1 for two ideal gases), and it is basically impossible to set up a stable interface in equilibrium. This problem becomes even more severe when one of the fluids has a stiff equation of state, as is the case e.g. in liquids or solids, because then density inaccuracies will lead to even larger errors in the pressure.

The second method for obtaining the mass density $`\rho _i`$ is by integrating $`\rho _i`$ in time via the time derivative of eqn. (7), leading to

$`{\displaystyle \frac{d}{dt}}\rho _i`$ $`=`$ $`{\displaystyle \underset{j}{}}m_j(𝐯_i𝐯_j)W_{ij}`$ (8)
where eqn. (5) has been used, and the abbreviation $`W_{ij}:=(W)(𝐱_i𝐱_j)`$ has been introduced. This method has the advantage that the initial data for $`\rho _i`$ can be chosen freely, so that density discontinuities can be modelled. This can be used to simulate surfaces of liquids and solids.

The problem with this method is similar to the problem encountered when calculating the density directly from the particle distribution by eqn. (7). The particles on each side of a phase interface “see”, via the term $`m_j`$ in eqn. (8), very different particle masses on the other side of the interface. The values of $`d/dt\rho _i`$ then contain large inaccuracies, leading to instabilities at phase interfaces.

### 2.2 Modified SPH

However, eqn. (7) is not cast in stone. An ansatz equivalent to but different from the one leading to this equation is not to smooth the mass density $`\rho _i`$, but rather the particle density $`n_i=1/V_i`$ . This is easily motivated by the fact that the mass density can be discontinuous over a phase interface, while the particle density is not, according to our assumption of similar spatial resolutions on both sides. Smoothing the particle density $`n_i:=\stackrel{~}{n}(𝐱_i)`$ via eqn. (2) leads to

$`n_i`$ $`=`$ $`{\displaystyle \underset{j}{}}W_{ij}`$ (9)
and by taking its time derivative, the equation

$`{\displaystyle \frac{d}{dt}}n_i`$ $`=`$ $`{\displaystyle \underset{j}{}}(𝐯_i𝐯_j)W_{ij}`$ (10)
is obtained after using eqn. (5). As it is customary in SPH to use $`\rho _i`$ instead of $`n_i`$, we apply eqns. (3) and (4) and arrive at

$`{\displaystyle \frac{d}{dt}}\rho _i`$ $`=`$ $`m_i{\displaystyle \underset{j}{}}(𝐯_i𝐯_j)W_{ij}\text{.}`$ (11)
This new formulation of the equation of continuity is the key element of our SPH approach. It should be noted that this equation is identical to eqn. (8) when all particle masses $`m_i`$ are the same, which is the case for many single-phase SPH simulations.

For the simulations presented in this text, we discretise the Euler and the internal energy equations in established ways:

$`{\displaystyle \frac{d}{dt}}𝐯_i`$ $`=`$ $`{\displaystyle \frac{1}{\rho _i}}{\displaystyle \underset{j}{}}{\displaystyle \frac{m_j}{\rho _j}}\left(p_j+p_i\right)W_{ij}`$ (12)
$`{\displaystyle \frac{d}{dt}}e_i`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{1}{\rho _i}}{\displaystyle \underset{j}{}}{\displaystyle \frac{m_j}{\rho _j}}\left(p_j+p_i\right)(𝐯_i𝐯_j)W_{ij}`$ (13)
where $`e_i`$ is the specific internal energy. The symmetrisations $`(p_j+p_i)`$ are e.g. explained in .

## 3 Tests

In the following we test the new SPH formulation and compare it to analytic solutions as well as simulations performed using standard SPH as described in . That is, the only difference between these two formulations is that we use eqn. (11) instead of eqn. (8). This also means that other parts of an SPH code such as time integration or artificial viscosity are not affected. As test cases, we use an advection problem, a sound wave encountering a discontinuous change in the sound speed, and a shock tube with a Diesel–air interface.

### 3.1 Advection equation

We compare the standard and the new SPH formulation by simulating a one-dimentional advection equation. That is, we solve the equation of continuity for the density $`\rho `$ while prescribing the velocity field $`v`$. The velocity profile (which is constant in time) and the initial density profile are given by

$`v(x)`$ $`=`$ $`{\displaystyle \frac{x}{1+qx^2}}`$ (14)
$`\rho _0(x)`$ $`=`$ $`Ax^2\mathrm{exp}\left\{\left({\displaystyle \frac{xx_0}{W}}\right)^2\right\}`$ (15)
with $`A=1.5`$, $`x_0=1`$, $`W=0.4`$, and $`q=0.2`$. We initially place the particles with equidistant spacings with a density of $`n=10`$ particles per unit length and use a smoothing length of $`h=0.25`$. The particle masses are chosen according to the initial density profile at the initial particle positions, i.e. they differ. Advection problems are particularly well suited test problems for Lagrangian methods, so we expect a high accuracy from this low resolution.

The results of simulating this equation with both the standard and the new SPH formulation are presented in figure 2, which shows the density $`\rho `$ at five different times. Both formulations track the analytic solution very nicely in spite of the coarse resolution. However, at later times, the standard SPH formulation underestimates the density near the peaks, while the new formulation stays much closer to the analytic solution.

### 3.2 Sound wave

The sound wave test case consists of two regions containing the same ideal gas, but with different densities and in pressure equilibrium. The density discontinuity is located at $`x=0`$ with a density ratio of $`10:1`$. These conditions also lead to different temperatures, and to sound speeds with a ratio of $`1:\sqrt{10}`$. Figure 3 shows an initially Gaussian-shaped sound wave at different times.

At $`t=0.5`$ the initial wave travels to the left. At $`t=0`$ the wave has reached the interface where it is partially transmitted and partially reflected. At $`t=+0.5`$ the wave consists of two packets, travelling in different directions with different speeds. The simulation was performed with $`n=200`$ particles per unit length and a smoothing length of $`h=0.1`$. The analytic solution is shown as dotted line underneath the simulation result. The SPH simulation with the new equation of continuity (11) tracks the analytic solution quite well. On the other hand, standard SPH using eqn. (8) performs rather poorly in this case, as can be seen in the graph in the lower right hand corner: the transmission and reflection coefficients are wrong, and the pressure develops spikes at the interface. We assume that the reason for this is just the one demonstrated in figure 1.

### 3.3 Shock tube

A further test case for our formulation is a shock tube containing a Diesel–air interface, shown in figure 4. The initial discontinuity is at $`x=0`$, with air to the left and (liquid) Diesel oil to the right. The initial pressure ratio is $`2:1`$, the density ratio about $`1:80`$. Here we have artificially decreased the air density and increased the air temperature by a factor of $`10`$ to create a more difficult test case. Table 1 lists the exact initial data for this test case. The shock wave in the Diesel oil travels to the right, the rarefaction wave in the air to the left. Because the Diesel oil is nearly incompressible, the final pressure is close to the initial pressure in the air phase. The equation of state for the Diesel oil was kindly provided to us by the Robert Bosch GmbH.

The simulation was performed with $`n=100`$ particles per metre and a smoothing length of $`h=5\times 10^2\mathrm{m}`$. We use the artificial viscosity presented in  with a viscosity coefficient of $`\alpha =0.5`$, because some artificial viscosity is necessary to produce entropy in the shock front. This spreads out the shock front over about 8 smoothing lengths, which is acceptable for our purposes.

The initial pressure discontinuity remains visible as spikes at the contact discontinuity. These spikes are caused by the numerical initial data, which have a discontinuity and hence contain high frequency modes that are not resolved in the simulation. According to eqn. (1), the initial data should be smoothed before the SPH formalism is applied. We skip this step because we want to show that these high frequency modes do not harm the simulations. They remain present, but are not amplified. The formulation is stable. Table 2 compares several important quantities of the simulation result to the analytic solution, showing very good agreement in the shock relations.

We did not manage to perform this simulation with standard SPH. As illustrated in figure 1, the error in the density near the contact discontinuity leads to an error in the pressure. The pressure error can be orders of magnitude larger than the pressure itself for fluids which have a stiff equation of state. This problem does not exist in the new formulation.

## 4 Conclusion

We describe a modification to the standard SPH formalism that smoothes the particle density instead of the mass density. As tested by simulating an advection equation, sound waves, and shock waves, the new formulation either yields more accurate results than standard SPH, or the equivalent simulation with standard SPH is not stable. We conclude that this modified SPH formulation is an effective method for simulating multi-phase flows, and conjecture that this modification is beneficial for all simulations where particles have different masses.

The authors would like to thank Prof. Hanns Ruder for encouraging this work. We would also like to thank the Robert Bosch GmbH for providing an equation of state for Diesel oil. Financial support was provided by the Robert Bosch GmbH and the Ministerium für Wissenschaft und Forschung in Baden-Württemberg.