text:
# DO WE REALLY OBSERVE A BOW SHOCK IN N157B…?

## A PULSAR WIND INSIDE A SUPERNOVA REMNANT

A supernova explosion marks the end of the evolutionary track of a massive star, but launches the beginning of the evolution of a supernova remnant (SNR). If the fossil of the progenitor star is an active pulsar, the associated pulsar wind drives a pulsar wind nebula (PWN) into the surrounding medium (Rees and Gunn 1974). The evolution of the PWN is coupled and determined by the evolution of the SNR, because the total energy release of the relativistic pulsar wind over the pulsar’s lifetime is small ($`10^{49}10^{50}`$ erg) compared with the total mechanical energy of the SNR ($`10^{51}`$ erg).

A young SNR ($`100\mathrm{1\hspace{0.33em}000}`$ years) consists of a twofold shock structure, i.e. a forward shock which propagates into the interstellar medium (ISM) and a reverse shock which forms due to the deceleration of the forward shock by the surrounding ISM (McKee 1974). When the forward shock has swept up a comparable amount of mass as was ejected in the explosion event, the reverse shock propagates back into the interior of the SNR, where the PWN is expanding through the freely expanding ejecta of the progenitor star. As the reverse shock encounters the PWN, the latter is crushed, due to the huge pressure downstream of the reverse shock relative to the pressure inside the PWN (van der Swaluw et al. (2001); Blondin et al. (2001)).

Van der Swaluw et al. (2001) have shown that in the case of a centered pulsar, the interaction between the reverse shock and the centrally powered PWN results in an unsteady expansion stage, during which the interior of the pulsar wind bubble adjusts itself to the pressure of the surrounding SNR. Ultimately the expansion of the PWN becomes subsonic. However, when the pulsar is born with a kick velocity, the position of the active pulsar has become off-center with respect to the SNR at this stage, and will ultimately overtake and break through the shell of the decelerating SNR. Before this break-through event the motion of the pulsar becomes supersonic which deforms the head of the PWN, containing the active pulsar, into a bow shock (van der Swaluw et al. 1998).

In this paper I present results from a hydrodynamical simulation, carried out with the Versatile Advection Code (VAC)111See http://www.phys.uu.nl/~toth/. The simulation describes the evolution of a pulsar wind nebula when the pulsar is moving at a high velocity ($`\mathrm{1\hspace{0.33em}000}`$ km/sec) through its associated supernova remnant. The simulation shows that the position of the pulsar is strongly off-centered with respect to its pulsar wind bubble, after the interaction between the reverse shock and the PWN. Furthermore the simulation confirms the formation of a bow shock at a later stage of the PWN, as the motion of the pulsar becomes supersonic. The pulsar inside the SNR N157B is strongly off-centered with respect to its PWN. However the age of the remnant is rather young to explain the PWN morphology as a bow shock. Therefore I propose a scenario for which the pulsar wind nebula inside N157B has been recently crushed by the reverse shock. The simulation shows that this explains the off-centered position of the pulsar with respect to its PWN, the age of the remnant and the overall morphology of the PWN in a natural way.

## IMPORTANT TIMESCALES FOR PWN EVOLUTION

Initially the PWN is bounded by a strong shock propagating into the freely expanding ejecta of the progenitor star. As the forward shock of the SNR is expanding into the ISM, it starts to decelerate. McKee and Truelove (1995) define a timescale $`t_{\mathrm{ST}}`$, which marks the age of the remnant when it has swept up roughly 1.61 times the ejected mass $`M_{\mathrm{ej}}`$. They show that when the age of the remnant $`t_{\mathrm{snr}}`$, approximately equals $`t_{\mathrm{snr}}5t_{\mathrm{ST}}`$, the reverse shock hits the center of the SNR. Therefore one can roughly equal this timescale with the age at which the PWN interacts with the reverse shock:

$$t_{\mathrm{rev}}\mathrm{\hspace{0.33em}5}t_{\mathrm{ST}}=\mathrm{\hspace{0.33em}1\hspace{0.33em}045}E_{51}^{1/2}\left(\frac{M_{\mathrm{ej}}}{M_{}}\right)^{5/6}n_0^{1/3}\mathrm{years},$$
(1)
here $`E_{51}`$ is the total mechanical energy of the SNR in units of $`10^{51}`$ erg and $`n_0`$ is the ambient hydrogen number density assuming an interstellar composition of 10 H : 1 He. The above timescale is very close to the one given by Reynolds and Chevalier (1984) and marks the end of the stage of the PWN during which it is expanding supersonically.

After the interaction of the PWN with the reverse shock, the pulsar is propagating through a medium which has been reheated due to the passage of the reverse shock. Therefore the motion of the pulsar and the expansion of the PWN are both subsonic. Van der Swaluw et al. (1998) have shown that the motion of the pulsar becomes supersonic and deforms the pulsar wind bubble into a bow shock, when the position of the pulsar $`R_{\mathrm{psr}}`$, with respect to the radius of the SNR shock $`R_{\mathrm{snr}}`$, equals:

$$R_{\mathrm{psr}}0.677R_{\mathrm{snr}},$$
(2)
assuming a Sedov-Taylor expansion rate for the SNR. This occurs at roughly half the crossing time, i.e. the age of the SNR when the pulsar overtakes the shell of the remnant:

$$t_{\mathrm{cr}}\mathrm{\hspace{0.33em}1.4}\times 10^4E_{51}^{1/3}V_{1000}^{5/3}n_0^{1/3}\mathrm{years},$$
(3)
here $`V_{1000}`$ is the velocity of the pulsar in units of $`\mathrm{1\hspace{0.33em}000}`$ km/sec. Using the above timescales, one can roughly summarise the PWN evolution inside a Sedov-Taylor SNR as:

$``$ a supersonically expansion stage ($`t_{\mathrm{snr}}<t_{\mathrm{rev}}`$)

$``$ a subsonically expansion stage ($`t_{\mathrm{rev}}<t_{\mathrm{snr}}<0.5t_{\mathrm{cr}}`$)

$``$ a bow shock stage ($`0.5t_{\mathrm{cr}}<t_{\mathrm{snr}}<t_{\mathrm{cr}}`$)

## HYDRODYNAMICAL SIMULATIONS

### Simulation Method

I performed a hydrodynamical simulation of a SNR with a mechanical energy of $`E_0=10^{51}`$ erg and a total ejecta mass of $`M_{\mathrm{ej}}=3M_{}`$ expanding into a uniform ISM with a hydrogen number density of $`n_0=0.43`$. The simulation is performed in the restframe of the pulsar (wind), placed at the (explosion) center of the remnant with a constant luminosity $`L_{\mathrm{pw}}=10^{38}`$ ergs/sec and a pulsar velocity $`V_{\mathrm{psr}}=\mathrm{1\hspace{0.33em}000}`$ km/sec. The pulsar wind luminosity is taken to be constant throughout the simulation. This is required in order to resolve the pulsar wind termination shock. The results are certainly not changed qualitatively, because the total integrated energy input by the pulsar wind is only a $`5`$ % of the total mechanical energy of the SNR at the end of the simulation. This insures that the evolution of the pulsar wind bubble is completely determined by its surrounding evolving SNR. I used the Versatile Advection Code (Tóth 1996) to solve the equations of gas dynamics on a uniform grid with axial symmetry, using a cylindrical coordinate system $`(R,z)`$ in the restframe of the pulsar. I simulate the pulsar wind by depositing mass $`\dot{M}`$ and energy $`L_{\mathrm{pw}}`$ continuously in a few grid cells concentrated around the position of the pulsar. The terminal velocity of the pulsar wind is determined from these two parameters, i.e. $`v_{\mathrm{}}=\sqrt{2L_{\mathrm{pw}}/\dot{M}}`$, and has a value much larger then all the other velocities of interest in the simulation. The current version of the VAC code does not include relativistic hydrodynamics, therefore the best approach available is to take an adiabatic index of the relativistic fluid $`\gamma _{\mathrm{pwn}}=5/3`$.

Fig. 1. Logarithmic gray-scale representation of the density distribution at an age $`t_{\mathrm{snr}}=\mathrm{2\hspace{0.33em}250}`$ years.

### The Supersonically Expansion Stage (Figure 1)

Figure 1 shows the density profile of the PWN/SNR system at an age of $`t_{\mathrm{snr}}=\mathrm{2\hspace{0.33em}250}`$ years. The PWN has been dragged towards the reverse shock by the motion of the pulsar. The front part of the PWN has already been crushed by the reverse shock (the part indicated in the figure by interaction regime), it will take another $`500`$ years before the whole PWN will be completely crushed by the reverse shock. One can clearly distinguish the forward shock, contact discontinuity, PWN shock and the pulsar wind (termination) shock. In the figure the cross corresponds with the pulsar position.

Fig. 2. Logarithmic gray-scale representation of the density distribution at an age $`t_{\mathrm{snr}}=\mathrm{3\hspace{0.33em}000}`$ years.

### The subsonically Expansion Stage (Figure 2)

Figure 2 shows the density profile of the PWN/SNR system at an age of $`t_{\mathrm{snr}}=\mathrm{3\hspace{0.33em}000}`$ years, at which time the whole PWN has been crushed by the reverse shock. One can clearly distinguish the crushed relic part of the bubble blown in the initial stage of the PWN, when it was bounded by a strong shock. The pulsar is positioned at the head of the PWN (indicated by a cross again), but the feature ahead of the pulsar is not a bow shock but is the reflected reverse shock. Notice that although there is no bow shock around the PWN, the pulsar is strongly off-centered with respect to the PWN, after its interaction with the reverse shock.

### The bow shock stage (Figure 3)

Figure 3 shows the density profile of the PWN/SNR system at an age of $`t_{\mathrm{snr}}=\mathrm{15\hspace{0.33em}750}`$ years, at the end of the simulation. The pulsar is approaching the shell of the SNR and the head of the pulsar wind bubble, containing the active pulsar, has been deformed into a bow shock.The simulation confirms the formation of a bow shock at the head of the PWN, at half the crossing time $`t_{\mathrm{cr}}`$, when the position of the pulsar $`R_{\mathrm{psr}}`$ with respect to the shell of the SNR $`R_{\mathrm{snr}}`$ equals $`R_{\mathrm{psr}}0.677R_{\mathrm{snr}}`$, as predicted by van der Swaluw et al. (1998). Notice that the lengthscale of the PWN is small compared with the SNR shock structure. This confirms the approximation made by van der Swaluw et al. (2003) to model the break-through event of a pulsar wind through the shell of its SNR, by neglecting the curvature of the SNR blast wave.

Fig. 3. Logarithmic gray-scale representation of the density distribution at an age $`t_{\mathrm{snr}}=\mathrm{15\hspace{0.33em}750}`$ years.

### Summary

The results of the simulation performed in this paper clearly show the three stages of an evolving PWN in an expanding SNR; a supersonically expansion stage, a subsonically expansion stage and a bow shock stage. The simulation bears out that after the interaction between the reverse shock and the PWN, the pulsar is strongly off-centered with respect to its PWN, while the motion of the pulsar is still subsonic. Furthermore the simulation shows the formation of a bow shock at the head of the PWN, containing the active pulsar, at a later stage, when the motion of the pulsar becomes supersonic. However, the expansion of the relic part of the PWN, blown in the supersonically expansion stage, remains subsonic.

## THE PULSAR WIND NEBULA INSIDE SNR N157B

The bright X-ray object N157B is a young SNR dominated by plerionic emission originating from the PWN inside this object (Dickel and Wang 2003). The age of the remnant is approximately $`\mathrm{5\hspace{0.33em}000}`$ years old (Wang and Gotthelf 1998) and contains a 16 ms pulsar (Marshall et al. 1998). The velocity of this pulsar is high ($`V_{\mathrm{psr}}\mathrm{1\hspace{0.33em}000}`$ km/sec), if one assumes that the pulsar was born in the central region of the bright radio emission (Lazendic et al. 2000). Wang and Gotthelf (1998) argue for a bow shock interpretation of the PWN in N157B. A classical example of a pulsar wind bow shock inside its associated SNR is the remnant CTB80, with an estimated age of $`t_{\mathrm{snr}}\mathrm{100\hspace{0.33em}000}`$ years and an associated pulsar with a spindown luminosity of $`L_{\mathrm{pw}}4.0\times 10^{36}`$ erg/sec. These values show a remarkable contrast with the parameters of the SNR N157B ($`t_{\mathrm{snr}}\mathrm{5\hspace{0.33em}000}`$ years and $`L_{\mathrm{pw}}4.8\times 10^{38}`$ erg/sec). Furthermore from Figure 2 of Wang et al. (2001) it seems like the position of the pulsar is more or less centered in the SNR. This is in contrast with the analysis performed by van der Swaluw et al. (1998), which predicts $`R_{\mathrm{psr}}/R_{\mathrm{snr}}0.677`$ for a bow shock PWN, a result confirmed by the simulation in this paper.

In order to solve the issues mentioned above, I propose the following evolutionary stage for the PWN inside SNR N157B: the PWN has recently been crushed by the reverse shock of its associated SNR. Due to the combination of the reverse shock interaction and the high velocity of the pulsar, the position of the pulsar is strongly off-centered with respect to its PWN, as shown by the results of the simulations performed in this paper. Figure 2 represents the current morphology of N157B, where the relic PWN represents the central bright parts of the radio and X-ray emission in N157B. This scenario automatically implies the absence of a bow shock in N157B. Furthermore the age of N157B is in rough agreement with this scenario ($`t_{\mathrm{rev}}<t_{\mathrm{N157B}}<0.5t_{\mathrm{cr}}`$).

## CONCLUSION

I have performed a hydrodynamical simulation of a PWN, when the pulsar is moving at a high velocity through its associated SNR. The simulation shows that one can roughly distinguish between a supersonically expanding PWN, a subsonically expanding PWN and a stage during which part of the PWN is bounded by a bow shock, due to the supersonic motion of the pulsar. After the interaction between the PWN and the reverse shock the position of the pulsar has become off-centered with respect to its PWN. Therefore one should be careful in distinguishing a young subsonically expanding PWN ($`t_{\mathrm{rev}}<t_{\mathrm{pwn}}<0.5t_{\mathrm{cr}}`$) and an older PWN with a bow shock around its pulsar ($`0.5t_{\mathrm{cr}}<t_{\mathrm{pwn}}<t_{\mathrm{cr}}`$). In the case of N157B, I have argued for a no as an answer to the title of this paper.

E-mail address of E. van der Swaluw swaluw@rijnh.nl


format: text

text:
# Compact 𝑈(1) gauge theories in 2+1 dimensions and the physics of low dimensional insulating materials

## 1 Introduction

Effective abelian gauge theories of doped Mott insulators are very popular in the condensed matter literature. Despite the enormous amount of papers published over the last ten years, there are still many points of disagreement among researchers, mostly related to poorly understood issues Kim ; Nayak ; Nagaosa ; Matsui . One key point is the interplay between the confinement/deconfinement properties of the many theories available, and spin-charge separation. It is usually believed that if deconfinement occurs, then spin-charge separation would also occur. However, it remains a controversial issue if deconfinement is possible for gauge fields coupled to matter with the fundamental charge. The problem here is that most $`U(1)`$ effective gauge theories are studied in $`d=2+1`$ space-time dimensions. For this dimensionality $`U(1)`$ gauge fields are strongly confining. For instance, if we neglect the coupling to matter fields we obtain that pure compact Maxwell theory permanently confines electric charges in $`d=2+1`$ dimensions Polyakov . This is in contrast with pure compact Maxwell theory in $`d=3+1`$ dimensions where a deconfinement transition occurs Peskin . A natural question to ask in this context is whether the coupling to matter fields is able to change the permanent confinement in $`2+1`$ dimensions. One example where the answer to this question is affirmative corresponds to the case where bosons are coupled to an abelian gauge field with integer charge $`q>1`$ FradShe . This is most easily understood by considering the lattice gauge theory version of such a model, whose action is given by

$$S=\beta \underset{i,\mu }{}[\mathrm{cos}(_\mu \theta _iqA_{i\mu })1]\kappa \underset{i,\mu ,\nu }{}[\mathrm{cos}(F_{i\mu \nu })1],$$
(1)
where we have $`_\mu \theta _i\theta _{i+\widehat{\mu }}\theta _i`$ and $`F_{i\mu \nu }_\mu A_{i\nu }_\nu A_{i\mu }`$. The above action corresponds to the compact lattice abelian Higgs (CLAH) model. When $`\beta \mathrm{}`$ the theory becomes a $`Z_q`$ gauge theory which is known to have a deconfining transition in $`2+1`$ dimensions. For $`\kappa \mathrm{}`$ we obtain the $`XY`$ model, which also undergoes a phase transition. The critical points corresponding to these two limiting cases are connected by a critical line separating the confining from the deconfining phase. There is no Coulomb phase in $`2+1`$ dimensions. The phase diagram and the critical properties of the action (1) were studied in detail using large scale Monte Carlo simulations in Ref. Sudbo . The CLAH model appears in many contexts in condensed matter physics. We shall cite only two important recent examples. In the $`q=2`$ version it arises as an effective theory for the two-dimensional quantum Heisenberg antiferromagnet (QHA) Sachdev . There the gauge field originates from the Berry phase. Another recent example corresponds to the strongly correlated limit of a bosonic model considered by Senthil and Motrunich SM . In this context the CLAH model describes a transition from an ordinary Mott insulating phase to a fractionalized insulating phase SM .

However, bosonic actions like the one in Eq. (1) are only part of the complete effective action associated with a doped Mott insulator. It usually contains also fermions coupled to a gauge field. IL ; LN  The traditional route starts with a slave boson representation of the $`tJ`$ model where the projected electron operator is written as a composite particle, $`c_{i\alpha }=b_i^{}f_{i\alpha }`$, where $`b_i`$ is an auxiliary boson and $`f_{i\alpha }`$ is an auxiliary fermion. The auxiliary fields obey the constraint $`b_i^{}b_i+_\alpha f_{i\alpha }^{}f_{i\alpha }=1`$ at each lattice site. In the resulting effective action the auxiliary fields interact through a compact gauge field. The most popular way of doing things attribute the charge of electron to the bosonic field. Thus, the fermion would have no charge and would only carry the spin degree of freedom. If the effective model undergoes a deconfinement transition fermions and bosons will have an independent dynamics and in this way we can say the spin and charge are separated. Nayak Nayak  pointed out recently that things are not so simple because the assignment of the whole electron charge to the auxiliary boson is completely arbitrary. We could well say that the boson carries charge $`\delta e`$ while the fermion carries charge $`(1\delta )e`$. The constraints of the theory, which are enforced by the coupling to the gauge field, would then ensure that the physics is unaffected by this arbitrary choice. This scenario implies that the auxiliary fields should not be associated with the physical charge and spin excitations of the theory. The auxiliary fields introduced this way are not part of the physical spectrum. This situation is reminiscent of the analysis made by Mudry and Fradkin some time ago. Mudry

In this paper, we review some recents results on the model (1). We start discussing in the Section 2 the finite-size scaling (FSS) analysis of the model in the $`q>1`$ case. We employ a new and very accurate method for extracting the critical exponents. In Section 3 we discuss the phase diagram, which is obtained for $`q=2,3,4,5`$. In Section 4 we discuss the $`q=1`$ case, whose phase diagram is presently a matter of intensive debate KNS ; Herbut ; Chernodub .

## 2 The universality class of the deconfinement transition for $`q>1`$

In 1981 Bhanot and Freedman Bhanot  made a finite-size scaling (FSS) analysis of the CLAH model for $`q=2`$ and $`d=2+1`$. They obtained the phase diagram of the model but the critical exponents were not calculated. In order to determine the universality class of the model it is necessary to compute the critical exponents. The accurate determination of the critical behavior is not an easy task. In principle, on the basis of the results of Fradkin and Shenker FradShe , we might think that the universality class will be the same as of the $`Z_q`$ spin model, except for the limit $`\kappa \mathrm{}`$ where the universality class is obviously that of the $`XY`$ model. If we assume that this is indeed so, we obtain that the critical exponents for the $`q=2`$ case have Ising values. However, we have recently shown through large scale Monte Carlo simulations that the situation is more complicated Sudbo . For example, for $`q=2`$ we have obtained that the exponents are those of the Ising model in a large part of the phase diagram in the $`\kappa \beta `$-plane, but there is a finite interval $`\kappa _1<\kappa <\kappa _2`$ where the exponents vary continuously, before reaching $`XY`$ values for $`\kappa >\kappa _2`$. Therefore, it seems that the model features a fixed line rather than a fixed point and belongs to a new universality class.

The FSS analysis of Bhanot and Freedman Bhanot  relies on the second moment of the free energy, i.e., the specific heat. In this case it is well known that very large system sizes are needed in order to identify a genuine critical behavior to high accuracy. We have shown that much better results are obtained by performing a FSS analysis based on the third moment of the free energy. Such a procedure is not only more accurate: it allows also for an independent determination of the exponents $`\nu `$ and $`\alpha `$, providing in this way a check of hyperscaling Sudbo . The reason for this lies in the double-peak structure of the third moment. The FSS in this case is such that the width between the negative and positive peak scales as $`L^{1/\nu }`$, while the height of the positive peak with respect to the negative one scales as $`L^{(1+\alpha )/\nu }`$. This scaling behavior is shown schematically for a generic third moment $`M_3`$ in Fig. 1.

The critical exponents as a function of $`\kappa /\beta `$ are shown in Fig. 2 for $`q=2`$. It is clearly seen that the Ising critical behavior connects to the $`XY`$ behavior through a regime of continuously varying critical exponents.

## 3 The phase diagram for $`q>1`$

We have obtained the phase diagram for $`q=2,3,4,5`$. Among these values of $`q`$, only the $`q=3`$ case exhibits a first-order phase transition up to some point where it changes to second-order. This behavior is consistent with the fact that the three-state Potts model exhibits a first-order transition in three dimensions. This model is dual to the $`Z_3`$ gauge theory. The phase diagram for $`q=3`$ case therefore features a tricritical point. The tricritical point in the phase diagram is estimated as $`(\beta _{\mathrm{tri}},\kappa _{\mathrm{tri}})=(1.23\pm 0.03,1.73\pm 0.03)`$, corresponding to a ratio $`\kappa _{\mathrm{tri}}/\beta _{\mathrm{tri}}=1.39\pm 0.06`$. At this point it is worth making the following remark. The non-compact version of the model (1) is known to exhibit a second-order phase transition for all values of $`q`$. As matter of fact, the value of $`q`$ is not important in that case, and can be absorbed into a redefinition of the gauge field. The universality class corresponds to the so called “inverted” $`XY`$ transition Dasgupta . However, if the radial part of the scalar field is allowed to fluctuate, thus generalizing the non-compact version of (1), then it is possible to show using duality arguments that the resulting model has a tricritical point approximately at $`\kappa _{\mathrm{tri}}/\beta _{\mathrm{tri}}0.64`$ Kleinert . This estimate was confirmed recently by large scale Monte Carlo simulations Mo . In the case of the $`q=3`$ CLAH model there is no need to consider the fluctuations of the radial component of the field to obtain a tricritical point. The origin of this tricritical point is completely different from the one of the non-compact model. Anyway, it is interesting to note that the ratio $`\kappa _{\mathrm{tri}}/\beta _{\mathrm{tri}}`$ in the $`q=3`$ CLAH model is more than twice larger than the corresponding ratio in the non-compact case.

In the context of the model discussed recently by Senthil and Motrunich SM , the phase diagram for the $`q=2`$ case has the following physical meaning. The confined phase corresponds to gap to integer charge excitations and is therefore interpreted as an ordinary Mott insulating phase. The deconfined phase, on the other hand, is the same as the Higgs phase. This gives a gap to half-integer charge excitations and is interpreted as a fractionalized insulator. The model can in principle be artificially realized with present day technology by building an array of Josephson junctions of a particular type SM .

## 4 The $`q=1`$ case

The $`q=1`$ case was generally believed to not exhibit any phase transition at all. The coupling to matter fields would not change the scenario of the $`d=2+1`$ compact Maxwell theory, where it is known that no phase transition takes place Polyakov : the electric charges would never deconfine. A research initiated with the study of the QHA showed that things seem not to be so simple IL . The main point is that by integrating out the matter fields we obtain to quadratic order in the gauge field fluctuations the effective action

$$S_{\mathrm{eff}}d^3xF_{\mu \nu }\frac{1}{\sqrt{^2}}F_{\mu \nu }.$$
(2)
Therefore, while in the absence of matter fields the monopoles interact through a $`1/r`$ potential, in the presence of matter fields this behavior changes to $`\mathrm{ln}r`$. In the former case it is well known that upon dualizing the theory a sine-Gordon action is obtained. In the latter case, however, an anomalous sine-Gordon action arises KNS :

$$S=\frac{1}{8\pi ^2K}d^3x[\phi (^2)^{3/2}\phi 2z\mathrm{cos}\phi ],$$
(3)
where $`K=1/g^2`$, with $`g`$ being the gauge coupling. Thus, we arrive at a behavior similar to the two-dimensional case. By simply looking at the scaling behavior of the fugacity of the monopole gas, one would be lead to conclude that a Kosterlitz-Thouless (KT) phase transition happens in three space-time dimensions. Note that this KT transition is at zero temperature and has nothing to do with other KT transitions obtained by dimensional reduction due to temperature effects. However, it is not enough to study the scaling of just the fugacity. Since we are in $`d=2+1`$ a more thorough analysis is needed in order to fully establish the KT behavior. To this end, it is necessary to also know the scaling behavior of the “stiffness” $`K`$, since the scaling of the fugacity and the scaling of the stiffness mutually influence each other in a subtle way.

For the usual sine-Gordon model in $`d`$ dimensions we obtain the following coupled recursion relations for the fugacity and the stiffness Kosterlitz

$$\frac{dK^1}{dl}=4\pi ^2y^2(2d)K^1,$$
(4)
$$\frac{dy}{dl}=\left[d2\pi ^2f(d)K\right]y,$$
(5)
where $`f(d)=(d2)\mathrm{\Gamma }[(d2)/2]/(4\pi )^{d/2}`$. For $`d=2`$ the above equations imply the existence of a line of fixed points characteristic of the KT transition. For $`d=3`$, however, no fixed point is found and therefore no phase transition takes place. This is consistent with Polyakov’s result for the compact Maxwell theory in $`d=2+1`$ Polyakov . The issue here is how the above recursion relations are modified in the case of the anomalous sine-Gordon model, Eq. (3). To investigate this in great generality we will consider screening in a Coulomb gas of monopoles with a propagator of the form $`1/|p|^\sigma `$ in $`2d<4`$ dimensions with $`\sigma >0`$. The particular case $`d=\sigma =3`$ corresponds to our anomalous sine-Gordon theory. The bare potential is given by $`U_0(r)=4\pi ^2KV(r)`$, where

$$V(r)=\frac{\mathrm{\Gamma }\left(\frac{d\sigma }{2}\right)}{2^\sigma \pi ^{d/2}\mathrm{\Gamma }(\sigma /2)}[(\mathrm{\Lambda }r)^{\sigma d}1],$$
(6)
with $`\mathrm{\Lambda }`$ being an ultraviolet cutoff. Note that for $`d=\sigma `$ we obtain a potential $`\mathrm{ln}r`$. The bare electric field is given by $`E_0(r)=4\pi ^2KA(d,\sigma )r^{\sigma d1}/r_0^{\sigma d}`$, where $`r_01/\mathrm{\Lambda }`$ and $`A(d,\sigma )=(d\sigma )\mathrm{\Gamma }[(d\sigma )/2]/[2^\sigma \pi ^{d/2}\mathrm{\Gamma }(\sigma /2)]`$. The bare electric field is renormalized by the other dipoles which are treated as a dielectric medium. The renormalized electric field is then given by

$$E(r)=\frac{4\pi ^2KA(d,\sigma )r^{\sigma d1}}{\epsilon (r)},$$
(7)
where $`\epsilon (r)`$ is the scale-dependent dielectric constant of the medium. This problem can be solved self-consistently for a renormalized potential $`U(r)`$ derived from the above electric field KNS . In this way we obtain the recursion relations

$$\frac{dK^1}{dl}=4\pi ^2y^2(\sigma d)K^1,$$
(8)
$$\frac{dy}{dl}=\left[d\eta _y2\pi ^2A(d,\sigma )K\right]y,$$
(9)
where the anomalous dimension of the fugacity is given by $`\eta _y=(\sigma 2)/2`$. The above equations reduce to Eqs. (4) and (5) when $`\sigma =2`$. The case relevant to the $`q=1`$ CLAH model is $`\sigma =d=3`$. In this case Eqs. (8) and (9) are very similar to the usual KT recursion relations, except for the presence of the anomalous scaling dimension of the fugacity, $`\eta _y`$ which is nonzero in our case and given by $`\eta _y=1/2`$. Thus, we obtain that the monopole gas undergoes a KT-like phase transition. This is in contrast with Polyakov’s compact Maxwell theory where the monopoles are always in the plasma phase. Using the usual duality arguments, we identify the plasma phase with the confined phase for the electric charges. The dielectric phase of the monopole gas is accordingly identified with the deconfined phase for the electric charges.

At this point an important remark is in order. The above screening analysis was made in real space, in the spirit of the original Kosterlitz and Thouless paper KT . Indeed, we have considered a local scale dependent dielectric constant. Screening arguments in momentum space usually involve a dielectric constant which is local in momentum space, leading to an electric displacement vector of the form $`𝐃(𝐫)=d^dr^{}ϵ(𝐫𝐫^{})𝐄(𝐫^{})`$, while in our case we have simply $`𝐃(𝐫)=ϵ(𝐫)𝐄(𝐫)`$. In the present context this may be understood by using a classical electrostatic argument. Let us consider the potential (6) in the case $`\sigma =d`$, such that we have a $`\mathrm{ln}r`$-potential in $`d`$ dimensions. Due to Gauss’s theorem in $`d`$-dimensions, the field equation is given by Dixit  $`(𝐄r^{2d})=S_d\rho (r)`$, where $`S_d`$ is the solid angle in $`d`$-dimensions. This can be cast in a more familiar form $`𝐃=S_dr_0^{d2}\rho (r)`$ by introducing the “dielectric constant” $`\epsilon _0(r)=(r/r_0)^{2d}`$. Note that for $`d>2`$ this dielectric constant vanishes at large distances. This is rather unusual since in classical electromagnetism it can be shown that the dielectric constant is always greater than one. However, such a situation corresponds precisely to the anti-screening mechanism discussed many years ago in theories of confinement Kogut . Thus, already at this level, even before taking the dipoles into account, we can define a kind of scale-dependent dielectric constant which is local in real space.

## 5 Conclusions

In this paper we discussed the deconfinement transition in the CLAH model, Eq. (1), in $`d=2+1`$. To each value of the charge $`q`$ corresponds a distinct universality class, in contrast to the non-compact model. While a deconfined transition is expected for $`q>1`$, the case $`q=1`$ is still controversial KNS ; Herbut ; Chernodub  and needs further study. One feature that also deserves further study is the possible existence of a fixed line in the $`q=2`$ case. The numerical results strongly suggest this possibility but multicritical fixed points are also possible.

###### Acknowledgements.

FSN would like to thank Boris Pioline and Christoph Berger for the invitation to give the talk in the EPS HEP2003 conference, upon which this paper is based. We thank the Deutsche Forschungsgemeinschaft (DFG) and the Norwegian Research Council, Project No. 131520/432, for financial support.


text:
# General Relativistic Magnetohydrodynamic Simulations of Collapsars

## 1 Introduction

Gamma ray bursts (GRBs) are one of the most enigmatic and most energetic events in the universe (e.g., Klebesadel, Strong, & Olson 1973; Fishman & Meegan 1995; van Paradijs, Kouveliotou, & Wijers 2000). GRBs and the afterglows are well described by the fireball model (e.g., Piran 1999; Mészáros 2002), in which a relativistic outflow is generated from a compact central engine. Rapid temporal decay of several afterglows is consistent with the evolution of a highly relativistic jet with bulk Lorentz factors $`10^210^3`$ (e.g., Sari, Piran, & Halphen 1999). The formation of relativistic jets from a compact central engine remains one of the major unsolved problems in GRB models.

What is the central engine of GRBs? From recent observations, some evidence was found for a connection between GRBs and the death of massive stars. Analyses of host galaxies show a correlation between star-forming regions and the position of GRBs inside the host galaxy (Bloom, Kulkarni, & Djorgovski 2002). A “bump” resembling that seen in light curves of Type Ic supernovae has also been detected in the optical afterglow of several GRBs (e.g., Bloom et al. 2002; Garnavich et al. 2003). GRB 980425 has been associated with an optical supernova, SN 1998bw (e.g., Iwamoto et al. 1998; Woosley, Eastman, & Schmidt 1999). From the early-phase observation of a very bright afterglow (e.g., Uemura et al 2003; Price et al. 2003) it has been revealed that GRB 030329 is likely to be in association with SN 2003dh (Stanek et al. 2003; Hjorth et al. 2003). Several authors (Frail et al. 2001; Panaitescu & Kumar 2001; Bloom, Frail, & Kulkarni 2003) have studied beaming angles and energies of a number of GRBs. They have found that central engines of GRBs release supernova-like energies ($`10^{51}`$ erg). It is thus probable that a major subclass of GRBs is a consequence of the collapse of a massive star.

One of the most attractive scenarios involving massive stars is the collapsar model (Woosley 1993; MacFadyen & Woosley 1999). A “collapsar” is a rotating massive star. The collapsar model is divided into two classes according to the formation history of the black hole. A “type I collapsar” is a failed supernova. The collapse of an iron core leads temporarily to a neutron star formation. A black hole, however, is formed quickly within a few seconds as a result of the accretion of matter through the stalled shock wave. In the mean time, infalling envelope matter is slowed by rotation in the equatorial plane and forms a disk. The formation and propagation of relativistic flows from a type I collapsar have been studied numerically in both Newtonian (MacFadyen & Woosley 1999) and special relativistic hydrodynamic simulations (Aloy et al. 2000; Zhang, Woosley, & MacFadyen 2003). The other model is a “type II collapsar,” in which a black hole is formed over a longer period of time accompanying a successful supernova. The supernova generates an outgoing shock and ejects all the helium and heavy elements outside the neutron star. However, some of the post-shock gas fails to reach escape velocity and is pulled back toward the proto-neutron star by gravity. When enough gas has fallen back, the neutron star collapses to a black hole. If the infalling matter has sufficient angular momentum, it forms a disk at the same time. The formation and propagation of relativistic flows from a type II collapsar have been studied numerically in Newtonian hydrodynamic simulations (MacFadyen, Woosley, & Heger 2001). However, these previous numerical simulations of collapsar models do not fully address the outflow mechanism. In these simulations, in fact, the authors estimate the energy of a jet assuming neutrino annihilation or magnetohydrodynamic (MHD) process as a source and input a jet with that energy from the inner boundary. There have been a lot of studies of astrophysical jets with Newtonian (e.g., Uchida & Shibata 1985; Shibata & Uchida 1986, 1990; Kudoh, Matsumoto, & Shibata 1998; Kato, Kudoh, & Shibata 2002) and relativistic MHD simulations (e.g., Koide, Shibata, & Kudoh 1998, 1999). They have fully addressed the formation, acceleration, and collimation of jets from accretion disks. We think that the formation of relativistic flow from collapsars should also be studied by MHD simulations. In fact, several authors (Cameron 2001; Wheeler, Meier, & Wilson 2002) proposed that the relativistic jets are generated by the MHD process when the massive star collapses to a rapidly rotating neutron star.

It is suspected that large-scale magnetic fields play an important role in the formation of a GRB. Magnetic fields are suitable for extracting energy on the burst time-scale from the debris disk in the collapsar scenario. If the field is initially weak, it is likely to be exponentially amplified as a result of differential rotation and/or dynamo. Large-scale magnetic fields can help guide and collimate the outflow and also contribute to its acceleration (e.g., Usov 1994; Thompson 1994; Meszaros & Rees 1997; Katz 1997; Kluźniak & Ruderman 1998; Lyutikov & Blackman 2001; Drenkhahn & Spruit 2002; Vlahakis & Königl 2001, 2003ab). Even if the flow is not magnetically driven, the field should be strong enough to account for the observed synchrotron emission (Spruit, Daigne & Drenkhahn 2001). If the outflow is largely Poynting flux dominated, it can solve the baryon contamination problem. Magnetic fields have been favored in this respect for driving GRB outflows.

The effect of stellar rotation and intrinsic magnetic fields on gravitational collapse of massive stars was numerically studied by several authors (LeBlanc & Wilson 1970; Symbalisty 1984; Ardeljan et al. 2000). Symbalisty (1984) computed the magnetorotational core collapse of a $`15M_{}`$ star by numerically solving the 2.5-dimensional MHD equations together with the neutrino transport. The simulations showed the formation of two oppositely directed, high-density, supersonic jets in the combination of a rapid rotation and a very strong dipole magnetic field. Khokhlov et al. (2001) assumed that the jets are generated by a magnetorotational mechanism when a stellar core collapses into a neutron star and simulated the process of the jet propagation through the star. Proga et al. (2003) simulated the collapsar model by using a pseudo-Newtonian MHD code including some essential microphysics (physical equation of state \[EOS\], photodisintegration of nuclei, and neutrino cooling) for GRBs.

In this study we perform 2.5-dimensional general relativistic MHD(GRMHD) simulations of the gravitational collapse of a rotating star with magnetic field as a model for a collapsar. The collapsar is in some sense an anisotropic supernova, and it is considered that relativistic jets from collapsars are launched by MHD processes in accreting matter and/or by neutrino annihilation. We investigate the physics of the formation of jets, the acceleration force on jets, and the dependence of the acceleration of jets on the initial magnetic field strengths and on the initial rotational velocity. We describe the numerical method in our simulations in section 2 and present our results in section 3. The summary and discussion are given in section 4.

## 2 Numerical Method

### 2.1 Basic Equations

In order to study the formation of relativistic jets from a collapsar, we use a 2.5-dimensional GRMHD code with spherical polar coordinates ($`r,\theta ,\varphi `$) (i.e., $`/\varphi =0`$, but $`v_\varphi `$ and $`B_\varphi `$ are nonzero; Koide 2003. The method is based on the general relativistic formulation of the conservation laws of particle-number and energy-momentum, Maxwell equations, and Ohm’s law with no electrical resistance (ideal MHD condition) on a curved space-time (Weinberg 1972; Koide et al. 1998, 1999; Koide 2003).

The space-time $`(x^0,x^1,x^2,x^3)=(ct,x^1,x^2,x^3)`$ is described by the metric $`g_{\mu \nu }`$, where the line element $`ds`$ is given by $`(ds)^2=g_{\mu \nu }dx^\mu dx^\nu `$. Here $`c`$ is the speed of light. The basic equations of GRMHD in four-dimensional space-time are

$$_\nu (\rho U^\nu )=\frac{1}{\sqrt{g}}\frac{}{x^\nu }(\sqrt{g}\rho U^\nu )=0,$$
(1)
$$_\nu T^{\nu \mu }=\frac{1}{\sqrt{g}}\frac{}{x^\nu }(\sqrt{g}T^{\mu \nu })+\mathrm{\Gamma }_{\sigma \nu }^\mu T^{\sigma \nu }=0,$$
(2)
$$_\mu F_{\nu \lambda }+_\nu F_{\lambda \mu }+_\lambda F_{\mu \nu }=0,$$
(3)
$$_\mu F^{\mu \nu }=\mu _0J^\nu ,$$
(4)
where $`_\nu `$ is the covariant derivative, $`g`$ is the determinant of $`g_{\mu \nu }`$ (as a matrix), and $`\mathrm{\Gamma }_{\mu \nu }^\lambda \frac{1}{2}g^{\lambda \sigma }\left(\frac{g_{\mu \nu }}{x^\sigma }+\frac{g_{\nu \sigma }}{x^\mu }+\frac{g_{\sigma \mu }}{x^\nu }\right)`$ are the Christoffel symbols. $`U^\nu `$ and $`J^\nu =(c\rho _e,J^1,J^2,J^3)`$ are four-velocity and four-current density, respectively ($`\rho _e`$ is the electric charge density). The general relativistic energy momentum tensor $`T^{\mu \nu }`$ is given as

$$T^{\mu \nu }=pg^{\mu \nu }+(e_{int}+p)U^\mu U^\nu +F_\sigma ^\mu F^{\nu \sigma }\frac{1}{4}g^{\mu \nu }F^{\lambda \kappa }F_{\lambda \kappa }$$
(5)
where $`F^{\mu \nu }`$ is the electromagnetic field-strength tensor, $`F_{\mu \nu }=_\mu A_\nu _\nu A_\mu `$ and $`A^\mu =(\mathrm{\Phi }_e/c,A^1,A^2,A^3)`$ is the four-vector potential ($`\mathrm{\Phi }_e`$ is the electro-static potential). The electric field $`E_i`$ and the magnetic field $`B_i`$ are given by $`E_i=cF_{i0}`$ ($`i=1,2,3`$) and $`B_1=F_{23}`$, $`B_2=F_{31}`$, $`B_3=F_{12}`$, respectively. The electric field $`E_i`$ and the magnetic field $`B_i`$ are normalized as $`E_i=E_i^{}/\sqrt{\mu }`$ and $`B_i=B_i^{}/\sqrt{\mu }`$, respectively, where $`\mu `$ is the magnetic permeability and quantities with asterisks are given in MKSA units (SI units). The scalars $`\rho `$, $`p`$, and $`e_{int}`$ are proper mass density, proper pressure, and proper internal energy density, respectively. We assume mainly for numerical simplicity that matter can be described as an ideal gas; hence, $`e_{int}=\rho c^2+p/(\mathrm{\Gamma }1)`$, where $`\mathrm{\Gamma }`$ is the specific heat ratio although we know that the gamma law EOS is not very appropriate for the gravitational collapse of massive stars. We further assume that $`\mathrm{\Gamma }=5/3`$ in the simulations for numerical reasons although it is more appropriate to adopt $`\mathrm{\Gamma }=4/3`$ for the gas of current interest. In addition to the equations, we assume the infinite electric conductivity condition

$$F_{\mu \nu }U^\nu =0.$$
(6)
Then, equations (1)-(3) together with equation (6) are a closed system of equations. Equation (4) is used only to calculate the current density $`J^\mu `$.

If we assume that the off-diagonal elements of the metric $`g_{\mu \nu }`$ vanish,

$$g_{\mu \nu }=0(\mu \nu ),$$
(7)
and we use the notation

$$g_{00}=h_0^2,g_{11}=h_1^2,g_{22}=h_2^2,g_{33}=h_3^2,$$
(8)
then the line element can be written as

$$(ds)^2=g_{\mu \nu }dx^\mu dx^\nu =h_0^2(cdt)^2+\underset{i=1}{\overset{3}{}}h_i^2(dx^i)^2.$$
(9)
We employ the 3+1 formalism (Thorne, Price, & Macdonald 1986; Koide et al. 1999; Koide 2003). The velocity of light is written explicitly. In the non-relativistic limit ($`c\mathrm{}`$), they reduce to the standard Newtonian MHD equations.

In the GRMHD code, a simplified total variation diminishing (TVD) method is employed (Davis 1984). This method is useful because it requires only the maximum speed of physical waves but not the eigenvalues of the Jacobian of the equations.

De Villiers & Hawley (2003) and Gammie, McKinney, & Tóth (2003) have developed a new GRMHD code independently. De Villiers & Hawley (2003) used a non-conservative scheme that is based on techniques developed by Hawley, Smarr, & Wilson (1984) for axisymmetric hydrodynamics around black holes. On the other hand, Gammie, McKinney, & Tóth (2003) used a conservative, shock-capturing scheme. They have performed some one- and two- dimensional tests and have obtained consistent results with each other.

We do not consider the evolution of the metric because the accreted mass is sufficiently small in the time-scale of the simulations. The current calculation focuses on GRMHD with simplified microphysics. Although photodisintegration of bound nuclei and cooling by neutrino emission (thermal and especially Urca processes) have been shown to be of critical importance for the dynamics of collapsar accretion disks (Popham, Woosley & Fryer 1998; MacFadyen & Woosley 1999), we ignore these effects in the current calculations for simplicity. It is hoped that future GRMHD simulations will include essential microphysics (neutrino cooling, physical EOS, photodisintegration) to realistically simulate collapsar accretion disks.

### 2.2 Initial Condition

As for the initial model, we have the collapsar model in mind. In principle, we should start calculations from a realistic progenitor model with rotation and magnetic field. Although some progress has been made lately on the precollapse evolution of massive stars with rotation and magnetic field (Heger et al. 2003), their models have a lot of uncertainty yet. Hence, we take the following pragmatic approach: Core collapse is assumed to fail in this star, and a few $`M_{}`$ black hole is inserted by hand into the calculation. The simulations presented in the paper follow the accretion of gas into a black hole assumed to have formed before the calculation begins. Furthermore, We use Bruenn’s realistic one-dimensional supernova model(Bruenn 1992) as a rough guide to what the density structure exterior to the black hole might be. The gas in Bruenn’s model must have collapsed before the present calculation began to form the black hole that is assumed to be present. However, we use the scaling in that model as a guide for the scale-free stellar structure. We employ only the profiles of the density, pressure, and radial velocity as our initial condition. In this way, we can discuss generic features of the dynamics.

We consider a non-rotating black hole. The metric is given by

$$g_{00}=1\frac{r_S}{R},g_{11}=\left(1\frac{r_S}{R}\right)^1,g_{22}=R^2,g_{33}=R^2\mathrm{sin}\theta ,$$
(10)
where $`r_S2GM_{}/c^2`$ is the Schwarzschild radius. We add the effect of stellar rotation and intrinsic magnetic field in the simulations as the initial condition. The initial rotational velocity distribution is assumed to be a function of the distance from the rotation axis,

$$r=R\mathrm{sin}\theta ,$$
(11)
only:

$$v_\varphi =v_0\frac{x_0^2}{r^2+x_0^2}r.$$
(12)
Here, $`v_0`$ is a model parameter for rotational velocity. We fix $`x_0=100r_S`$ in this paper. The initial magnetic field is assumed to be uniform and parallel to the rotational axis. This is known as the Wald solution for a non-rotating black hole: $`B_R=B_0\mathrm{cos}\theta ,B_\theta =\alpha B_0\mathrm{sin}\theta `$. $`B_0`$ is the overall magnetic field strength and $`\alpha =g_{00}`$ is the lapse function. Figure 1 shows the schematic picture of our simulation. The distributions of various physical quantities (density $`\rho `$, pressure $`p`$, each velocity component $`v_r,v_\varphi ,v_z`$, each magnetic field component $`B_r,B_\varphi ,B_z`$, and plasma beta $`\beta =p_{gas}/p_{mag}`$) on the equatorial plane for the initial state are shown in Figure 2. In the distribution of density and pressure, we can see that the weak shock is sitting at about 18 $`r_S`$. We input the stellar rotation up to 18 $`r_S`$ because the accreted mass for stellar matter in the pre-shock region is sufficiently small in the time scale. We neglect rotation in the pre-shock region for numerical simplicity. The rotation is almost uniform (see Eq. ). This initial rotation profile is similar to the previous simulations of rotational core collapse (Mönchmeyer & Müller 1989; Yamada & Sato 1994).

We want to emphasize that our simulations are scale free. The normalization units are summarized in Table 1. We use typical values for normalizations. If we set a black hole mass $`M_{BH}3M_{}`$, the Schwarzschild radius $`r_S`$ is about $`8\times 10^5`$ cm and the time unit $`\tau _S(r_S/c)`$ is $`3.0\times 10^5`$ s. The units of magnetic field strength and pressure depend on the normalization of the density. If we take, for example, the density unit $`\rho _0=10^{10}\mathrm{g}/\mathrm{cm}^3`$, the magnetic field strength unit is $`3.0\times 10^{14}`$ G and the pressure unit is $`P=10^{31}\mathrm{dyn}/\mathrm{cm}^2`$. The models computed in this paper are summarized in Table 2. These units are the values at the last stable Keplerian orbit of a non-rotating black hole on the equatorial plane \[$`(r,z)=(r_0,0)`$\]. The variables with subscript “0” represent the values at the last stable orbit throughout this paper.

We also use the following nondimensional parameters:

$$E_{th}=\frac{V_{s0}^2}{V_{K0}^2},$$
(13)
$$E_{mag}=\frac{V_{A0}^2}{V_{K0}^2},$$
(14)
$$E_{rot}=\frac{V_\varphi ^2}{V_{K0}^2},$$
(15)
where $`V_{s0}`$ is the relativistic sound speed,

$$V_{s0}=c\sqrt{\frac{\mathrm{\Gamma }p_0}{\rho _0c^2+\mathrm{\Gamma }p_0/(\mathrm{\Gamma }1)}},$$
(16)
$`V_{A0}`$ is the Alfvén velocity,

$$V_{A0}=c\frac{B_0}{\sqrt{\rho _0c^2+\mathrm{\Gamma }p_0/(\mathrm{\Gamma }1)+B_0^2}},$$
(17)
and $`V_{K0}`$ is Keplerian velocity,

$$V_{K0}=c/[2(r_0/r_S1)]^{1/2}.$$
(18)
They are all evaluated at the last stable orbit. The values of the parameters adopted in this paper are summarized in Table 2.

The simulations are done in the region $`2r_SR60r_S`$, $`0\theta \pi /2`$, with $`120\times 120`$ mesh points. We assume axisymmetry with respect to $`z`$-axis and mirror symmetry with respect to the equatorial plane. We employ a free boundary condition, which waves, fluids and magnetic fields can pass through freely, at $`R=2r_S`$ and $`R=60r_S`$ (see Koide et al. 1999).

## 3 Numerical Results

### 3.1 Formation of Jet

We shall first discuss the results of the case A3 (Table 2; from now on referred to as the standard case). The model parameters are $`B_0=0.05`$, $`v_0=0.01`$, $`E_{th}=2.36\times 10^3`$, $`E_{mag}=1.68\times 10^3`$, and $`E_{rot}=5.36\times 10^2`$. Figure 3 shows the time evolutions of density. The stellar matter falls onto the central black hole at first. This collapse is anisotropic as a result of the effects of rotation and magnetic field. The accreting matter falls more slowly on the equatorial plane than on the rotational axis. The matter piles up on the equatorial plane, and a disklike structure is formed near the central black hole. At $`t/\tau _S=120`$, the disk density is about 2 orders of magnitude higher than the density of matter above the disk. Since the magnetic field is frozen into the plasma, it is dragged by the accreting matter and deformed to be more radial. At about $`t/\tau _S=60`$, a shock wave is formed near the central black hole and propagates outward anisotropically, because of the geometry of the magnetic field. A vortex pattern is generated inside the shock front, and moves outwardly with the propagating shock front, and produces a jetlike outflow. As this jetlike outflow moves up along the magnetic field lines, it becomes collimated. The simulation is terminated at $`t/\tau _S=175`$, at which time the jet is still propagating outward.

The jetlike outflow is generated behind the shock wave. Figure 4 shows the time evolutions of plasma beta ($`=P_{gas}/P_{mag}`$) and toroidal magnetic field for the standard case A3. The shock wave expands with twisted magnetic fields, and the outflow has high magnetic pressures. The jetlike outflow formed in the simulation is thus magnetically driven.

The magnetic field plays an important role also in the collimation. The plasma beta is initially low in the whole region, indicating that the magnetic pressure is dominant over the gas pressure (see Fig.(4)). Not only the pinching force of the toroidal magnetic field but also the geometry of the poloidal magnetic field (i.e. poloidal magnetic pressure) plays a crucial role.

### 3.2 Mechanism of Shock

We now investigate in more detail the mechanism of shock formation in the simulation. Figure 5 shows the time evolution of the spatial distributions of plasma beta ($`\beta =P_{gas}/P_{mag}`$) and toroidal magnetic field near the black hole in the standard case A2. At $`t/\tau _S=60.0`$, the accreting stellar matter rotates faster than earlier around the central black hole. As time goes on, the magnetic field is more and more twisted as a result of the differential rotation of accreting matter after the poloidal magnetic field is deformed to be more radial by the radial accretion. The toroidal magnetic field is amplified significantly near the central black hole. The amplified magnetic field expands outward as Alfvén waves and launches an outgoing shock wave. We can see that the expansion of the toroidal magnetic field takes place at the same position as the shock front. The plasma beta ($`\beta =P_{gas}/P_{mag}`$) inside the shock wave is low ($`\beta <1`$), implying that the shock wave is generated mainly by the magnetic field.

Figure 6 and 7 show the time evolutions of various physical quantities (density, gas and magnetic pressures, magnetic field, velocity, and forces) on the cylindrical surface with $`x/r_S=4`$. It shows the propagation of the outgoing shock. It is also found that when the shock wave is generated near the central black hole, the magnetic field is amplified by about a factor of 6 from the initial strength of the magnetic field near the black hole. The magnetic force $`F_{mag}`$ is comparable to the centrifugal force $`F_{cen}`$. These are about an order of magnitude higher than the pressure gradient force $`F_p`$. This confirms the claim that the magnetic force, together with the magnetocentrifugal force produces the shock wave rather than the pressure gradient force.

### 3.3 Properties of Jet

We discuss here the properties of the jet found in our simulation. Figure 8 shows the distributions of various physical quantities along the jetlike outflow at $`t/\tau _S=175`$. By this time, the shock wave has already passed through the upper boundary ($`z/r_S=30`$) of this figure.

It is easily seen in the velocity distribution that the jet has a mildly relativistic velocity, $`0.3`$ c (the poloidal velocity is $`0.1`$ c). It is larger than the sound velocity, as well as the Keplerian velocity. It also clearly exceeds the escape velocity at that point, and the jet is likely to get out of the stellar remnant. The toroidal velocity is larger than the poloidal velocity in the jet-forming region. The density of the jet is about twice as high as the density of the surrounding matter, $`3\times 10^8\mathrm{g}/\mathrm{cm}^3`$ if the density unit is $`\rho _0=10^{10}\mathrm{g}/\mathrm{cm}^3`$. The toroidal magnetic field is the dominant component of the magnetic field in the jet-forming region. The total magnetic field strength there is $`3.0\times 10^{13}\mathrm{G}`$.

On the other hand, we see in the energy distribution that the magnetic energy is comparable with the rotational energy and larger than the gravitational energy, which is consistent with the previous statement that the jet is accelerated by the magnetic and centrifugal forces.

Figure 9 shows the time variation of the mass flux of the jet and the accretion rate, the kinetic energy and the Poynting flux, and the jet and Alfvén wave components of the Poynting flux at $`z/r_S12`$ for the standard case A2. Here, the jet and Alfvén components are defined as

$`S_{EMjet}`$ $`=`$ $`{\displaystyle \frac{B_\varphi ^2}{4\pi }}v_z,`$ (19)
$`S_{EMAlf}`$ $`=`$ $`{\displaystyle \frac{B_\varphi B_z}{4\pi }}v_\varphi .`$ (20)
The Poynting flux is an order of magnitude larger than the kinetic energy flux. In the Poynting flux, the Alfvén wave component dominates over the jet components. This means that the Alfvén wave transports more energy outward. Using the simulation data, we can estimate that the kinetic energy of the jet, $`E_{jet}`$, is $`10^{49}\mathrm{ergs}`$ if the density of the jet is $`10^{10}\mathrm{g}/\mathrm{cm}^3`$. This is 2 orders of magnitude lower than the standard energy of GRBs ($`10^{51}\mathrm{ergs}`$). If Poynting flux is considered, the total energy becomes $`10^{50}\mathrm{ergs}`$, still less than the energy of GRBs.

### 3.4 Dependence on the Initial Magnetic Field Strength

The dependence of the jet properties on the initial magnetic field strength has also been investigated. Figure 10 shows the density and plasma beta distributions for case A3 in which the initial magnetic field is greater than the standard case A2 ($`B_0=0.07,v_0=0.01`$). The jetlike outflow is weaker and fainter than in the standard case A2. The maximum poloidal velocity of the jet is less than 0.05 c, although the jet velocity is still larger than the escape velocity. Thus, even in this case, a jet will emerge. The shock wave itself, as well as the magnetic twist, is also weaker.

The opposite is true for weaker magnetic fields. Shown in Figure 11 is model A1 ($`B_0=0.03,v_0=0.01`$). The jetlike outflow is clearly stronger than in the standard case A2. The maximum poloidal velocity of the jet becomes about 0.1c. The total velocity exceeds the escape velocity inside the shock. The magnetic twist is significantly stronger, and as a result the propagating shock wave is stronger as well.

Figure 12 summarizes these dependenices on the initial magnetic field strength. As the initial magnetic field strength increases, the jet velocity ($`v_z`$, $`v_\varphi `$) increases and the magnetic twist decreases. However, for stronger magnetic field ($`B_0>0.055`$) the jet velocity decreases with increasing $`B_0`$ and the magnetic twist still continues to decrease. In order to produce a strong shock, the magnetic field has to be twisted significantly so that it can store enough magnetic energy. If the initial magnetic field is strong, the magnetic field cannot be twisted significantly because Alfvén waves propagate as soon as the magnetic field is twisted a little bit. As a result, the jet velocity does not rise up and the magnetic twist remains weak. Therefore, weaker initial magnetic fields are favorable for a stronger jet.

Let us now consider the physical reason of the above results. In the Newtonian case, the time evolution of the toroidal magnetic field is given as

$$\frac{B_\varphi }{t}=\omega B_p,$$
(21)
where

$$\omega \frac{1}{t_{rot}}$$
(22)
and $`t_{rot}`$ is the rotation time scale. From this we obtain

$$\frac{B_\varphi }{B_p}\omega t.$$
(23)
In this situation, the comparing time scale is the propagation time scale of the Alfvén wave

$$t\frac{L}{v_{Az}}\frac{L\sqrt{4\pi \rho }}{B_z},$$
(24)
where $`L`$ is the typical length and $`v_{Az}`$ is the vertical component of the Alfvén wave velocity vector. Substituting equation (24) for equation (23), we have

$$\frac{B_\varphi }{B_p}\frac{\omega L\sqrt{4\pi \rho }}{B_z}\frac{1}{B_z}.$$
(25)
In this simulation, $`B_zB_0`$. Therefore,

$$\frac{B_\varphi }{B_p}\frac{1}{B_0}.$$
(26)
This can explain the dependence of the magnetic twist on the initial magnetic field strength, especially in the case of the Alfvén wave component (see Fig. 11c, dashed line).

The upward motion of the fluid behind the twist wave front is induced by the $`𝐉\times 𝐁`$ force. If we neglect other forces, the equation of motion for the fluid element in the z-direction becomes

$$\rho \frac{v_z}{t}\left(\frac{B_\varphi ^2}{8\pi }\right),$$
(27)
which can be rewritten as

$$v_z\frac{1}{\rho }\frac{t}{z}\left(\frac{B_\varphi ^2}{8\pi }\right).$$
(28)
In this situation the time scale is determined by the propagation time scale of the Alfvén wave. Thus, $`t/z1/v_{Az}`$. Since in the region near the wave front, $`B_pB_zB_0`$, equation (28) can be rewritten as

$$v_z\frac{1}{\rho }\frac{1}{v_{Az}}\left(\frac{B_\varphi ^2}{8\pi }\right)B_0^1.$$
(29)
This explains the dependence of $`v_z`$ on $`B_0`$ for $`B_0>0.055`$ in Figure 11a; i.e., as the initial magnetic field strength increases, the vertical component of the jet velocity decreases. However, the jet velocity has a limit when the initial magnetic field strength becomes even weaker. It can be understood as follows: the energy transport by the magnetic field is the Poynting flux

$$\frac{c}{4\pi }𝐄\times 𝐁\frac{B_\varphi ^2}{4\pi }v_z+\frac{B_zB_\varphi }{4\pi }v_\varphi .$$
(30)
The first term on the right-hand side of equation (30) is the jet component, whereas the second term is the Alfvén wave component. This equation also can explain the dependence of the jet velocity on the initial magnetic field strength. When the initial magnetic field strength becomes stronger, $`B_z`$ increases. This means that the Alfvén wave component takes away more energy than the jet component. Thus, a less magnetic energy is converted to the kinetic energy of the jet, so that the jet becomes weaker. When the initial magnetic field strength becomes weaker, the magnetic twist becomes stronger and $`B_\varphi `$ increases. Thus, much of the magnetic energy is converted to the kinetic energy of the jet. However, the converted magnetic energy has a limit. When the magnetic twist ($`B_p/B_\varphi `$) of the Alfvén wave becomes 1, it becomes a limit. This threshold is determined by the balance between the jet component and the Alfvén wave component

$$\frac{B_\varphi ^2}{4\pi }v_z=\frac{B_zB_\varphi }{4\pi }v_\varphi ,$$
(31)
which can be rewritten as

$$\frac{B_\varphi }{B_z}=\frac{v_\varphi }{v_z}1.$$
(32)
Thus, we can understand that when the magnetic twist in the Alfvén wave becomes 1, the jet velocity becomes maximum.

On the other hand, when the initial magnetic field is zero (case D), the shock wave is not generated and the jetlike outflow is not formed either. This clearly demonstrates that the generation of the shock wave, as well as the jetlike outflow, requires magnetic fields.

### 3.5 Dependence on the Initial Rotational Velocity

Next we show the dependence of the jet properties on the initial rotational velocity. Figure 13 shows the results (density and plasma beta distributions) of the fast rotation case B2 ($`B_0=0.05,v_0=0.015`$). The jetlike outflow is stronger and wider than in the standard case A2. The maximum poloidal velocity of the jet is about 0.1 c. The jet velocity exceeds the local escape velocity in the whole region. The magnetic twist is stronger accordingly.

The slow rotation case B1 ($`B_0=0.05,v_0=0.005`$) is shown in Figure 14. The jetlike outflow is weaker and fainter. The maximum poloidal velocity of the jet is less than 0.01 c. The jet velocity is barely larger than the escape velocity in some regions. The magnetic twist is weaker correspondingly.

Figure 15 shows the dependence of various physical quantities of the jet on the initial rotational velocity. As the initial rotational velocity becomes faster, the jet velocity ($`v_z`$, $`v_\varphi `$) becomes faster and the magnetic twist becomes stronger up to a certain value ($`v_0=0.01`$). For faster initial rotational velocity ($`v_0>0.01`$), the vertical component of the jet velocity ($`v_z`$) and the magnetic twist are almost constant and the toroidal component of the jet velocity ($`v_\varphi `$) becomes a little bit slower. If the initial rotation is sufficiently fast, the magnetic field is twisted significantly and stores enough energy to produce a fast jet. However, the stored magnetic energy has a limit. It is determined by the competition between the propagation time of Alfvén waves and the rotation time of the disk (or the twisting time). This can be understood by the same reason as the dependence on the initial magnetic field strength. The limit is determined when the magnetic twist of the Alfvén wave becomes 1. Provided with the results of this and the previous sections, the magnetic field strength is the dominant factor to determine the maximum jet velocity. It should be noted, however, that when the initial rotational velocity is zero (case C), the shock wave is not generated and the jetlike outflow is not formed either. Hence, the magnetic twist is crucial to produce the shock wave and the jetlike outflow.

## 4 Summary and Discussion

We have studied the generation of a jet from gravitational collapse of a rotating star with magnetic fields by using the 2.5-dimensional GRMHD simulation code. Our results are summarized as follows:

1. When the stellar matter falls onto the black hole, the collapse is anisotropic because of the effect of rotation and magnetic field and a disklike structure is formed near the black hole. A shock wave is launched near the black hole and propagates outward. The jetlike outflow is produced inside the shock. The maximum velocity of the jet is mildly relativistic ($`0.3`$ c). This result is consistent with the pseudo-Newtonian case (Proga et al. 2003). The kinetic energy of the jet is $`10^{49}\mathrm{ergs}`$.
2. The shock wave is generated mainly by magnetic and centrifugal forces. The frozen magnetic field is twisted and amplified by the rotation of stellar matter particularly near the black hole. It starts to expand outwardly with Alfvén waves and leads to the production of the shock wave.
3. The acceleration of the jet is also driven by the magnetic pressure and centrifugal forces. The jet is collimated in the course of propagation by the pinching force of the toroidal magnetic field and the geometry of the poloidal magnetic field.
4. The magnetic twist is large in the jet-forming region. The jet has a helical structure. This is similar to other astrophysical jets formed by MHD processes.
5. The Poynting flux is an order of magnitude larger than the kinetic energy flux. Alfvén waves as the Poynting flux transports more energy outward than the jet.
6. As the initial magnetic field strength increases, the jet velocity increases and the magnetic twist decreases. However, for stronger magnetic field the jet velocity decreases with increasing initial magnetic field strength and the magnetic twist goes on decreasing. The dependence of the jet properties on the initial rotational velocity is similar to the dependence on the initial magnetic field strength. As the initial rotational velocity becomes faster, the jet velocity becomes faster and the magnetic twist becomes stronger up to a certain value. For faster initial rotational velocity, the jet velocity and the magnetic twist are almost constant. The magnetic energy is converted to the kinetic energy of the jet. However, the converted magnetic energy has a limit. It is determined by the competition between the propagation time of Alfvén waves and the rotation time of the disk (or the twisting time). Hence the jet velocity is saturated at some point.

The maximum jet velocity in the current simulations was about $`0.3`$ c and the typical kinetic energy of the jet was $`10^{49}\mathrm{ergs}`$ . The jet is too slow and too weak for the jet of GRBs. However, the baryon mass in the jet is small. This gives good results for the GRB model. We have to consider other acceleration mechanisms. We expect that disk jets formed from an accretion disk become faster than the jet in the current simulations because they can release more gravitational energy and have enough energy to be applied to the GRB model. We believe that, if the current simulations can be continued for longer times, a disk-jet may be formed by the MHD processes. If the jet contains enough energy to convert the kinetic energy of the jet, then the break-out of the jet through the stellar surface is a possibility. When the jet goes through the stellar surface, the strong density gradient may accelerate the jet. In fact, some authors (Aloy et al. 2000; Zhang, Woosley, & MacFadyen 2003) have shown numerically that a significant acceleration of the jet occurs and the terminal Lorentz factor becomes as high as $`\mathrm{\Gamma }50`$. Since they used special relativistic hydrodynamics codes and neglected magnetic fields altogether, we think that it is important to simulate the propagation of the jet outside the stellar surface by a GRMHD code properly evaluating the importance of magnetic fields for the dynamics and further propagation of the jet.

On the other hand, our results may apply to baryon-rich outflows associated with the so-called failed GRBs. It is supposed to be a high baryon-load fireball attaining mildly relativistic velocities and not producing GRBs. Such failed GRBs may have event rates higher than successful GRBs (Woosley et al. 2002; Huang, Dai, & Lu 2002). Some failed GRBs can be observed as “hypernovae.” SN 2002ap may be such an event. It shows broad-line spectral features like the famous SN 1998bw (Kinugasa et al. 2002), yet it is not associated with a GRB. Kawabata et al. (2002) proposed a jet model for SN 2002ap consistent with their spectropolarimetric observations. Totani (2003) estimated that the jet should be moving at about 0.23 c in the direction perpendicular to the line of sight and that the kinetic energy of the jet was $`5\times 10^{50}\mathrm{ergs}`$ .

We have assumed a uniform global magnetic field in the simulations. However, this assumption may be unrealistic. From observations, it is inferred that rotating compact stars in general have a dipole-like magnetic field. Hence, it may be more likely that the rotating stars collapse with a dipole-like magnetic field. On the other hand, it is possible that the magnetic field near the central black hole is radial because of radial accretion of matter. We will study these possibilities in a forthcoming paper. There are a lot of works on astrophysical jets using non-relativistic MHD simulations with a dipole magnetic field configuration (e.g., Hayashi, Shibata, & Matumoto 1996; Goodson, Winglee, & Böhm 1997; Goodson, Böhm, & Winglee 199; Goodson & Winglee 1999; Romanova et al. 2002; Kato et al. 2003). These simulations show a high-velocity wind and lead the ejection of hot plasmoids. If the plasmoids are ejected intermittently, they will collide with each other, forming “internal shocks,” which may explain the time variability of GRBs (Shibata & Aoki 2003; Aoki, Yashiro, & Shibata. 2003).

Y. M. appreciates many helpful conversations on GRBs with S. Aoki, M. Uemura, and T. Totani. He also thanks K. Watarai, Y. Kato, K. Uehara, K. Nishikawa, and T. Haugboulle for useful discussions. This work was partially supported by Japan Science and Technology Cooperation (ACT-JST), Grant-in-Aid for the 21st Century COE “Center for Diversity and Universality in Physics” and Grants-in-Aid for the scientific research from the Ministry of Education, Science, Sports, Technology, and Culture of Japan through 14079202, 14540226 (PI: K. Shibata) and 14740166. The numerical computations were partly carried out on the VPP5000 at the Astronomical Data Analysis Center of the National Astronomical Observatory, Japan (yym17b), and partly on the Alpha Server ES40 at the Yukawa Institute for Theoretical Physics of Kyoto University, Japan.


text:
# Mechanical versus thermodynamical melting in pressure–induced amorphization: the role of defects

## I Introduction

Since the first reports on pressure–induced amorphization (PIA), there has been an increasing interest in the phenomenon that has been observed in a large class of systems, both upon compression and decompression (see Ref.  for a review and references on particular materials). The experimental interest stems from the fact that PIA occurs in some extremely widespread and important materials, namely water mishima  and quartz and its polymorphs silica , and also because PIA provides a novel route to the synthesis of amorphous materials, in addition to the traditional technique of quenching from the melt. ¿From a basic point of view, there are controversies about the mechanism of PIA and the nature of the amorphous phase, in particular regarding whether it is equivalent or not to the one obtained by quenching sio2anis ; iniguez . Since in most cases crystalline phases that undergo PIA show reentrant melting in the $`P`$$`T`$ plane, it has been suggested mishima  that PIA is nothing but melting at temperatures below the glass temperature of the supercooled fluid. Other studies (specially numerical ones softening1 ; softening2 ; grande ; sciortino ) have instead emphasized the relationship between PIA and mechanical instabilities. It has been observed in fact, that in many cases PIA is triggered by the softening of a shear phonon mode binggeli94 . In addition, some memory effects memory  and anisotropic properties sio2anis , although controversial contro , show that many times what appears to be an amorphous phase preserves within its structure signatures of the parent crystalline phase.

In view of the broad phenomenology briefly stated above, it is highly desirable to take advantage of model systems in which amorphization can be studied in a transparent way, allowing to look in detail into the mechanisms and characteristics of the transformation. This motivates the present work, in which we study a simple two–dimensional system of identical point–like particles interacting through a specially devised two–body potential. The simplicity of the model allows us to study large systems with some amount of defects, and observe directly the crucial role they play in the transformation.

In the model we study, there are a few different crystalline ground states depending on the applied pressure. We present here the results of the evolution of the most compact structure (stable at the highest pressures) upon pressure release decompress . As we will see, PIA at temperatures in which particle diffusion is negligible is always related to mechanical instabilities. For perfect lattices, PIA is reflected in the softening of a shear phonon mode. This leads to local distortions that produce the destabilization of new vibrational modes, leading to a cascade of instabilities grande , which drives the system toward an amorphous structure. However, in the presence of defects, localized vibrational modes exist, that may become unstable before any extended vibrational mode does. This favors the nucleation of the amorphous phase at pressures much closer to the thermodynamic equilibrium value between the crystalline and amorphous phases.

The work is divided as follows. In Section II we present the model. The results at zero temperature are contained in Section III. In Section IV we give evidence that the disordered samples obtained can be called trully amorphous. In Section V we present results at finite temperatures, and Section VI contains some discussion and conclusions.

## II The model

An isotropic, purely repulsive interparticle potential is used, with a strict hard core at a distance $`r_0`$ plus an almost triangular repulsive shoulder. The pairwise interaction potential $`V(r)`$ between two particles separated a distance $`r`$ is given by jagla2001

$`V(r)`$ $`=`$ $`\mathrm{}\text{ for }r<r_0`$
$`V(r)`$ $`=`$ $`\epsilon _0\left[1.22.8125(r/r_01.08)^2+{\displaystyle \frac{0.008}{r/r_01}}\right]^2`$ (1)
$`\text{for }r_0<r<r_1`$
$`V(r)`$ $`=`$ $`0\text{ for }r>r_1`$
where $`r_1=1.17315r_0`$, and $`r_0`$ sets the length unit. The potential is plotted in the inset of Fig. 1. This kind of potential has been previously used to systematize the anomalous properties of tetrahedrally coordinated materials jagla98 ; jagla2001 ; sadr . Then we expect they are also appropriate to study amorphization under pressure, since this phenomenon occurs for most of these materials. Its crucial characteristic is the existence of two possible equilibrium distances between particles.

The system is simulated by standard molecular dynamics in the NVT ensemble with periodic boundary conditions. At a given volume, the quantities of interest are evaluated, and the volume is changed in steps of the order of $`0.01\%`$ by rescaling all coordinates of the particles and size of the simulation box. Temperature is fixed by rescaling the velocities of the particles whenever necessary. As the amorphization process implies the existence of mechanical instabilities, we typically observe that kinetic energy tends to increase during amorphization. In a real situation this energy transforms into heat. Here we simply eliminate it by the mentioned rescaling procedure. Pressure is calculated by a direct evaluation in terms of the interparticle forces. We choose to model the system at constant volume in order to survey all regions of the volume–pressure curve, including those that would be unstable in constant pressure simulations. The results to be presented correspond to a two–dimensional system, to facilitate a direct visualization of the particle configurations. We should mention however, that the same phenomenology was observed in three–dimensional samples with the same interparticle potential.

## III Zero temperature results

We take as the starting configuration the one corresponding to the triangular lattice with lattice parameter $`r_0`$, which is stable at high pressures. In Fig. 1 we present the results of simulations at $`T=0`$ in a system of 2800 particles. The thin continuous line corresponds to the $`v`$$`P`$ relation assuming the system remains always triangular. The reentrance of this line is a consequence of the particular form of the interaction potential. The numerical results for a perfect lattice (full circles) follow this line up to some maximum volume, at which they abruptly depart from it. This is the pressure in which a phonon with vanishing energy first appears in the triangular structure (see Fig. 2). Analytical evaluation shows that this instability occurs when $`3V^{}(a)/a+V^{\prime \prime }(a)=0`$, where primes denotes derivatives of the potential, and $`a`$ is the lattice parameter. This expression is valid as long as the interaction between next nearest neighbors is zero, as it is in the present case. The critical volume $`v^c`$ and pressure $`P_0^c`$ are indicated by dotted lines in Fig. 1, and they are fully compatible with the numerical results. The unstable phonons turn out to be shear phonons with the wave vector $`𝐤`$ oriented perpendicularly to one of the three most compact directions in the lattice. Note that for the interparticle potential we use, the longitudinal phonon branch along these directions is of the type $`\mathrm{sin}(k)`$, and all shear phonons –irrespective of the absolute value of $`k`$– become zero energy at the same point. Then, in the present case the instability pressure is macroscopically signaled by the vanishing of the stress dependent shear modulus of the material gprimap ; zhou96 , $`\mu =(C_{11}C_{22})/2P=0`$.

At $`v^c`$ there is a sharp and abrupt increase of pressure singular . Snapshots of the particle configurations give clues of what happens in the system. In Fig. 3(a) we see the triangular configuration just before the instability. Figure 3(b) shows that after changing the volume a very small quantity, an instability has propagated in the system. We stress that the evolution of the system from the configuration in Fig. 3(a) to that in 3(b) is triggered by a very small volume change, and is just the steepest descend evolution toward a local energy minimum of the energy landscape. In this evolution, pressure recovers to a higher value indicated as $`P^{}`$ (dashed line) in Fig. 1. This value is roughly maintained tbp  upon further volume increase up to $`v/r_0^21.27`$, which is close to the volume at which the disordered regions have taken over the whole system. If volume is increased further, pressure decreases noticeably.

The thermodynamic equilibrium pressure $`P^e`$ at $`T=0`$ between the crystalline and disordered structures is indicated by the arrow in Fig. 1. The value of $`P^e`$ was calculated by enthalpy evaluations of the triangular lattice and the amorphous structure obtained by quenching from the melt. Thus, $`P^e`$ represents the natural extrapolation of the melting line to zero temperature. We see that $`P^e`$ is close to (though a bit higher than) $`P^{}`$.

It must be noticed that constant pressure simulations would have shown at $`P_0^c`$ an abrupt transition between a completely ordered and a completely disordered system, while $`P^{}`$ would have remained hidden. However, $`P^{}`$ has a clear physical meaning as the transition pressure once the disordered phase has been nucleated. This suggests that if nucleation centers are present in the system, $`P^{}`$ will be experimentally accessible as the actual transition pressure. In fact, Fig. 1 shows also the $`v`$$`P`$ evolution for a system with a single vacancy and for a system with grain boundaries. In the former case, the instability occurs before $`v^c`$ is reached, and in the latter case it even occurs without any pressure reentrance nota . Note that the three curves tend to coincide after the first destabilization of the original lattice. Two snapshots of the system for the poly–crystalline case \[Fig. 3(e) and (f)\] show how the disordered phase grows from grain boundaries, that provide nucleation centers for the transformation. These results are a clear indication that defects are very effective in lowering the energy barrier for the transition, and also explain why, in experiments, transformation pressures close to the thermodynamic values are usually observed. The lowering of the energy barrier due to defects and its relation to thermodynamic melting has been discussed by Mizushima et al. mizu94  in the context of crystal–crystal pressure–induced transformations. Note that in the case in which the transformation is triggered by defects, the phonon spectrum of the system (and in particular the elastic constants) gives no indication of the instability that is about to occur. This is particularly obvious in the case of a single vacancy: a single defect cannot modify the phonon spectrum of an infinite sample, and then its effect it is not seen in the elastic constants, but it produces a finite change in the value of the critical pressure $`P_0^c`$single .

It was already mentioned that for the interaction potential used, the whole branch of shear phonons become zero energy at $`P_0^c`$. We want to emphasize however that amorphization it is not related to this degeneracy, as it can also be triggered by a single phonon becoming unstable. In fact, we did simulations with a modified potential in which an attractive term $`V_a`$ reaching up to second nearest neighbors was included. The actual potential used is the one of the previous section plus an attractive term given by

$`V_a(r)`$ $`=`$ $`2\epsilon _0\left(rr_1\right)^2\text{ for }r_1<r<\left(r_1+2r_0\right)/2,`$
$`V_a(r)`$ $`=`$ $`2\epsilon _0\left[\left(r2r_0\right)^2\left(r_12r_0\right)^2/2\right]`$ (2)
$`\text{for }\left(r_1+2r_0\right)/2<r<2r_0.`$
The instability corresponds now to the vanishing of the energy of edge zone shear phonons only, as indicated in Fig. 2. The amorphization of the system occurs precisely at the point where this phonon becomes zero energy, and it was observed to be similar to the previous case where the full shear phonon branch becomes unstable at the same pressure. This rules out an amorphization mechanism in which at the instability, arbitrary combinations of the unstable phonons generate disordered movements of the atoms, and favors a cascade mechanism as described in the discussion section below.

## IV Comparison between decompressed and quenched amorphous

The claim that the disordered structure we obtain upon decompression is really amorphous is supported by the following facts. We find no evidence of orientational order in the structure factor $`S(𝐤)`$, even in the case in which we start with a mono–crystalline sample. Moreover, a quantitative comparison with the structure factor of a system quenched from the liquid (prototypic of an amorphous sample) shows that they are indistinguishable (Fig. 4). In addition, we have failed to find any single systematic difference between the two amorphous.

## V Finite temperature effects

If temperature is not strictly zero, even the perfect lattice amorphizes before the ultimate mechanical instability limit is reached. Simulations at finite temperature with perfect samples show that the amorphization occurs at a temperature dependent critical pressure $`P^c(T)`$ (such that $`P^c(T=0)=P_0^c`$). The temperature dependence of the critical pressure originates in the fact that a thermally activated process may destabilize the soft phonon before it actually becomes zero energy.

It is instructive to see in some more detail how this destabilization occurs. Consider the normal modes (phonons) of the system. We will study the case in which only a single phonon with wave vector $`𝐤`$ perpendicular to one of the high density planes has non-zero amplitude. We expand the energy of the lattice in powers of the amplitude of the phonon oscillation, going one order beyond the harmonic approximation. Let us call $`\delta _l`$ and $`\delta _t`$ the longitudinal and transversal oscillation amplitude of that phonon. As we will see there is a nontrivial coupling between the two at the instability.

Taking into account that our potential produces interactions only among first neighbors in the lattice, the energy per particle $`e_0`$ of the system up to third order in the amplitudes can be written after some lengthy but straightforward calculation as

$`e_0`$ $`=`$ $`{\displaystyle \frac{1}{8}}\mathrm{sin}^2\left({\displaystyle \frac{\sqrt{3}ka}{4}}\right)`$
$`\times \left[\delta _t^2\left({\displaystyle \frac{3V^{}}{a}}+V^{\prime \prime }\right)+\delta _l^2({\displaystyle \frac{V^{}}{a}}+3V^{\prime \prime })\right]`$
$`+{\displaystyle \frac{\sqrt{3}}{8}}\mathrm{sin}^3\left({\displaystyle \frac{\sqrt{3}ka}{4}}\right)\left[\delta _t^2\delta _l\left({\displaystyle \frac{4V^{}}{a^2}}+{\displaystyle \frac{V^{\prime \prime }}{a}}+V^{\prime \prime \prime }\right)\right].`$
The first line in (V) is the harmonic contribution, the second is a cubic term in the displacements (there is an additional cubic term proportional to $`\delta _l^3`$, but this can be shown to be not relevant for the analysis below). The system becomes unstable when the coefficient of $`\delta _t^2`$ vanishes, providing again the condition $`3V^{}/a+V^{\prime \prime }=0`$, which defines the $`T=0`$ values $`P_0^c`$ and $`v^c`$. Close enough to this instability point, $`3V^{}/a+V^{\prime \prime }`$ is proportional to $`PP_0^c`$ and to lowest order all other coefficients can be considered to be constants. Then generically, the energy can be written as

$`e_0`$ $`=`$ $`\mathrm{sin}^2\left({\displaystyle \frac{\sqrt{3}ka}{4}}\right)\left[A(PP^c)\delta _t^2+B\delta _l^2\right]`$ (4)
$`+`$ $`\mathrm{sin}^3\left({\displaystyle \frac{\sqrt{3}ka}{4}}\right)\left[D\delta _t^2\delta _l\right].`$
The existence of a $`\delta _t^2\delta _l`$ term is very important, as it indicates that for $`P`$ slightly larger than $`P_0^c`$ there is a saddle point (actually two, one with positive and the other with negative $`\delta _t`$) close to $`\delta _t=0`$, $`\delta _l=0`$, than can be determined requiring stationarity of (4). The result is

$`\delta _t^S`$ $`=`$ $`\pm {\displaystyle \frac{\sqrt{2AB(PP^c)}}{D\mathrm{sin}(\frac{\sqrt{3}ka}{4})}},`$
$`\delta _l^S`$ $`=`$ $`{\displaystyle \frac{A(PP^c)}{D\mathrm{sin}(\frac{\sqrt{3}ka}{4})}},`$ (5)
where the upperscript indicates the values at the saddle. If the system reaches this saddle it can escape from the local minimum at the origin, namely, the system destabilizes. The energy barrier $`h`$ for this process is obtained by reinserting (5) into (4). The result is:

$$h=\frac{4A^2B(PP^c)^2}{D^2}.$$
(6)
Note that the barrier is exactly the same whatever the wave vector of the phonon considered (this is no longer true if further neighbors interactions are included, but the possibility to escape through the jump of a barrier remains). The instability mechanism is then driven by the $`\delta _t^2\delta _l`$ term in the energy, what implies a not trivial coupling of transverse modes (the one actually having vanishing frequency) and longitudinal ones. At finite temperatures the barrier $`h`$ can be surmounted. Nucleation theory tells that escape time $`\tau `$ is proportional to $`\mathrm{exp}(h/T)`$. Assuming the preexponential factor is a constant $`\tau _0`$ (this is certainly not true, but dependences on $`\tau _0`$ become weak in the final result) we obtain that if a time $`t_0`$ is given, the system will overcame the energy barrier if $`t_0\tau `$, and from here we obtain the formula for the temperature necessary to escape a given barrier, namely:

$$T>h\mathrm{ln}^1(t_0/\tau _0).$$
(7)
For practical purposes the logarithmic factor can be usually taken to be approximately 0.1. As $`h(PP_0^c)^2`$, we obtain that the critical pressure $`P^c(T)`$ increases as $`T^{1/2}`$ at finite temperatures, a behavior that is well reproduced in the simulations (Fig. 5). The $`P^c(T)`$ line is then the pressure amorphization line for a sample without defects, and can be properly interpreted as the ‘mechanical melting’ line of this system. In the absence of defects it marks the limit on which the crystalline phase destabilizes (purely mechanically at $`T=0`$, and by a thermally activated process at $`T0`$). Note however, that the elastic constants of the system, in particular the shear modulus $`\mu `$, do not extrapolate to zero at $`P^c(T)`$, except at $`T=0`$ (see Fig. 6). In fact, the numerical evaluation by means of fluctuation formulas zhou96 ; constelast  shows that the effect of temperature in the elastic constants is very small. At $`T=0`$ the stress dependent shear modulus $`\mu `$ vanishes at $`P_0^c`$, as we already knew from the phonon dispersion relation, but at $`T0`$ extrapolation is consistent with a vanishing of $`\mu `$ very close to $`P_0^c`$, and certainly not at $`P^c(T)`$. The thermally activated process that at finite temperatures leads to the instability at the higher pressure $`P^c(T)`$ is a rare event that it is enough to occur once to completely destabilize the system. This process is not captured in the value of the elastic constants.

Up to here the results for samples without defects. For defective samples the amorphization line $`P^{}(T)`$ is only weakly dependent on $`T`$. In addition, it is very close to the extrapolation of the melting line obtained from simulations at higher temperatures (see Fig. 5). $`P^{}(T)`$ and $`P^c(T)`$ are well different from each other, particularly at very low temperatures. Defects are responsible for this difference, and provide the link between the ‘mechanical’ and ‘thermodynamical’ scenarios for pressure amorphization.

## VI Discussion and Conclusions

The experimental observation that in some materials PIA occurs roughly along the extrapolated melting line has motivated the suggestion mishima  that ‘thermodynamic melting’ is the underlying driving force for amorphization. Our results show that amorphization at zero temperature is always a phenomenon related to mechanical instabilities. In perfect lattices, mechanical instabilities reflect necessarily in the phonon spectrum of the material, and amorphization can be easily related to a mechanical process. However, for defective lattices, the mechanical instability that triggers amorphization is associated with localized vibrational modes around defects, and this instability occurs typically without any noticeably signature in the phonon spectrum, i.e., the mechanical nature of the process is more subtle. We emphasize that a single defect is sufficient to trigger the amorphization of a whole macroscopic samplesingle .

How this first instability is able to transform the highly ordered original lattice into a fully amorphous configuration is not easy to understand, specially considering that this evolution is fully deterministic from a mechanical point of view, and can be described as the steepest descendant path of the configuration point of the system onto its energy landscape. A deep understanding of this process will surely shed light also onto the very definition of what an amorphous material is, typically characterized by what it lacks, instead of what it possesses. Here we make only the following considerations grande . The first instability of the lattice (at $`T=0`$) can always be described as the vanishing of the frequency of one of the normal modes of the system. This normal mode may be a phonon for a perfect lattice, or a localized mode for a defective lattice. When the coordinate characterizing the destabilizing mode grows, it typically goes beyond the applicability limit of the linear theory, and the whole mechanical stability of the system has to be reanalyzed. Then, although the analysis of Section V is appropriate to determine the barrier for the system to scape from its metastable equilibrium, it cannot be continued when the system actually overcomes the saddle without including the interaction with other modes. For instance, it occurs in our system that once the coordinate of the unstable phonon starts to increase, it couples to other phonons and the results is that a shuffle of two pieces of the material develops. In this way, the instability generated by an extended object (the unstable phonon) produced a localized (one–dimensional) defect on the lattice.fonon . But it turns out that, in the present case, the new configuration of the lattice is not stable either, and a new distortion occurs spatially close to the region already distorted (see Fig. 7). In the present case of amorphization, it is likely that a cascade of these processes occurs that leads to an amorphous structure. The spatially localized nature of the process reflects in the fact that for partially transformed samples, transformed and untransformed regions are spatially separated, as seen in Fig. 3. We however emphasize that another possibility could be that after the first phonon instability and the first shuffle, the system reaches a mechanically stable configuration, then further increases of the volume will produce independent shuffles in the samples, along the three directions related by symmetry. If this is the case, then this single phonon will drive a transformation (which has the character of a martensitic transformation) to a new (poly–)crystalline phase. This and other intermediate cases, that in our model occur for different parameters, will be discussed separately.

In summary, we have arrived at a picture of the mechanism of PIA, reconciling mechanical and thermodynamical melting in a unifying scheme. The present model exhibits a one–step crystal–to–amorphous transition to a truly amorphous state that does not differ from the disordered structure obtained by quenching from the melt. Once the transformation is triggered, it is followed by a cascade process that produces the amorphous structure. For a perfect lattice, and vanishingly small temperatures, the transformation takes place through a shear phonon instability at a pressure far apart the thermodynamical equilibrium pressure. For defective lattices, the mechanical instability shows up by the destabilization of localized vibrational modes, close to the thermodynamical equilibrium pressure. In either case, the transformation is always triggered by mechanical instabilities. Temperature can produce the transition to occur before the mechanical stability limit is reached, due to a process of thermally activated jump over a barrier. In real samples with appropriate defects, experimental transition pressures will be seen close to the extrapolation of the melting line. These facts reconcile different claims in the literature regarding the nature of the PIA of materials. We remark that, although we presented results for a two–dimensional model system, preliminary results in three–dimensional samples qualitatively agree with the phenomenology stated above, and will be presented elsewhere. We finally stress that these considerations have potential applications in the study of melting mechanisms of crystalline materialsjin .

## VII Acknowledgments

We thank N. Binggeli for useful comments on the manuscript. S. B. is financially supported by CONICET (Argentina).