text:
# Three-dimensional MHD simulations of radiatively cooling, pulsed jets

## 1 Introduction

The bright Herbig-Haro (HH) objects associated with low-mass young stellar objects, are radiating shock fronts immersed in collimated optical jets that may extend from a few 1000 AU (e.g., the microjets in the Orion Nebula; Bally, O’Dell & McCaughrean 2000), to very large parsec-length scales, like the $`giant`$ HH jets (Heathcote et al. 1998; Bally & Devine 1994; see also Reipurth 1999 for a catalogue of such objects).

Currently, the most accepted model for the production of the discrete HH knots is by time variability in the driving sources of the jets. Their bow shock morphology, high spatial velocity, and symmetric pattern with respect to the star (e.g., Reipurth 1989a, 1989b; Bürke, Mundt & Ray 1988; Zinnecker, McCaughrean & Rayner, 1997, 1998) all indicate that the HH shocks generally arise through the steepening of velocity fluctuations in the underlying, supersonic outflow. Strong support for this conjecture has been given by theoretical studies which have confirmed that traveling shocks created in this manner reproduce the essential properties of the observed knots (e.g., Raga et al. 1990; Raga & Kofman 1992; Kofman & Raga 1992; Stone & Norman 1993; de Gouveia Dal Pino & Benz 1994; Völker et al. 1999; de Gouveia Dal Pino 2001). 111We note that although pinch-modes of the Kelvin-Helmholtz (K-H) were originally suggested as the source of the knotty pattern in HH jets, it has been found that these modes have their strength highly reduced in the presence of radiative cooling (e.g., Blondin, Fryxell and Königl 1990; de Gouveia Dal Pino & Benz 1993). Besides, they typically lead to a too fast periodicity to match the observations (e.g., Stone, Xu, and Hardee 1997). Nevertheless, nonaxisymmetric K-H modes appear to be responsible for the gentle wandering or wiggling observed in many jets (e.g., de Gouveia Dal Pino, Birkinshaw, & Benz 1996; Stone, Xu, & Hardee 1997).

Magnetic fields seem to play a fundamental role in several phases of the jet launching and propagation processes. They seem to be relevant to: (i) operate locally in an accretion disk, in order to regulate the angular momentum process (Balbus & Hawley 1998) and the accretion rate; (ii) launch outflows through magneto-centrifugal forces from the disk itself (e.g., Blandford & Payne 1982; Kudoh, Matumoto & Shibata 1998; Königl & Pudritz 2000), or from the disk-star boundary region (the so-called X-point; e.g., Shu et al. 2000); and (iii) collimate the outflow once launched (e.g., Ouyed, Pudritz & Stone 1997). On the larger scales, since the outflowing material does not seem to be highly resistive (e.g., Frank et al. 1999), the magnetic field may also play a relevant role on its dynamics and evolution. Previously, Morse et al. (1992) have inferred an upper limit for the magnetic field in the ambient medium upstream the HH34 bow shock of about $`1020\times 10^6`$ G. More recently, Ray et al. (1997), based on polarization measurements, have detected magnetic field strengths of the order 1 G at a distance of few tens of AU from T Tau S. Using magnetic flux conservation, this latter result would give magnetic fields of the order of $`10^3`$ G for a toroidal field configuration, and $`10^6`$ G for a longitudinal field at distances $`0.1`$ pc. For typical jet parameters, this would imply a plasma parameter $`\beta =p_{j,gas}/(B^2/8\pi )\mathrm{\hspace{0.17em}10}^3`$ to $`10^3`$ at those distances. Although Ray et al. have concluded that the high values measured for B at $`10`$ AU probably come from regions of strong magnetic field amplification (such as behind shocks), the figures above suggest that magnetic fields must play some role in the overall jet dynamics and propagation even on the outer scales where local amplification (particularly behind the leading and the internal bow shocks) is expected to occur (e.g., Cerqueira, de Gouveia Dal Pino & Herant 1997; hereafter, Paper I).

Lately, there has been increasing focusing on MHD studies of the propagation of overdense, radiative cooling jets in an attempt to look for possible signatures of magnetic fields on the large scales of the HH outflows (see, e.g., de Gouveia Dal Pino & Cerqueira 2001 for a review). Frank et al. (1997, 1998), for example, in axisymmetric studies of steady-state jets with toroidal magnetic fields (and $`\beta =8\pi p_{th}/B^2<1`$), have reported the formation of pinch modes driven by magnetic tension, and the appearance of a magnetically confined nose cone structure of compressed gas between the bow shock and the Mach disk at the jet head.

In previous investigations (de Gouveia Dal Pino & Cerqueira 1996; Paper I; Cerqueira & de Gouveia Dal Pino 1999, hereafter Paper II), we have employed a three-dimensional Smoothed Particle Hydrodynamic (SPH) code that has been implemented to compute the effects of magnetic fields to investigate steady-state jets. Assuming different initial magnetic field configurations (in approximate equipartition with the gas) permeating both the jet and the ambient medium, we have found that there is an increase in the jet collimation in comparison to purely hydrodynamic calculations, which is caused by the amplification and reorientation of the magnetic fields behind the shocks (particularly in the presence of a helical field). Further, we have found that the presence of a helical magnetic field may inhibit the fragmentation of the dense shell formed by cooling and Rayleigh-Taylor instability at the head of the jet (Paper I), while a longitudinal field retains the fragmentation. This result could be an indication that the latter configuration could dominate near the jet head as the observed jets have a clumpy structure there. We have also evidenced from the simulations the development of the MHD K-H modes, but the weakness of the pinch modes makes it doubtful that they could play an important role in the production of the bright knots (Paper II). Where there is overlap, these results have been confirmed by recent axisymmetric calculations (e.g., Gardiner et al. 2000; Stone & Hardee 2000; O’Sullivan & Ray 2000). In particular, Gardiner et al. (2000), assuming longitudinal magnetic fields (with $`\beta =0.1\mathrm{\hspace{0.17em}10}^7`$), have found that the global characteristics of the flow are not strongly affected by the $`𝐁`$-field strength, and a predominantly axial field tends to increase the collimation and order, and inhibit instabilities in the flow, but the presence of internal pulses increases the likelihood of magnetic reconnection.

The 2-D calculations of steady-state and pulsed MHD jets of Stone & Hardee (2000; hereafter, SH) propagating into unmagnetized ambient media have also shown that the magnetic fields affect the fragmentation of the shell at the jet head, and a strong toroidal field peaking near the jet surface develops a nose cone which should be unstable in three-dimensions. They have also found that the radial hoop stresses due to the toroidal field confine the shocked jet material in the internal pulses resulting in higher densities in the pulses which are strongly peaked towards the jet axis in comparison to purely hydrodynamic calculations. A similar line of numerical studies has also been conducted by O’Sullivan & Ray (2000; hereafter OR). Most of their results essentially confirm previous hydrodynamic (e.g., Blondin, Fryxell, & Königl 1990; de Gouveia Dal Pino & Benz 1993, hereafter GB93; Stone & Norman 1994) and the MHD results above. However, unlike these works, they have assumed jets with initially lower densities with respect to the ambient medium (corresponding to a density ratio $`\eta `$ = 1), and which are highly over-pressurized with respect to the ambient medium. This causes the development of more complex cocoons surrounding the beams and the formation of crossing shocks which help to refocus the beam and the internal pulses. In particular, they find that the hoop stresses associated with toroidal fields can cause the $`disruption`$ of the internal knots, even for $`\beta `$1.

In a recent study, we have presented the first results of 3-D MHD simulations of the early evolution of pulsed jets (Cerqueira & de Gouveia Dal Pino 2001, hereafter, Paper III). In the present paper, we attempt to extend these previous investigations by exploring in three-dimensions the role played by three different initial magnetic field configurations on the evolution and emission structure of radiatively cooling, pulsed jets, considering a more extensive range of parameters. As before (Paper I, II and III), together with a baseline non-magnetic calculation, we employ a modified version of our 3-D SPH magnetized code and consider a longitudinal and a helical field topology (in approximate equipartition with the gas), both permeating the jet and the ambient medium. Also, in order to make a closer comparison with previous studies, a third geometry involving a purely toroidal field permeating the jet is considered.

The paper is planned as follows. In §2, we will briefly outline the numerical technique and setup. In §3, we present the results of the simulations for pulsed jets; in §4, we compare the results with previous 2-D MHD calculation and discuss some observational implications of them, and present our final remarks and conclusions.

## 2 Numerical method

We solve the equations that describe the evolution of a fluid using the Smoothed Particle Hydrodynamic (SPH) technique (e.g.; Benz 1990; Monaghan 1999). SPH is a Lagrangean scheme that makes use of particles to follow fluid parcels, and the evaluation of a given physical quantity is given by an interpolation scheme, using for this a kernel function ($`W`$). The kernel has an adaptative smoothing length ($`h`$) value, which defines the space volume of the SPH particles. Like any other physical quantities, magnetic fields can also be handled with SPH. In order to take into account the presence of the magnetic field, we have modified a pure HD, SPH code. In Paper II, we give the details of such modification, as well as the results of some tests in the MHD limit. The application of the SPH technique to pure hydrodynamic jet simulations can be seen in several papers (e.g., GB93; de Gouveia Dal Pino & Benz 1994; Chernin et al. 1994; de Gouveia Dal Pino & Birkinshaw 1996; de Gouveia Dal Pino, Birkinshaw, & Benz 1996; de Gouveia Dal Pino 1999, 2001), so that our code has been continuously tested in both HD and MHD regimes. Such an MHD, SPH based code is often referred as “SPMHD” (e.g., Stellingwerf & Peterkin 1990; Meglicki 1994).

The solved system of MHD equations, in its ideal approximation, is given by:

$$\frac{d\rho }{dt}=\rho 𝐯$$
$`(1a)`$
$$\frac{d𝐯}{dt}=\frac{p}{\rho }+\frac{1}{4\pi \rho }(\times 𝐁)\times 𝐁$$
$`(1b)`$
$$\frac{du}{dt}=\frac{p}{\rho }(𝐯)$$
$`(1c)`$
$$\frac{d𝐁}{dt}=𝐁(𝐯)+(𝐁)𝐯$$
$`(1d)`$
where $`\rho `$ is the density; $`𝐁`$ is the magnetic field; $`u`$ is the specific internal energy and $``$ is the radiative cooling rate. As in our previous works (see Paper II and references therein), we have adopted the cooling function given by Katz (1989) for a gas of cosmic abundances cooling from $`T10^6`$K to $`10^4`$K. An ideal equation of state is used to close the above system of equations:

$$p=(\gamma 1)\rho u$$
$`(1e)`$
with $`\gamma =5/3`$.

### 2.1 Initial and boundary conditions

Our computational domain is a rectangular box with dimensions $`30R_j`$ x $`30R_j`$, and $`10R_j`$ y,z $`10R_j`$, where $`R_j`$ is the initial jet radius (which is also the code distance unit). The Cartesian coordinate system has its origin at the center of the box and the jet flows through the x-axis. The jet is continuously injected into the bottom of the box \[at $`𝐫=(30R_j,0,0)`$\]. Inside the box, the SPH particles are initially distributed on a cubic lattice. An outflow boundary condition is assumed for the boundaries of the box. The particles are smoothed out by a spherically symmetric kernel function of width $`h`$, and the initial values of $`h`$ were chosen to be $`0.4R_j`$ and $`0.2R_j`$ for the ambient and jet particles, respectively, so that we have up to 400,000 SPH particles at the beginning of the calculation.

We adopt a sinusoidal velocity profile for the pulsing jet at the inlet:

$$v_o(t)=v_j[1+A\mathrm{sin}(\frac{2\pi }{P}t)]$$
$`(2)`$
where $`v_j`$ is the mean jet speed, and $`A`$ and $`P`$ are the amplitude and the period of the velocity oscillations, respectively. We have chosen two periods for our models, namely, $`P=1`$ and $`P0.5`$ (in code units, where the unit of time in our calculations corresponds to the ratio $`t_d=R_j/c_a38`$ years, for the set of parameters adopted here; see below). The number density in the flow at injection, $`\rho _j`$, is assumed to be constant, which implies that $`\dot{m}=\rho _jv_j(t)`$ is not constant. These models are equivalent to the constant density injection models recently investigated analytically by Cantó et al. (2000).

We consider different initial magnetic field profiles. As in papers I, II and III, we adopt 1) an inilial constant longitudinal magnetic field \[$`𝐁=(B_0,0,0)`$\], both inside and outside the jet beam; and 2) a helical force-free field that also extends to the ambient medium which is described by the following equations (see Todo et al. 1993, and Fig. 1 of Paper II):

$$B_r=0$$
$`(3a)`$
$$B_\varphi (r)=B_0\left[\frac{0.5Cdr^2}{(r+0.5d)^3}\right]^{1/2}$$
$`(3b)`$
$$B_x(r)=B_0\left[1\frac{Cr^2(r+d)}{(r+0.5d)^3}\right]^{1/2}$$
$`(3c)`$
where $`r=\sqrt{y^2+z^2}`$ is the radial distance from the jet axis and the $`C`$ and $`d`$ constants are given by 0.99 and $`3R_j`$, respectively. In these equations, $`B_0`$ is the maximum strength of the magnetic field and corresponds to the magnitude of the longitudinal component at the jet axis. For these configurations, the jet is assumed to have an initially constant gas pressure ($`p_j`$) which is in equilibrium with the ambient gas pressure ($`\kappa =p_j/p_a=1`$).

The third adopted configuration is a purely toroidal magnetic field ($`B_\varphi `$) whose functional form is given by Fig. 1 and equation (5) of SH. Only the jet beam is initially magnetized in this case, and for the sake of clarification, we here rewrite SH’s equation for $`B_\varphi `$:

$$B_\varphi (r)=\{\begin{array}{cc}B_{\varphi ,m}\frac{r}{r_m}& 0rr_m\\ B_{\varphi ,m}\frac{R_jr}{R_jr_m}& r_mrR_j\\ 0& R_j<r\end{array}$$
$`(4)`$
where $`r_m`$ is a free parameters that determines the position where the magnetic field has its maximum intensity (we adopt here $`r_m0.9R_j`$). Lind et al. (1989), Frank et al. (1997, 1998) and OR, have also used a similar toroidal profile. In this case, in order to ensure initial magnetostatic equilibrium, the jet gas pressure has a radial profile with a maximum at the jet axis $`p_j(0)2.07p_j(R_j)`$, and $`p_j(R_j)=p_a`$ at the jet surface (see Figure 1 of SH).

The physical conditions at the jet inlet are parameterized by the following nondimensional numbers: $`\eta =n_j/n_a`$, the jet-to-ambient number density ratio; $`M_a=v_j/c_a`$, the average jet velocity to the ambient sound velocity ratio; $`\kappa =p_j/p_a`$, the jet-to-ambient thermal pressure ratio; $`\beta =8\pi p_{th}/B^2`$, the thermal-to-magnetic pressure ratio. As we will also study jets in which the magnetic field is a function of the radial distance $`r`$, we can define a mean magnetossonic Mach number as being $`M_{ms}=M(r)r𝑑r/r𝑑r`$, where $`M(r)=v_j/[c+v_A(r)]^{1/2}`$ (e.g., SH), and $`v_A`$ is the Alfven speed (see Table 1).

### 2.2  Magnetic Field Reversals

We should make some remarks on the late evolution of radiatively cooling magnetized jets containing longitudinal components. In Paper II, we have shown examples with the development of magnetic field reversals at the contact discontinuity between the jet and the cocoon with intensities up to 5 times their initial magnitude. As stressed in Paper II, field reversals of the longitudinal component occur in both purely longitudinal and helical magnetic field configurations, because the field lines are amplified by compression in the nonparallel shocks at the jet head, and are enforced to flow backward with the shocked plasma into the cocoon. In this process, the lines are reoriented and sometimes have their polarization reversed. The increasing strength of the reversed fields due to shear at the contact discontinuity could lead to the development of strong pinching regions in the late evolution of some flows that ultimately could cause jet disruption (see Fig. 16 of Paper II). Although shear and compression are expected to enhance B, we have speculated in Paper II that very large amplifications of the reversed fields could be due either to inappropriate computation of field dissipation under our ideal-MHD treatment, or to numerical effects.

In order to test the first of these hypotheses, we have carried out numerical simulations under a non-ideal MHD approximation, using different values for the magnetic resistivity. Although a complete analysis of the results of these tests is out of the scope of the present paper, we can briefly comment their consequences to the present study. We have found that a non-null magnetic resistivity ($`\eta _M`$) is able only to postpone the effects above at the jet/cocoon interface \[tested for a broad range of values of $`\eta _M`$ from the thermal value (which is $`T^{3/2}`$; e.g., Spitzer 1956), to values many orders of magnitude higher\].

These findings have suggested a possible numerical origin for the large field reversals found in the earlier calculations. Performing then, extremely high resolution SPMHD calculations, using up to 106 SPH particles and the same conditions as in Fig. 16 of Paper II, we could eliminate the anomalous amplification of the reversed fields and the disruption of the beam, while the other features that appeared in the original lower resolution model were left unaffected both at the head and along the jet beam. The anomalous amplification of the fields was, therefore, being caused essentially by a poorer interpolation of the physical quantities at the contact discontinuity in the lower resolution models. As the extremely high resolution simulations expend too much computational time and memory without significantly improving the overall results, we have instead, presently introduced in our code a numerical $`cleaning`$ $`switch`$, in order to eliminate the anomalous magnetic fields at the jet/cocoon interface, in calculations involving logitudinal magnetic fields. The cleaning switch allows the reduction of the strength of the magnetic fields at the jet/cocoon interface by an amount of about 1% whenever the reversed $`𝐁`$-field component exceeds a critical value that was chosen to be twice as large as the maximum value of the field in the head (where $`𝐁`$ regularly reaches its maximum physical amplification behind the strong shocks). In average, less than 0.01% of the SPH particles have their magnetic field vectors modified by such switch, which is applied in only less than 1% of the total number of time steps. The tests have indicated that the system is stabilized against anomalous field amplification and jet disruption with the employment of this switch and, more important, the results of these models are comparable with those obtained in the extremely high resolution calculations.

## 3 The simulations

In this section, we present the results of our 3-D, MHD numerical simulations of radiatively cooling pulsed jets. The physical parameters of our models are summarized in Table 1. All the parameters are suitable for the study of protostellar jets, namely: $`R_j2.5\times 10^{15}`$ cm, $`M_a15`$, or $`v_j15\times c_a250`$ km s-1; $`\eta =n_j/n_a=5`$; $`n_a=200`$ cm-1 (e.g., Reipurth & Raga 1999). We have adopted in most of the simulated models, a maximum intensity for the magnetic field as given by the equipartition condition: $`\beta =8\pi p_{th}/B^21`$, which implies a maximum intensity of $``$ 80 $`\mu `$G. In some cases, we have adopted $`\beta 0.1`$, which implies a maximum $`B260\mu `$G. The MHD models are all compared with the purely HD counterparts for which $`\beta =\mathrm{}`$.

### 3.1 Effects of different magnetic field configurations

Figure 1 displays the midiplane density contours (left) and the velocity field distributions (right) for four supermagnetosonic, radiatively cooling, pulsed jets with initial $`\beta 1`$ after they have propagated over a distance $`35R_j`$, at $`t/t_d3`$. The top jet is purely hydrodynamical (model HD1 in Table 1); the second jet (from top to bottom) has an initial constant longitudinal magnetic field (model ML1 in Table 1); the third jet has an initial helical magnetic field (MH1), and the bottom jet has an initial toroidal magnetic field (MT1). The early evolution of these jets was discussed in Paper III.

At the time depicted in Figure 1, the leading working surface at the jet head is followed by 4 other features. Each new pulse at the left forms within a time of about $`0.2t_d`$. This value agrees with that predicted by linear analysis of pulsed jets with sinusoidal injection profile (Raga et al. 1990; Raga & Cantó 1998). Like the leading working surface, each internal feature consists of a double-shock structure, an upstream reverse shock that decelerates the high velocity material entering the pulse, and a downstream forward shock sweeping up the low velocity material ahead of the pulse. Each of these internal working surfaces (IWS) or knots, widens and broadens as it propagates downstream and squeezes shocked material sideways into the cocoon. Since $`P0.5t_d`$, there has been time for 5 IWSs to form through steepenning of the input sinusoidal profile. The first of these IWSs has already overtaken the leading working surface and merged with it. With the impact, it is disrupted and its debris are partially deposited into the cocoon thus providing a complex fragmented structure at the head as required by the observations (as is the case, for example, of HH 1/2; e.g., Jeff Hester, Stapelfeldt & Scowen 1998).

As previously reported (see Paper III, and de Gouveia Dal Pino & Cerqueira 2001), Fig. 1 suggests that the overall morphology of the 3-D pulsed jet is not very much affected by the presence of the different magnetic field configurations in equipartition with the gas. Compared with the purely HD jet (top panel), the introduction of the distinct $`𝐁`$-profiles tends, instead, to alter essentially the detailed structure and the emission properties (see below) behind the shocks at both the head and internal knots, particularly when helical or toroidal fields are present. For example, like in steady-state calculations (see Paper I and II), while the MHD jet with initial longitudinal $`𝐁`$-field (second panel) exhibits a fragmented dense shell at the head which is very similar to that in the pure HD jet, the MHD jets with helical and toroidal fields present a more collimated head structure due to the action of the tension forces associated with the toroidal component of the $`𝐁`$-field (third and fourth panels; see also Paper III).

Several physical quantities along the beam axis of the jets of Fig. 1 with initial longitudinal (ML1), and toroidal (MT1) fields are displayed in Figures 2a and 2b, respectively. As previously reported for the helical jet of Fig. 1 (see Paper III), we see (top panel, Figs. 2a and b) that the initial sinusoidal velocity profile impressed on the flow at the inlet steepens into the familiar sawtooth pattern (e.g., Raga & Kofman 1992) as the faster material catches up with the slower, upstream gas in each pulse. The longitudinal field ($`B_x`$) depicted in Fig. 2a (bottom panel) is anti-correlated with the sharp density peaks within each IWS (second panel), where it attains a minimum value, and amplifies between them. On the other hand, the toroidal magnetic field ($`B_\varphi `$) depicted in Fig. 2b (bottom panel) sharpens within the knots and rarefies between them. These results are in agreement with those found in Paper III for the jet with initial helical geometry (see also Gardiner & Frank 2000).

Figures 3a and b depict the same HD jet (HD1, Fig. 3a), and the MHD jet with helical field (MH1, Fig. 3b) of Figure 1, but both more evolved, i.e., after they have propagated over a distance $`60R_j`$ ($`t/t_d5`$). We note in both the development of a gentle wandering or wiggling along the jet axis which is much more pronounced in the MHD jet. This is excited by the kink mode of the Kelvin-Helmholtz (K-H) instability. From the simulation we obtain a wavelength for the oscillation $`\lambda 16R_j`$, which is in good agreement with that derived for the resonant wavelength of the fundamental kink mode of the K-H instability from the linear theory (e.g., Hardee, Clarke & Rosen 1997; Hardee & Stone 1997, Hardee & Norman 1988). This effect could provide a potential explanation for the wiggling features often detected in HH jets (e.g., HH 111 and HH 46/47; Reipurth et al. 1997; Heathcote et al. 1996).

Figure 4 displays the temporal evolution of the density of the leading working surface (left panels) and of the second IWS (right panels) for the four jets of Figure 1: HD1 (top), ML1 (second panel), MH1 (third panel), and MT1 (bottom). As reported in previous 3-D calculations of steady-state jets (e.g., GB93; Paper I and II), the density of the leading working surface in the HD1 model (top panel, left) shows a temporal variability that is caused by global thermal instabilities (GTI) of the strong, high velocity radiative shocks (e.g., Gaetz, Edgar & Chevalier 1988; GB93; Paper II). Consistently with GTI theory (GB93), the density of the shell oscillates with a period which is of the order of the cooling time behind the shocks, $`\tau 2.2t_{cool}`$ (where $`t_{cool}0.77t_d`$ for the HD1 model as obtained from the initial conditions, or $`\tau 1.7t_d`$). The introduction of the different magnetic field configurations alters both the amplitude and the pattern of the density oscillations, particularly in the presence of the toroidal field (bottom panel, left), where the delay in the oscillations can be interpreted as being due to the increase of the cooling distance behind the shock caused by the amplification of the toroidal component (e.g., Paper I and II). In the right panels of Figure 4, we see that the density of the IWSs also varies in time due to the global thermal instabilities, with a period $`\tau 2t_d`$, and the oscillation pattern is also more affected in the presence of the toroidal field. In all cases, the oscillation peak decays with time as a consequence of both the sideways deposition of the shocked material into the cocoon and the damping of the GTI (GB93). \[This density decay is less obvious in the leading working surface (left panels) because of the continuous accumulation of material in the head due to the merging of the IWSs.\]

### 3.2 Effects of different pulsation periods

Figure 5 displays the midiplane density contours (left) and the velocity field distributions (right) for four pulsed jets after they have propagated over a distance $`40R_j`$, at $`t/t_d=3`$. The top jet is purely hydrodynamic (model HD2 in Table 1); the second (from top to bottom) has an initial constant longitudinal magnetic field (model ML2 in Table 1); the third has an initial helical magnetic field (MH2), and the bottom has an initial toroidal magnetic field (MT2). The initial conditions are the same as in Figure 1, but here the period of the injected sinusoidal velocity profile is twice as large ($`P=1t_d`$; see Table 1). The four jets in Figure 5 are slightly advanced in the computational domain when compared to those of Figure 1. This is a natural consequence of the larger period $`P`$ of the injected velocity variability (see, e.g., de Gouveia Dal Pino 2001). As the average propagation velocity of the leading bow shock ($`v_{bs}`$) is larger than that of the jets of Figure 1, a smaller amount of shocked material is accumulated behind the jet shock or Mach disk (whose velocity is $`v_{Md}v_jv_{bs}`$). Therefore, the dense, cold shell that develops at the head due to the radiative cooling of the shocked material is thinner and less clumpy in the larger period jets.

As before, the overall jet morphology is not very much affected by the presence of the different magnetic field configurations, although the tension forces associated with the toroidal magnetic field component in the MHD jets with initially helical and toroidal geometries (third and fourth panels) increases the jet collimation with respect to the HD jet. Global thermal instabilities also occur in these models. As in Figure 4, Figure 6 shows that the leading working surface (top-left panel) and the IWSs (here represented by the second IWS, top-right panel) of the purely HD jet also undergo temporal density oscillations due to the GTI. Likewise, the shape and amplitude of the oscillations are affected by the introduction of the magnetic fields, particularly in the helical (third panels) and toroidal (fourth panels) cases.

The H$`\alpha `$ intensity is known to be strong behind the knots of the observed HH jets. As discussed in Paper III, we can estimate the H$`\alpha `$ emissivity ($`I_{\mathrm{H}\alpha }`$) behind the IWSs using the results of our simulations and the relation given by Raga & Cantó (1998), $`I_{\mathrm{H}\alpha }\rho _dv_s^{3.8}`$, where $`\rho _d`$ is the downstream pre-shock density, and $`v_s`$ is the shock speed. In Paper III (see also de Gouveia Dal Pino & Cerqueira 2001), we have found for the MHD jets of Fig. 1 that the H$`\alpha `$ intensity behind the internal knots is amplified up to factors $``$ 3, and 4 in the presence of the helical and toroidal field geometries, respectively, relative to the HD jet. The same evaluation perfomed for the jets of longer period of Figure 5, results in amplifications of the H$`\alpha `$ intensity up to factors $``$ 5, and 4, for the MHD jets with helical and toroidal fields, respectively, with respect to the H$`\alpha `$ intensity in the HD jet, while for the MHD jet with initial longitudinal field, the H$`\alpha `$ intensity is essentially the same as that of purely HD jet (see Figure 7).

### 3.3 Effects of different values of $`\beta `$

Figures 8 and 9 compare the early evolution of jets with different magnetic field strengths. Figure 8 displays two MHD jets with initial longitudinal magnetic field geometry: a $`\beta 1`$ (ML1, middle panel), and a $`\beta 0.1`$ jet (ML1b, bottom panel), which are compared with a HD jet (HD1, $`\beta =\mathrm{}`$, top panel). Likewise, Figure 9 displays two MHD jets with initial helical magnetic field geometries: a $`\beta 1`$ (MH1, middle panel), and a $`\beta 0.1`$ jet (MH1b, bottom), which are also compared with the HD jet (HD1, top panel). The remaining initial conditions in these figures are the same as in Figure 1. In the $`\beta 0.1`$ models, the maximum initial magnetic field is $`260\mu `$G. For both magnetic field configurations, we note a considerable increase in the jet collimation with larger $`𝐁`$-field ($`\beta 0.1`$) that is obviously caused by the action of the much stronger confining magnetic tension and pressure forces. Nonetheless, the global characteristics of the flow found in the HD, and $`\beta 1`$ jets are retained in the $`\beta 0.1`$ flows, particularly in the longitudinal field configuration (Figure 8). For example, the clumpy structures that develop at the head of the HD, and $`\beta 1`$ jets as a consequence of the combined effects of the radiative cooling and Rayleigh-Taylor (R-T) instabilities (see Paper I and II), are still present in the $`\beta 0.1`$ jet. This result is somewhat in contradiction with the recent 2-D results of OR, who found that the growth of the R-T instability and, therefore, the development of fragmentation at the jet head, could be inhibited in the presence of longitudinal fields even at equipartition conditions. This could be due to the differences in the assumed input conditions of their model with respect to the present ones. We note, however, that according to R-T instability theory an enhancement of the instability is expected in the presence of longitudinal $`𝐁`$-field (e.g., Jun, Norman & Stone 1995).

Figure 10 depicts the velocity (top panel), density (second panel), toroidal (third panel), and longitudinal (bottom panel) magnetic field components along the beam axis 222The toroidal magnetic field component, $`B_\varphi `$, has been actually taken close to the jet surface other than along the jet’s central axis where its value is too small (see eq. 3b). for the $`\beta 1`$ jet of Figure 9 (model MH1b in Table 1), at $`t/t_d=2`$. As before, the toroidal component of the magnetic field (third panel) sharpens within the knots and rarefies between them, while the longitudinal component (forth panel) is stronger between the knots. For the $`\beta 1`$ models with helical and toroidal configurations (see Figs. 2b, and 9), we find that the toroidal magnetic field components behind the IWSs are amplified by factors of $``$ 1.5 with respect to the initial values. For the $`\beta 0.1`$ model (Figs. 9 and 10), we find only a slightly larger amplification of $`B_\varphi `$ behind the IWSs $`2`$.

As in previous 2-D calculations (SH, OR), we also find that the density behind the IWSs increases with increasing $`𝐁`$-strength (or decreasing $`\beta `$), particularly in the presence of toroidal field components (Figure 9). This density enhancement is caused again by the confining tension forces associated with the toroidal fields which are larger for lower $`\beta `$-jets. However, in our 3-D calculations, we find smaller density amplifications. For example, for the $`\beta 1`$, and 0.1 jets with helical fields (Figure 9), we find a maximum density amplification behind the IWSs of a factor $`1.4`$, and $`4`$, respectively, relative to the HD jet, and for the $`\beta 1`$ jet with purely toroidal field (Fig. 1), we find a maximum density amplification of a factor $`1.6`$. On the other hand, SH have found larger density amplifications by factors $`240`$ in their 2-D simulations of jets with purely toroidal fields and $`\beta 1`$, and $`0.25`$, respectively. In the 2-D simulations of OR, the magnetic compression of the IWSs by the toroidal fields is more drastic, causing the disruption of them, even in $`\beta 1`$ jets. Although these differences can be partially attributed to differences in the input conditions of each model (see §1), the results above suggest that the effects of the magnetic fields upon the jet structure tend to be smoothed out in 3-D flows (see discussion below).

## 4 Discussion and Conclusions

We have presented the results of fully three-dimensional (3-D), magnetohydrodynamic (MHD) numerical simulations of radiatively cooling, overdense, pulsed jets, using the Smoothed Particle Magnetohydrodynamic (SPMHD) technique. We have explored the role of magnetic fields in pulsed jets considering different magnetic field strengths ($`\beta 0.1,1`$, and $`\mathrm{}`$) and configurations: longitudinal (aligned with the jet axis), helical, and purely toroidal magnetic fields. Our results could be summarized as follows:

1 - Magnetic fields with intensities of a few 10$`\mu `$G, which are of the order of those inferred for Herbig-Haro (HH) jets (e.g., Morse et al. 1992, 1994), and which are in close equipartition of energy density with the thermal gas are able to improve collimation of the jet with respect to a non-magnetic jet. Such collimation effect is, as expected, more pronounced for magnetized jets with initial helical and toroidal fields.

2 - Although of relevance for the formation of structures in the shell at the jet head, the combined effect of non-uniform cooling and the Rayleigh-Taylor (R-T) instability (e.g., GB93; Paper I), which is particularly evident in the case of steady-state jets (see Paper I), becomes less important in pulsed jets. In these cases, the continuous impact of the internal working surfaces (IWSs) with the jet head can also cause shell fragmentation. While the presence of toroidal magnetic field components can modify the morphology close to the jet head inhibiting its R-T fragmentation in the early evolution of the jet, as previously reported (paper I, III), the later impact of the pulsed-induced internal knots causes the appearance of a clumpy, complex morphology at the jet head (as required by the observations of Herbig-Haro jets) even in MHD models with helical or toroidal configurations.

3 - The overall morphology of the jet does not seem to be very much affected by the presence of the different magnetic field configurations, or different magnetic field strengths (see also paper III). Nonetheless, the role of the magnetic fields in the determination of the detailed structure and the emission properties behind the shocks in the head and the IWSs cannot be ruled out. The 3-D simulations show that the density of the shell and IWSs undergoes oscillations with time (with periods of the order of the cooling time behind the shocks) that are caused by global thermal instabilities of the strong radiative shocks. The amplitude and shape of these oscillations are affected by the introduction of magnetic fields, particularly by toroidal configurations.

4 - The 3-D calculations also indicate that the the H$`\alpha `$ emissivity can be enhanced by factors up to $`4`$, and 5 in the MHD models with toroidal field components, and $`\beta `$ 1, and 0.1, respectively, when compared with the non-magnetic calculation. Using the same data from Hartigan, Raymond & Hartman (1987; see also Hartigan, Morse & Raymond 1994) that have allowed to obtain $`I_{\mathrm{H}\alpha }\rho _dv_s^{3.8}`$ for the H$`\alpha `$ intensity (Raga & Kofman 1992), we find that the intensity of the \[S II\] doublet is wealky dependent on the shock velocity ($`I_{[\mathrm{S}\mathrm{II}]}v_s^{0.2}`$). We can thus infer from our results and from these two relations, that the $`I_{[\mathrm{S}\mathrm{II}]}/I_{\mathrm{H}\alpha }`$ ratio (which is frequently used to determine the excitation level of the HH objects) can decrease by factors up to $`45`$ in the MHD models with helical or toroidal fields (for $`\beta 10.1`$, respectively), with respect to the HD jet. Thus if the HH jets are really imbedded in magnetic fields with strengths not far from the equipartition with the gas, the evaluations of the physical quantities (such as, ionization fraction, mass loss ratio, etc.), which are presently made without taking into account the magnetic fields, could be significantly modified. We note, however, that the estimates above of $`I_{[\mathrm{S}\mathrm{II}]}/I_{\mathrm{H}\alpha }`$ are very crude and preliminary. Future calculations coupling the (magneto)hydrodynamic calculation and the evolution of the chemical species will allow to obtain more precise estimates of the intensity ratios above. Also, it’s interesting to note that in the regime of IWS shock velocities, and preshock densities and $`𝐁`$-fields examined here ($`v_s30`$ km s-1, $`B<300\mu `$G, $`n_j1000`$ cm-3), previous “planar” shock calculations by Hartigan, Morse & Raymond (1994) have indicated that the $`I_{[\mathrm{S}\mathrm{II}]}/I_{\mathrm{H}\alpha }`$ ratio is less sensitive to the magnetic field intensity than in our calculations.

5 - The differences that arise in 2-D simulations with different magnetic field geometry and strengths seem to diminish in the 3-D flows. Like in 2-D models, the density behind the IWSs also increases with increasing $`𝐁`$-strength (or decreasing $`\beta `$) in the 3-D calculations, but at a lower rate. For jets with initial helical field and $`\beta =1`$, and 0.1, we have found maximum density amplifications behind the IWSs of factors $`1.4`$, and $`4`$ respectively, and for a $`\beta 1`$ jet with pure toroidal field, we have found a maximum density amplification of a factor $`1.6`$, relative to the HD jet. These amplification factors are substantially smaller than those found in the 2-D models (SH; OR).

6- The amplification of the toroidal field components behind the IWSs is also smaller in 3-D flows. We have found field amplification factors of $``$ 1.5 for MHD models with $`\beta 1`$ and helical or toroidal configurations, and amplification factors $`2`$ in the $`\beta 0.1`$ model with helical configuration, with respect to their initial magnetic fields. These values are much smaller than those detected, for instance, by SH, who found field amplifications up to $`7`$ and 60 in their $`\beta 1`$, and 0.25 models, respectively.

7 - We have detected no signatures for the formation of nose cones. Such elongated structures are often found to develop at the head of magnetized jets in 2-D simulations involving toroidal magnetic fields. They are formed as the shocked gas is confined between the Mach disk and the bow shock by the tension forces associated with the toroidal field (e.g., Lind et al. 1989; Köss, Müller & Hillebrandt 1990; OR; SH). Although our 3-D calculations are probably limited by lower numerical resolution, the typical extensions of the nose cones in the 2-D calculations are few times the jet radius (SH; OR), so that they should be apparent even in a low resolution calculation. Our results thus indicate that nose cones are unstable in 3-D. Perhaps more convincing is the fact that the jet shock (or Mach disk) is found to be very close to the bow shock in most of the HH jets, as is the case of the HH 47A which, following the interpretation of Morse et al. (1994), may have a transverse magnetic field of the order of $`150\mu `$G.

8 - We speculate that the lower density amplification and the absence of nose cones in the 3-D calculations could be, in part, explained by the fact that, as the magnetic forces are intrinsically three-dimensional, they cause part of the material to be deflected in a third direction, therefore, smoothing out the strong focusing of the shocked material that is otherwise detected in the 2-D calculations (see also de Gouveia Dal Pino & Cerqueira 2001).

9 - Our simulations involving helical magnetic field geometry reveal the development of a wiggling structure at the late evolution that is due to the excitation of the kink mode of the Kelvin-Helmholtz instability. This could potentially explain the gentle wandering often displayed by HH jets (e.g., HH 111, e.g., Reipurth et al. 1997; HH 46/47, e.g., Heathcote et al. 1996). The appearance of this mode whose development is not possible in 2-D axisymmetric models also reinforces the importance of fully 3-D calculations.

We are thankful to Alex Raga for his advice and suggestions. A.H.C. would like to thank the Brazilian agency FAPESP, that have supported this work under a Ph.D. Fellowship program (process 96/08182-3). E.M.G.D.P. acknowledges the Brazilian agencies FAPESP and CNPq for partial support. The simulations were performed on a cluster of Linux-based PC’s, whose purchase was made possible by FAPESP. Also, partial support from the PRONEX/FINEP (41.96.0908.00) project is acknowledged.


format: text

text:
# 1 Introduction

## 1 Introduction

Levy-stable laws are a rich class of probability distributions that allow skewness and fat tails and have many intriguing mathematical properties . They have been proposed as models for many types of physical and economic systems. There are several reasons for using Levy-stable laws to describe complex systems. First of all, in some cases there are solid theoretical reasons for expecting a non-Gaussian Levy-stable model, e.g. reflection off a rotating mirror yields a Cauchy distribution ($`\alpha =1`$), hitting times for a Brownian motion yield a Levy distribution ($`\alpha =0.5,\beta =1`$), the gravitational field of stars yields the Holtsmark distribution ($`\alpha =1.5`$) . The second reason is the Generalized Central Limit Theorem which states that the only possible non-trivial limit of normalized sums of independent identically distributed terms is Levy-stable . It is argued that some observed quantities are the sum of many small terms – asset prices, noise in communication systems, etc. – and hence a Levy-stable model should be used to describe such systems. The third argument for modeling with Levy-stable distributions is empirical: many large data sets exhibit fat tails (or heavy tails, as they are called in the mathematical literature) and skewness, for a review see . Such data sets are poorly described by a Gaussian model and usually can be quite well described by a Levy-stable distribution.

Recently, in a series of economic and econophysics articles the Levy-stability of returns has been rejected based on the log-log linear regression of the cumulative distribution function or the Hill estimator . In this paper we show that the cited estimation methods can give exponents well above the asymptotic limit for Levy-stable distributions with $`\alpha `$ close to 2, which results in overestimation of the tail exponent in finite samples. As a consequence, the reported value of the tail exponent $`\alpha `$ around 3 may suggest a Levy-stable distribution with $`\alpha 1.8`$.

## 2 Levy-stable distributions

Levy-stable laws were introduced by Paul Levy during his investigations of the behavior of sums of independent random variables in the early 1920’s . The lack of closed form formulas for probability density functions for all but three Levy-stable distributions (Gaussian, Cauchy and Levy), has been a major drawback to the use of Levy-stable distributions by practitioners. However, now there are reliable computer programs to compute Levy-stable densities, distribution functions and quantiles . With these programs, it is possible to use Levy-stable models in a variety of practical problems.

The Levy-stable distribution requires four parameters to describe: an index of stability (tail index, tail exponent or characteristic exponent) $`\alpha (0,2]`$, a skewness parameter $`\beta [1,1]`$, a scale parameter $`\sigma >0`$ and a location parameter $`\mu R`$. The tail exponent $`\alpha `$ determines the rate at which the tails of the distribution taper off, see Fig. 1. When $`\alpha =2`$, a Gaussian distribution results. When $`\alpha <2`$, the variance is infinite. When $`\alpha >1`$, the mean of the distribution exists and is equal to $`\mu `$. In general, the $`p`$–th moment of a Levy-stable random variable is finite if and only if $`p<\alpha `$. When the skewness parameter $`\beta `$ is positive, the distribution is skewed to the right. When it is negative, it is skewed to the left. When $`\beta =0`$, the distribution is symmetric about $`\mu `$. As $`\alpha `$ approaches 2, $`\beta `$ loses its effect and the distribution approaches the Gaussian distribution regardless of $`\beta `$. The last two parameters, $`\sigma `$ and $`\mu `$, are the usual scale and location parameters, i.e. $`\sigma `$ determines the width and $`\mu `$ the shift of the mode (the peak) of the distribution.

### 2.1 Characteristic function representation

Due to the lack of closed form formulas for densities, the Levy-stable distribution can be most conveniently described by its characteristic function $`\varphi (t)`$ – the inverse Fourier transform of the probability density function. However, there are multiple parameterizations for Levy-stable laws and much confusion has been caused by these different representations . The variety of formulas is caused by a combination of historical evolution and the numerous problems that have been analyzed using specialized forms of the Levy-stable distributions. The most popular parameterization of the characteristic function of $`XS_\alpha (\sigma ,\beta ,\mu )`$, i.e. a Levy-stable random variable with parameters $`\alpha `$, $`\sigma `$, $`\beta `$ and $`\mu `$, is given by :

$$\mathrm{log}\varphi (t)=\{\begin{array}{cc}\sigma ^\alpha |t|^\alpha \{1i\beta \mathrm{sign}(t)\mathrm{tan}\frac{\pi \alpha }{2}\}+i\mu t,\hfill & \alpha 1,\hfill \\ & \\ \sigma |t|\{1+i\beta \mathrm{sign}(t)\frac{2}{\pi }\mathrm{log}|t|\}+i\mu t,\hfill & \alpha =1.\hfill \end{array}$$
(1)
For numerical purposes, it is often useful  to use a different parameterization:

$$\mathrm{log}\varphi _0(t)=\{\begin{array}{cc}\sigma ^\alpha |t|^\alpha \{1+i\beta \mathrm{sign}(t)\mathrm{tan}\frac{\pi \alpha }{2}[(\sigma |t|)^{1\alpha }1]\}+i\mu _0t,\hfill & \alpha 1,\hfill \\ & \\ \sigma |t|\{1+i\beta \mathrm{sign}(t)\frac{2}{\pi }\mathrm{log}(\sigma |t|)\}+i\mu _0t,\hfill & \alpha =1.\hfill \end{array}$$
(2)
The $`S_\alpha ^0(\sigma ,\beta ,\mu _0)`$ parameterization is a variant of Zolotariev’s  (M)-parameterization, with the characteristic function and hence the density and the distribution function jointly continuous in all four parameters. In particular, percentiles and convergence to the power-law tail vary in a continuous way as $`\alpha `$ and $`\beta `$ vary. The location parameters of the two representations are related by $`\mu =\mu _0\beta \sigma \mathrm{tan}\frac{\pi \alpha }{2}`$ for $`\alpha 1`$ and $`\mu =\mu _0\beta \sigma \frac{2}{\pi }\mathrm{log}\sigma `$ for $`\alpha =1`$.

For simplicity, in Section 3 we will analyze only non-skewed ($`\beta =0`$) Levy-stable laws. This is not a very restrictive assumption, since most financial asset returns exhibit only slight skewness. For $`\beta =0`$ both representations are equivalent, however, in the general case the $`S^0`$ representation is preferred.

### 2.2 Simulation of Levy-stable variables

The complexity of the problem of simulating sequences of Levy-stable random variables results from the fact that there are no analytic expressions for the inverse $`F^1`$ of the cumulative distribution function. The first breakthrough was made by Kanter , who gave a direct method for simulating $`S_\alpha (1,1,0)`$ random variables, for $`\alpha <1`$. It turned out that this method could be easily adapted to the general case. Chambers, Mallows and Stuck  were the first to give the formulas.

The algorithm for constructing a random variable $`XS_\alpha (1,\beta ,0)`$, in representation (1), is the following :

* generate a random variable $`V`$ uniformly distributed on $`(\frac{\pi }{2},\frac{\pi }{2})`$ and an independent exponential random variable $`W`$ with mean 1;
* for $`\alpha 1`$ compute:

$$X=S_{\alpha ,\beta }\times \frac{\mathrm{sin}(\alpha (V+B_{\alpha ,\beta }))}{(\mathrm{cos}(V))^{1/\alpha }}\times \left(\frac{\mathrm{cos}(V\alpha (V+B_{\alpha ,\beta }))}{W}\right)^{(1\alpha )/\alpha },$$
(3)
where

$`B_{\alpha ,\beta }`$ $`=`$ $`{\displaystyle \frac{\mathrm{arctan}(\beta \mathrm{tan}\frac{\pi \alpha }{2})}{\alpha }},`$
$`S_{\alpha ,\beta }`$ $`=`$ $`\left[1+\beta ^2\mathrm{tan}^2{\displaystyle \frac{\pi \alpha }{2}}\right]^{1/(2\alpha )};`$
* for $`\alpha =1`$ compute:

$$X=\frac{2}{\pi }\left[\left(\frac{\pi }{2}+\beta V\right)\mathrm{tan}V\beta \mathrm{log}\left(\frac{\frac{\pi }{2}W\mathrm{cos}V}{\frac{\pi }{2}+\beta V}\right)\right].$$
(4)

Given the formulas for simulation of standard Levy-stable random variables, we can easily simulate a Levy-stable random variable for all admissible values of the parameters $`\alpha `$, $`\sigma `$, $`\beta `$ and $`\mu `$ using the following property: if $`XS_\alpha (1,\beta ,0)`$ then

$`Y=\{\begin{array}{cc}\sigma X+\mu ,\hfill & \alpha 1\text{,}\hfill \\ & \\ \sigma X+\frac{2}{\pi }\beta \sigma \mathrm{log}\sigma +\mu ,\hfill & \alpha =1\text{,}\hfill \end{array}`$
is $`S_\alpha (\sigma ,\beta ,\mu )`$. The presented method is regarded as the fastest and the best one known. It is widely used in many software packages, including S-plus and STABLE .

### 2.3 Tail behavior

Levy  has shown that when $`\alpha <2`$ the tails of Levy-stable distributions are asymptotically equivalent to a Pareto law. Namely, if $`XS_{\alpha <2}(1,\beta ,0)`$ then as $`x\mathrm{}`$:

$`P(X>x)=1F(x)`$ $``$ $`C_\alpha (1+\beta )x^\alpha ,`$
$`P(X<x)=F(x)`$ $``$ $`C_\alpha (1\beta )x^\alpha ,`$
where

$$C_\alpha =\left(2_0^{\mathrm{}}x^\alpha \mathrm{sin}xdx\right)^1=\frac{1}{\pi }\mathrm{\Gamma }(\alpha )\mathrm{sin}\frac{\pi \alpha }{2}.$$

The convergence to a power-law tail varies for different $`\alpha `$’s and, as can be seen in Fig. 2, is slower for larger values of the tail index. Moreover, the tails of Levy-stable distribution functions exhibit a crossover from a power decay with exponent $`\alpha >2`$ to the true tail with exponent $`\alpha `$. This phenomenon is more visible for large $`\alpha `$’s and will be investigated further in the next Section.

## 3 Estimation of the tail index

The problem of estimating the tail index (as well as other parameters) is in general severely hampered by the lack of known closed–form density functions for all but a few members of the Levy-stable family. Fortunately, there are numerical methods that have been found useful in practice.

Currently there exist three estimation procedures for estimating Levy-stable law parameters worth recalling. The first one is based on a numerical approximation of the Levy-stable likelihood function. The ML method, as it is called, was originally developed by DuMouchel  and recently optimized by Nolan . It is the slowest of the three but possesses well known asymptotic properties.

The second method uses a regression on the sample characteristic function. It is both fast and accurate (as long as we are dealing with a sample generated by a Levy-stable law). The regression procedure was developed in the early 1980’s by Koutrouvelis  and recently improved by Kogon and Williams .

The last, but not least, is the quantile method of McCulloch . It is the fastest of the three, because it is based on tabulated quantiles of Levy-stable laws. Yet it lacks the universality of the other two, since it is restricted to $`\alpha 0.6`$.

All presented methods work pretty well assuming that the sample under consideration is indeed Levy-stable. However, if the data comes from a different distribution, these procedures may mislead more then the Hill and direct tail estimation methods. And since there are no formal tests for assessing the Levy-stability of a data set we suggest to first apply the ”visual inspection” or non-parametric tests to see whether the empirical densities resemble those of Levy-stable laws.

### 3.1 Log-log linear regression

The simplest and most straightforward method of estimating the tail index is to plot the right tail of the (empirical) cumulative distribution function (i.e. $`1F(x)`$) on a double logarithmic paper, as in Figs. 2-7 (see also Refs. ). The slope of the linear regression for large values of $`x`$ yields the estimate of the tail index $`\alpha `$, through the relation $`\alpha =\text{slope}`$.

This method is very sensitive to the sample size and the choice of the number of observations used in the regression. Moreover, the slope around $`3`$ may indicate a non-Levy-stable power-law decay in the tails or the contrary – a Levy-stable distribution with $`\alpha 1.8`$. To illustrate this, we simulated (using Eq. (3)) samples of size $`N=10^4`$ and $`10^6`$ of standard symmetric ($`\beta =\mu =0`$, $`\sigma =1`$) Levy-stable distributed variables with $`\alpha =1.95`$ and $`1.8`$. Next, we plotted the right tails of the empirical distribution functions on a double logarithmic paper. For $`\alpha `$ close to 2 the true tail behavior (2.3) is observed only for very large (also for very small, i.e. the negative tail) observations, see Figs. 3 and 5, after a crossover from a temporary power-like decay. Moreover, the obtained estimates still have a slight positive bias: $`0.03`$ for $`\alpha =1.8`$ and $`0.12`$ for $`\alpha =1.95`$, which suggests that perhaps even larger samples than $`10^6`$ observations should be used.

To test the method for extremely large data sets we also simulated samples of size $`N=10^8`$ of standard symmetric Levy-stable distributed variables with $`\alpha =1.95`$ and $`1.8`$. As can be seen in Fig. 7 the improvement over one million records samples is not substantial. The tail index estimate (2.02) is closer to the true value of $`\alpha =1.95`$, but still outside the Levy-stable regime. Similarly, for $`\alpha =1.8`$ the method overestimated the tail exponent and returned $`\alpha =1.82`$.

If a typical size data set is used, i.e. $`10^4`$ observations or less as in Refs.  and some data sets in Refs. , the plot may be quite misleading. An empirical cumulative distribution function of a Levy-stable sample of size $`10^4`$ and $`\alpha >1.5`$ does not exhibit the true tail behavior ($`x^\alpha `$ decay), but a temporary power-like decay with the slope (more precisely: absolute value of the slope) significantly greater then 2: $`|\text{slope}|`$ $`=4.99`$ for $`\alpha =1.95`$ and 2.89 for $`\alpha =1.8`$; see Figs. 4 and 6 and compare with Figs. 3 and 5, respectively. Slight differences in the slope of the initial temporary power-like decay (4.91 in Fig. 3 compared to 4.99 in Fig. 4 and 2.95 in Fig. 5 compared to 2.89 in Fig. 6) are caused by inaccurate estimation for large observations in the smaller samples.

Figures 3-7 clearly illustrate that the true tail behavior of Levy-stable laws is visible only for extremely large data sets. In practice, this means that in order to estimate $`\alpha `$ we must use high-frequency asset returns and restrict ourselves to the most ”outlying” observations. Otherwise, inference of the tail index may be strongly misleading and rejection of the Levy-stable regime unfounded. In Figures 3, 5 and 7 we used only the upper 0.5% or less of the records to estimate the true tail exponent. In general, the choice of observations used in the regression is subjective and can yield large estimation errors.

### 3.2 Hill estimator

Hill  proposed a method for estimating the tail index that does not assume a parametric form for the entire distribution function, but focuses only on the tail behavior. The Hill estimator is used to estimate the (Pareto) tail index $`\alpha `$, when the upper tail  of the distribution is of the form: $`1F(x)=Cx^\alpha `$. If $`X_{(1)},X_{(2)},\mathrm{},X_{(N)}`$ is the order statistics, i.e. original sample ordered so that $`X_{(1)}X_{(2)}\mathrm{}X_{(N)}`$, drawn from a population with law $`F`$ then the Hill estimate of $`\alpha `$ based on the $`k`$ largest order statistics is:

$$\alpha _{Hill}(k)=\left(\frac{1}{k}\underset{n=1}{\overset{k}{}}\mathrm{log}\frac{X_{(n)}}{X_{(k+1)}}\right)^1.$$
(6)
Unfortunately, it is difficult to choose the right value of $`k`$. In practice, $`\alpha _{Hill}(k)`$ is plotted against $`k`$ and one looks for a region where the plot levels off to identify the correct order statistic . Algorithms for choosing ”optimal” $`k`$ have been proposed in the literature (see e.g. Refs. ), but usually give numbers in the vicinity of the plateau.

To illustrate the performance of the Hill estimator for Levy-stable laws, we simulated (using Eq. (3)) samples of size $`N=10^4,10^6`$ and $`10^8`$ of standard symmetric ($`\beta =\mu =0`$, $`\sigma =1`$) Levy-stable distributed variables with $`\alpha =1.95,1.8,1.5`$ and 1. Next, we plotted the Hill statistic $`\alpha _{Hill}(k)`$ vs. $`k`$ and compared estimated and true values of $`\alpha `$.

Figure 8 presents the results of the analysis for Levy-stable samples of size $`10^4`$. As has been reported in the literature , for $`\alpha 1.5`$ the estimation is within reasonable limits but as $`\alpha `$ approaches 2 the Hill estimate is well above the Levy-stable regime. In the ”extreme” case of $`\alpha =1.95`$ there is no single value of $`k`$ that gives the right value of the tail exponent! For the $`1.8`$-stable sample Hill estimates are close to the true value of $`\alpha `$ only for $`k<50`$, whereas the plateau and ”optimal” order statistic  are in the range $`k(200,400)`$ yielding $`\alpha _{Hill}2.85`$. Actually the results are very close to those of the log-log linear regression, see Figs. 4 and 6 and Table 1, since in both methods estimates are obtained from the largest order statistics.

Figures 9 and 10 present a detailed study of the Hill estimator for samples of size $`N=10^6`$ and $`10^8`$. Since for $`\alpha 1.5`$ the results were satisfactory even for much smaller samples, in this study we restricted ourselves to $`\alpha =1.95`$ and $`1.8`$. When looking at the left panels of Fig. 9, i.e. for $`k<10\%`$ of the sample size (as in Fig. 8), the flat regions suggest similar values of $`\alpha _{Hill}`$ as did the Hill plots in Fig. 8. However, when we enlarge the pictures and plot the Hill statistics only for $`k<0.1\%`$ of the sample size we can observe a much better fit. For the $`1.8`$-stable sample Hill estimates are very close to the true value of $`\alpha `$ for $`k(150,400)`$. But for the ”hopeless” case of $`\alpha =1.95`$ again there is no single value of $`k`$ that gives the right value of the tail exponent. Yet this time the estimates are very close to the Levy-stable regime. Like for smaller samples, the Hill estimates agree quite well with log-log regression estimates of the true tail, see Figs. 3 and 5. For the extreme case of $`10^8`$ observations the results are similar. The flat regions suggest almost the same values of $`\alpha _{Hill}`$ as did the Hill plots in Fig. 9. In fact the Hill plots for $`k<10\%`$ of the sample size look so much alike, that in Fig. 10 we plotted only the counterparts of the right plots of Fig. 9, i.e. the enlarged plots for $`k<0.1\%`$ of the sample size. We can see that for $`\alpha =1.8`$ a good estimate is obtained only for $`k=25000,\mathrm{},60000`$ (i.e. for $`k<0.06\%`$ of sample size), whereas for $`\alpha =1.95`$ the estimate is in the Levy-stable regime only for $`k<20000`$ (i.e. for $`k<0.02\%`$ of sample size).

## 4 Conclusions

In the previous Section we showed that widely used tail index estimates (log-log linear regression and Hill) can give exponents well above the asymptotic limit for Levy-stable distributions with $`\alpha `$ close to 2. As a result, tail indices are significantly overestimated in samples of typical size. Only very large data sets ($`10^6`$ observations or more) exhibit the true tail behavior and decay as $`x^\alpha `$. In practice, this means that in order to estimate $`\alpha `$ we must use high-frequency asset returns and analyze only the most outlying values. Otherwise, inference of the tail index may be strongly misleading and rejection of the Levy-stable regime unfounded.

Recently extremely large data sets have become available to the researchers. The largest ones studied in the literature so far (in the context of tail behavior) are (i) the 40 million data points record of 5 minute increments for 1000 U.S. companies during the two year period 1994-95  and (ii) the 1 million data points record of 1 minute increments for the Standard&Poor’s 500 index during the 13-year period 1984-96 . However, the former data set is not homogeneous. It is formed out of 1000 data sets of 40000 records for individual companies. Therefore strong correlations, which can conceal the true nature of asset returns, are present in the data. The estimated tail indices for individual companies (see Fig. 1(b) in Ref. ) were found to range from $`\alpha =1.5`$ to 5.5. As we have shown in the previous Section these values can be easily obtained for samples of Levy-stable distributed variables.

On the other hand, the second (ii) data set may be regarded as homogeneous . The tails of the distribution of 1 minute returns were reported to decay in a power-law fashion with an exponent $`\alpha =2.95`$ for positive and $`\alpha =2.75`$ for negative observations . In an earlier paper by the same authors  the tail exponents were reported to be 2.93 and 3.02, respectively, which shows that the log-log linear regression method is very sensitive to the choice of observations used. Moreover, in both papers the tails of the empirical cumulative distribution function curl upward for extreme returns, see Fig. 1(c) in  and Fig. 4(a) in . This suggests that for very large and very small observations the distribution could be fitted by a power-law with a much smaller exponent. For example, if in Fig. 4(a) of  we plot an approximate regression line for negative returns in the region $`20g100`$ we find that the power-law exponent is less then 2. This may indicate that the tail exponents reported in both papers did not refer to the true tail behavior, but to the initial power-like decay (see Figs. 3 and 5) and that the rejection of the Levy-stable regime was not fully justified.

In this paper we have shown that the reported estimated tail exponent around 3 may very well indicate a Levy-stable distribution with $`\alpha 1.8`$. This is consistent with earlier findings (for a review see ) where the returns of numerous financial assets (individual stocks, indices, FX rates, etc.) were reported to be Levy-stable distributed with $`\alpha (1.65,2]`$. However, nothing we have said demonstrates that asset returns are indeed Levy-stably distributed. Although the analyzed tail index estimates are not sufficient to reject Levy-stability, by no means can we rule out a leptokurtic non-Levy-stable distribution that has power-law tails with $`\alpha >2`$.

## 5 Acknowledgments

The paper greatly benefited from Dietrich Stauffer’s critical reading and helpful comments. The research of the author was partially supported by KBN Grant no. PBZ 16/P03/99.


text:
# Accurate Calculations of the Peierls Stress in Small Periodic Cells

## 1 Introduction

Long, low-mobility -screw dislocations control low temperature plastic behavior in bcc metals. Unlike their fcc counter parts, bcc metals violate Schmid behavior and have many active slip planes. The microscopic origins of such behavior is crucial in the understanding of their plastic behavior. Therefore detailed and accurate first principle calculations are invaluable.

The Peierls stress, the zero temperature limit of the critical stress for slip, has been calculated in a variety of ways using various types of boundary conditions . These calculations, in general, rely on empirical potentials to simulate various configurations and incorporate a relatively large number of atoms. In general three types of boundary conditions are used, fixed cylindrical boundary conditions , lattice greens functions  or periodic boundary conditions . For these approaches a large number of atoms is generally taken in order to minimize any artificial effects the boundary conditions may impose. Generally, this number is too large to be suited for ab initio calculations, even on the fastest and largest super-computers. Applying corrections to these effects is therefore important in order to have a reliable ab initio calculation.

The simplest approach is to use fixed cylindrical boundary conditions. However, care must be taken due to the mismatch of the boundary with the dislocation, particularly when the dislocation moves. This mismatch can be minimized by performing calculations in very large systems, effectively extrapolating to the limit of infinite cylinder size. The size of the system needed to extract accurate values can be greatly reduced by applying leading order corrections due to the mismatch in the boundary . However, even with these corrections, accurate calculations can generally be made only for a cylinder with a radius of $``$30$`\AA `$ or greater, which corresponds to $``$700 atoms when using periodic boundary conditions along the dislocation line. Such fixed cylindrical boundary conditions are generally ill-suited for an ab initio calculation due to both this large number of atoms and the artificial effects of the free surface at the boundary on the electrons. An alternative approach to account for boundary effects, and thereby reduce the size of the system needed for accurate calculation, is to use lattice greens function techniques . This technique, although quite elegant, still faces the issue of free surface effects when applied to first principles quantum mechanical calculations.

Periodic boundary conditions are in general the most natural and straight forward way to calculate the Peierls stress in a density functional theory calculation. In periodic boundary conditions one is forced to have a net zero Burgers vector per unit cell. This is generally accomplished by using either a dipole or quadrupole array of dislocations . Peierls stress calculations have been done in periodic boundary conditions using empirical potentials, however, to avoid dislocation-dislocation interactions, generally requires the use of unit cells  which would be impracticable for ab initio calculations. Wang et al. have introduced a promising method for use in small, periodic cells, but this approach requires the decomposition of the total energy into individual atomic contributions, which while feasible for inter-atomic potential approaches, is not well-defined in quantum mechanical calculations.

In this paper it is shown how an accurate Peierls stress can be obtained in very small periodic cells. It is shown how leading order effects, due to the array of closely packed dislocations, can be accurately accounted. It is further shown how an accurate Peierls stress can be obtained by simply applying a pure shear to the system and minimizing the energy. This is a tremendous advantage in ab initio calculations as it minimizes the number of atomic configurations which need be explored. As an application of these ideas, we focus specifically on the calculation of the Peierls stress in bcc tantalum for a -screw dislocation in which the maximum resolved shear stress is along a $`\{110\}`$. This is among the most relevant geometries in understanding the plastic response of this material.

The paper is organized as follows. Section 2 presents the results of Peierls stress calculations in very large cylindrical cells, as a reference against which we shall compare our new approach. Section 3 presents the results of calculations using periodic boundary conditions and underscores the problems that arise. Section 4 shows how to overcome these problems, and, finally, Section 5 presents several techniques for accelerating the calculation. Finally, Section 6 concludes the paper.

## 2 Reference Calculations

All calculations presented throughout this manuscript employ a quantum-based, embedded atom force field (qEAMFF) developed in . To provide the reference value of the Peierls stress for our calculations within periodic boundary conditions, we first perform calculations for isolated dislocations in a cylindrical geometry with fixed boundary conditions.

To form the dislocation, we proceed as follows. First, we define coordinates for our calculations as follows: the $`x`$-axis lies along the $`[1\overline{1}0]`$ direction, the $`y`$ lies along the $`[11\overline{2}]`$, and the $`z`$ lies along the $`[111]`$. A cylinder with radius $`R2`$ of Ta is then taken (Figure 1) which has two regions, a fixed region and a relaxation region. In the relaxation region ($`r<R1`$), where $`r`$ is the radial distance from the center, the atoms are allowed to relax in response to their interatomic forces. The atoms in the fixed region ($`R1<r<R2`$) are held fixed to their positions according to the solution of anisotropic elasticity theory. Finally, periodic boundary conditions (with period of one Burger’s vector) are employed along the direction of the dislocation. Thus, these calculations consider only straight, infinite dislocations. A -screw dislocation then is obtained by displacing all atoms in the cylinder according to the solution of anisotropic elasticity theory  and then optimizing the positions of the atoms in the relaxation region.

To calculate the Peierls stress, we then apply a strain to the system which ensures that the resulting stress has only one component, $`\sigma _{xz}`$, which generates a force on the dislocation line in the $`[11\overline{2}]`$ direction . We then optimize the locations of the atoms in the relaxation region subject to this external strain. Figure 2 shows the resulting dislocation structure in a large cyllinder ($`R1120\AA `$) for a series of different strains using a differential displacement (DD) map . In these maps, the circles represent columns of atoms viewed along the $`[111]`$ direction. The arrows indicates the change in relative displacement that neighboring atomic columns make, relative to the bulk, due to the presence of the dislocation. The lengths of the arrows are normalized so that a displacement of 1/3 of a Burgers vector corresponds to an arrow of full length. The first panel (a) shows the ground-state structure of the dislocation at zero stress. The center of the dislocation is located at a diamond encased in triad of arrows. Going around this triad makes three displacements of one-third of a Burgers vector for a net displacement of one full Burgers vector relative to the bulk. The ground state structure of the core is seen to break the symmetry of the lattice by extending outward along three $`\{110\}`$ planes. This “degenerate split core” is consistent with molecular dynamic results, using periodic boundary conditions, found in reference  when using the same interatomic potential.

As the strain is applied to the ground state structure (Figures 2b-d), the dislocation feels force along the $`(1\overline{1}0)`$-plane. As the applied strain obtains a critical value, the dislocation center moves one lattice spacing along the $`(1\overline{1}0)`$-plane (Figure 2b). As the strain is increased further, the dislocation then glides along the $`(2\overline{1}\overline{1})`$ (Figures 2 c,d) in two steps, first along $`(10\overline{1})`$ and then spontaneously along $`(1\overline{1}0)`$. There are consequentially two Peierls barriers that the dislocation must overcome, the first along $`(1\overline{1}0)`$ and the second along $`(10\overline{1})`$. This leaves some ambiguity for the definition of the Peierls stress. Different authors have used different definitions . For clarity we will consider the value of the $`\sigma _{xz}`$ stress for the first jump (Figure 2a-b) as the first Peierls stress ($`P1`$) and the value of the $`\sigma _{xz}`$ stress for the second jump (Figure 2b-c) as the second Peierls stress ($`P2`$).

To extract the limit of these critical stresses for an isolated dislocation in an infinite crystal, we have repeated the above calculations in cylinders of various sizes ranging from $`R1=30\AA `$ to $`R1=120\AA `$ and extracted the Peierls stress as a function of radius. Figure 3 summarizes our results. In the limit, we find values of 0.74 GPa and 0.91 GPa for $`P1`$ and $`P2`$, respectively.

As a consistency check that we indeed expect finite-size effects to be small for our largest cylinders, one can employ the method of Shenoy and Phillips . This methods estimates the unaccounted restoring stress that the boundary applies to a displaced dislocation to be

$$\frac{K_sb}{2\pi }A\frac{d}{R^2},$$

where $`K_s=(S_{11}/(S_{11}S_{44}S_{15}^2))^{1/2}`$ in terms of the modified elastic compliances $`S_{ij}`$ , $`b`$ is the Burgers vector, $`d`$ is the distance the dislocation has moved from the center of the cylinder, $`R`$ is the radius of the cylinder, and $`A`$ is a dimensionless constant that can be calculated through elasticity theory or computationally. In our case, $`A2`$, $`K_s=68`$GPa and $`b=2.9\AA `$. Finally, for $`P1`$ we have $`d.1\AA `$ and for $`P2`$ we have $`d2.7\AA `$. For our largest cylinder ($`R1=120\AA `$), we calculate a restoring stress of $`0.0004`$GPa and $`0.012`$GPa, respectively. Both values are quite small, well within the uncertainties in Figure 3.

## 3 Periodic Boundary Conditions

In periodic boundary conditions, the unit cell must contain a net zero Burgers vector. In practice this is generally accomplished through the use of either a dipolar or quadrupolar array of dislocations . Our calculations employ a quadrupole array, which has been shown to be the more appropriate choice for screw-dislocations . Figure 4 illustrates a cell containing 270 atoms.

As our first attempt to calculate the Peierls stress in periodic boundary conditions, we take a quadrupole array with lattice vectors fixed at the values corresponding to the perfect bulk material, that is the appropriate lattice vectors prior the insertion of dislocations. Below, we refer to this choice as “unrelaxed” lattice vectors, as they generally do not correspond to the lattice vectors of an unstrained quadrupole array.

To extract an estimate of the Peierls stress, we then follow a procedure analogous to that in Section 2, applying strain until the dislocations move. Here, rather than applying the strain to the fixed region ($`R1<r<R2`$), we strain the unrelaxed lattice vectors. Then, we then compute the critical stress from the strain through the elastic constant matrix,

$`\sigma `$ $`=`$ $`𝐂ϵ,`$ (1)
where $`\sigma `$ is a column vector of the stresses, $`ϵ`$ is the applied strain vector, and $`𝐂`$ is the elastic constant matrix, assumed to equal that of the bulk material. The advantage of this approach is that it requires exploration of a minimal number of configurations and bulk lattice and elastic constants are relatively easy to obtain from first principles calculations. The disadvantage of this approach is that its underlying assumptions cast doubt of the accuracy of the results for unit cells of modest size.

Figure 5 explores the convergence of this approach with increasing cell size for calculation of the first Peierls stress ($`P1`$). The two smallest cells contain 90 atoms ($`2.9\AA \times 24\AA \times 24\AA `$) and 270 atoms ($`2.9\AA \times 42\AA \times 42\AA `$) (or 45 and 135 atoms, if symmetry is exploited), and are the only cells convenient for detailed ab initio studies. The results in these cells, however, are extremely poor, with errors of $`200\%`$ and $`53\%`$, respectively. Thus, great care must be taken when working with such small cells and a method is needed to correct for these finite-size effects.

## 4 Corrections

The preceding calculations make two major assumptions: (1) that the quadrupole array of dislocations does not change the elastic constant matrix from that of bulk material, and (2) that the use of unrelaxed lattice vectors is appropriate. To explore the impact of these assumptions, we repeat the above calculations while using both the relaxed lattice vectors and elastic constants of the quadrupolar array.

Numerical calculations show that the symmetry of the elastic constant matrix of the quadrupolar array ($`𝐂^{}`$) is the same as that of the bulk, although the individual components may differ. In particular, the stress-strain relation takes the form

$`\sigma `$ $`=`$ $`𝐂^{}ϵ`$ (2)
$`\left[\begin{array}{c}\sigma _{xx}\\ \sigma _{yy}\\ \sigma _{zz}\\ \sigma _{yz}\\ \sigma _{xz}\\ \sigma _{xy}\end{array}\right]`$ $`=`$ $`\left[\begin{array}{cc}\left(4\times 4\right)& 0\\ 0& \left(2\times 2\right)\end{array}\right]\left[\begin{array}{c}ϵ_{xx}\\ ϵ_{yy}\\ ϵ_{zz}\\ ϵ_{yz}\\ ϵ_{xz}\\ ϵ_{xy}\end{array}\right]`$ (17)
In this equation, the lower $`2\times 2`$ sub-block couples the $`\sigma _{xz}`$ and $`\sigma _{xy}`$ stresses only to the corresponding strains, while the upper $`4\times 4`$ block couples the $`\sigma _{xx},\sigma _{yy},\sigma _{zz}`$ and $`\sigma _{yz}`$ stresses only to their corresponding strains. Note, therefore, that for the present Peierls stress calculations only the lower $`2\times 2`$ block is relevant.

Figure 6 shows the results (diamonds in the figure) of the extraction of the Peierls stress when using relaxed lattice vectors and the elastic constants for the quadrupole. The figure shows that the Peierls stress now converges much more quickly with cell size. (The figure does not include results for the smallest cell, which proved unstable to the relaxation of the lattice vectors.) For the smallest stable cell (135 atoms including symmetry), the error is reduced from $`53\%`$ to only $`18\%`$, indicating the possibility of extracting reasonable results from cells of size suited to ab initio calculations. Although results for reasonably sized cells are accurate, this approach, however, is not necessarily well suited for ab initio calculations because relaxation of the lattice vectors and the calculation of the elastic constant matrix requires the exploration of many new atomic configurations.

## 5 Accurate and Efficient Peierls Stress Calculation

Having found an accurate approach, we now explore how to minimize that calculations associated with the above corrections. We begin by considering the necessity of relaxing the lattice vectors and then consider computation of the relevant components of the elastic constant matrix $`𝐂^{}`$.

### 5.1 Benefit of fixed lattice vectors

Working with unrelaxed lattice vectors generates spurious strains in the unit cell. The question is whether these strains create spurious stresses which confound the extraction of the Peierls stress. The elastic force which any such stresses would generate on the dislocation take the form :

$$𝐅_L=𝐛\sigma \times \eta .$$
(18)
Here, $`𝐅_L`$ is the force per unit length on the dislocation, $`𝐛`$ is the Burgers vector, $`\sigma `$ is the stress tensor written as a matrix, and $`\eta `$ is the sense vector, the direction along which the dislocation runs. For a $`[111]`$-screw dislocation $`𝐛`$ and $`\eta `$ both lie along the $`z`$axis, and, therefore, only two components of stress can generate a force on the dislocation, $`\sigma _{xz}`$ and $`\sigma _{yz}`$, and confound the extraction of the Peierls stress. We now consider whether the strains associated with the unrelaxed lattice vectors can generate such stresses. We distinguish two types of strain, dilation and shear.

The presence of the quadrupolar array tends to dilate the unit cell in the (111) plane. Due to symmetry, this dilation tends to be uniform, as we have confirmed by direct numerical calculation on the unit cells. The above form of the elastic constant matrix $`𝐂^{}`$ prevents such dilation from generating any $`\sigma _{xz}`$ component. The second confounding component $`\sigma _{yz}`$ would also vanish were the $`4\times 4`$ subblock of the elastic constant matrix to have precisely the same symmetry as that of the bulk. Numerically, we find that that this is almost the case, and that the dilation contribution to $`\sigma _{yz}`$ ranges from $`0.02\%`$ (in our largest cell) to only $`4\%`$ (in the 135 atom cell) of the stress experienced in the unrelaxed cell. The diagonal components of strain $`ϵ_{ii}`$ do not significantly affect the extraction of the Peierls stress, and therefore need not be relaxed to extract meaningful Peierls stresses.

Turning to shear strains, the equal spacing of the dislocation array leads to zero shear within linear elasticity theory. From symmetry, however, core-core interaction (non-elastic) can generate a force in the $`[1\overline{1}0]`$-direction only, Figure 7, which results in one non-zero component of strain, $`ϵ_{yz}`$. From direct calculations it is shown that this is in fact the only non-zero component of shear strain. Relaxing this component of strain, therefore, creates a artificial material environment, different from what isolated dislocations would experience. To more quickly approach the limit of isolated dislocations, one should therefore not relax this component of the strain. This advantage of using of unrelaxed lattice vectors is well known and commonly exploited in the materials literature in the study of two-dimensional defects such as grain boundaries or surfaces. We therefore expect to be able to extract accurate Peierls stresses without the need for relaxing either off-diagonal or diagonal components of strain.

Figure 6 shows the results (circles in the figure) of extracting the Peierls stress when using unrelaxed lattice vectors, but while still computing the stresses with the appropriate elastic constant matrix $`𝐂^{}`$. The results converge even more quickly than those obtained by relaxing the lattice vectors, thus supporting our analysis. Therefore, extremely good results (already within $`2\%`$ in the 135 atom cell) can be obtained by not relaxing the lattice vectors. This not only produces results of far superior quality, it also reduces the computational effort.

### 5.2 Extraction of elastic constants

The results of Section 5.1, while quite impressive, still require calculation of the elastic constant matrix $`𝐂^{}`$ of the quadrupolar array. Significant savings can be gained with the simple realization that only the $`2\times 2`$ sub-block of the elastic constant matrix enters calculation of $`P1`$ and $`P2`$ through Eq. 17. Moreover, within linear elasticity theory, these components can be extracted equally well at the relaxed or unrelaxed lattice vectors, thereby again mitigating the need for relaxation of the lattice vectors.

Not all of the components of the $`2\times 2`$ sub-block,

$$𝐂_{}^{}{}_{2}{}^{}\left[\begin{array}{cc}C_{44}^{}& C_{16}^{}\\ C_{16}^{}& C_{55}^{}\end{array}\right],$$
(19)
need be computed independently.

To see this, consider the application of a pure, shear strain $`ϵ_{xz}`$ to the system. As a result of the symmetry of $`𝐂^{}`$, this generates only two components of stress, $`\sigma _{xz}`$ and $`\sigma _{xy}`$. From equation 18, however, $`\sigma _{xy}`$ does not generate any forces on the dislocation and thus, within linear elasticity, does not effect our calculation of the Peierls stress. The remaining stress, whose critical value is the Peierls stress, can be computed from just one component of $`𝐂_{}^{}{}_{2}{}^{}`$,

$$\sigma _{xz}=C_{44}^{}ϵ_{xz}.$$
(20)
Finally, we note that this pre-factor can be extracted without any additional calculation. Figure 8 shows the total energy of the quadrupole array plotted as a function of strain during the extraction of the critical Peierls stress. Prior to the first dislocation glide event at $`ϵ_{xz}0.015`$, the energy increases quadratically according to

$$\mathrm{\Delta }E=\frac{1}{2}C_{44}^{}ϵ_{xz}^2,$$
(21)
which contains precisely the same pre-factor as in Eq. 20.

Figure 6 shows the result (triangles in the figure) of the extraction of the Peierls stress from the application of a pure, shear strain $`ϵ_{xz}`$ to the unrelaxed lattice vecotrs and extracting the relevant elastic constant from the curvature of the energy prior to the glide events. The results are of nearly the same quality as working with unrelaxed lattice vectors and using the full elastic constant matrix $`𝐂^{}`$.

### 5.3 Second Peierls stress: $`P2`$

Figure 9 shows our preliminary results for the the extraction of the second Peierls stress $`P2`$. These results are complicated by the fact that after the first transition $`P1`$, the quadrupole array is now distorted and no longer a perfect quadrupole, thereby generating non-negligible dislocation-dislocation forces and further modifying the elastic constant matrix of the cell. As a possible correction to this effect, we are presently exploring moving the distorted cores back to their original quadrupole locations before further increasing the stress. Nonetheless, we find that significant improvements can be made by using the elastic constant matrix of a perfect quadrupole array and ignoring the dislocation-dislocation interactions. Finally, results for pure shear calculations are quite encouraging.

## 6 Conclusion

This paper present an accurate and effective way to calculate the $`\{110\}`$ Peierls stress in a $`[111]`$ screw dislocation in a bcc material. The results show that accurate results can be obtained even for the smallest cells while using unrelaxed lattice vectors and extracting the elastic constants directly from the calculations. The method most importantly appears to make ab initio calculations of Peierls stresses viable in periodic boundary conditions for the first time.