辛方法的校正公式

1996年Wisdom等提出了对辛方法进行校正的概念和实践，现在继续对辛校正进行详尽讨论和数值比较，尤其对哈密顿函数可分解为一个主要部分和多个次要部分的一般情形，用Lie级数推导任意阶的各种辛算法的一次和二次辛校正公式并对一些算法给出具体的辛校正公式。又以日、木、土三体问题为模型进行数值实验，结果表明一次辛校正能提高精度，改善数值稳定性。计算效率也比较高，因而值得推荐使用，辛方法通常用大步长数值积分，这时二次辛校正并没有显著提高结果的精度，却大大增加了计算时间，不应予以推荐。     

Propagation of solar generated disturbances through the solar wind critical points: One-dimensional analysis

We investigate the proper method for mathematically simulating the formation of an interplanetary disturbance (IPD) in the subsonic, sub-Alfvénic region near the solar surface within the constraints of one-dimensional hydrodynamic and magnetohydrodynamic (MHD) analyses. We then numerically simulate the subsequent propagation of the IPD through the solar wind critical points in the equatorial plane to the outer corona. We show that, if the IPD is initiated outside the critical points, it always contains both a forward and reverse shock (a shock pair). This result contrasts with observations indicating that shock pairs at 1 AU which can be associated with solar events are rare occurrences in the solar wind. On the other hand, IPDs initiated inside the critical points contain only a forward shock at the leading edge. When the magnetic field is included in the simulation and the IPD is originated inside the critical points, the IPD contains a forward shock at its leading edge followed by large-amplitude, nonlinear, MHD waves which are convected outward by the solar wind. Unlike shock pairs, MHD waves are often observed in the solar wind. Hence, we conclude that physically realistic studies of the propagation of IPD which are assumed to originate near the solar surface must (1) initiate the IPD inside the critical points and (2) include the magnetic field. Although this conclusion is based on a one-dimensional analysis, we speculate that it would be equally valid in multi-dimensions.     

Parallel/Vector Integration Methods for Dynamical Astronomy

This paper reviews three recent works on the numerical methods to integrate ordinary differential equations (ODE), which are specially designed for parallel, vector, and/or multi-processor-unit(PU) computers. The first is the Picard-Chebyshev method (Fukushima, 1997a). It obtains a global solution of ODE in the form of Chebyshev polynomial of large (> 1000) degree by applying the Picard iteration repeatedly. The iteration converges for smooth problems and/or perturbed dynamics. The method runs around 100-1000 times faster in the vector mode than in the scalar mode of a certain computer with vector processors (Fukushima, 1997b). The second is a parallelization of a symplectic integrator (Saha et al., 1997). It regards the implicit midpoint rules covering thousands of timesteps as large-scale nonlinear equations and solves them by the fixed-point iteration. The method is applicable to Hamiltonian systems and is expected to lead an acceleration factor of around 50 in parallel computers with more than 1000 PUs. The last is a parallelization of the extrapolation method (Ito and Fukushima, 1997). It performs trial integrations in parallel. Also the trial integrations are further accelerated by balancing computational load among PUs by the technique of folding. The method is all-purpose and achieves an acceleration factor of around 3.5 by using several PUs. Finally, we give a perspective on the parallelization of some implicit integrators which require multiple corrections in solving implicit formulas like the implicit Hermitian integrators (Makino and Aarseth, 1992), (Hut et al., 1995) or the implicit symmetric multistep methods (Fukushima, 1998), (Fukushima, 1999). This revised version was published online in July 2006 with corrections to the Cover Date.     

Symplectic Integrators: Rotations and Roundoff Errors

We investigate the numerical implementation of a symplectic integrator combined with a rotation (as in the case of an elongated rotating primary). We show that a straightforward implementation of the rotation as a matrix multiplication destroys the conservative property of the global integrator, due to roundoff errors. According to Blank et al. (1997), there exists a KAM-like theorem for twist maps, where the angle of rotation is a function of the radius. This theorem proves the existence of invariant tori which confine the orbit and prevent shifts in radius. We replace the rotation by a twist map or a combination of shears that display the same kind of behaviour and show that we are able not only to recover the conservative properties of the rotation, but also make it more efficient in term of computing time. Next we test the shear combination together with symplectic integrator of order 2, 4, and 6 on a Keplerian orbit. The resulting integrator is conservative down to the roundoff errors. No linear drift of the energy remains, only a divergence as the square root of the number of iterations is to be seen, as in a random walk. We finally test the three symplectic integrators on a real case problem of the orbit of a satellite around an elongated irregular fast rotating primary. We compare these integrators to the well-known general purpose, self-adaptative Bulirsch–Stoer integrator. The sixth order symplectic integrator is more accurate and faster than the Bulirsch–Stoer integrator. The second- and fourth- order integrators are faster, but of interest only when extreme speed is mandatory. This revised version was published online in August 2006 with corrections to the Cover Date.     

Basic Properties of Compressible MHD Turbulence: Implications for Molecular Clouds

Recent advances in understanding of the basic properties of compressible Magnetohydrodynamic (MHD) turbulence call for revisions of some of the generally accepted concepts. First, the MHD turbulence is not so messy as it is usually believed. In fact, the notion of strong nonlinear coupling of compressible and incompressible motions is not tenable. Alfven, slow and fast modes of MHD turbulence follow their own cascades and exhibit degrees of anisotropy consistent with theoretical expectations. Second, the fast decay of turbulence is not related to the compressibility of fluid. Rates of decay of compressible and incompressible motions are very similar. Third, the viscosity by neutrals does not suppress MHD turbulence in a partially ionized gas. Instead, MHD turbulence develops magnetic cascade at scales below the scale at which neutrals damp ordinary hydrodynamic motions. The implications of those changes of MHD turbulence paradigm for molecular clouds require further studies. Those studies can benefit from testing of theoretical predictions using new statistical techniques that utilize spectroscopic data. We briefly discuss advances in development of tools using which the statistics of turbulent velocity can be recovered from observations.     

Accretion and Ejections in Magnetized Disks

In this contribution I review the main challenges for theory and numerical simulation of accretion turbulence in disks. I then present briefly a solution we have elaborated in recent years to a part of these questions: we have found an MHD instability which occurs in the inner region of a disk in the configuration (poloidal field of the order of equipartition with the gas pressure) used for MHD jet models. This instability has the unique characteristic that it re-emits toward the corona a fraction of the energy and angular momentum it extracts from the disk. It is a good candidate to explain the low-frequency Quasi-Periodic Oscillation observed in X-ray binaries, and we believe that it might occur also in YSO.     

Reflection Properties of Gravito-MHD Waves in an Inhomogeneous Horizontal Magnetic Field

We derive the dispersion equation for gravito-magnetohydrodynamical (MHD) waves in an isothermal, gravitationally stratified plasma with a horizontal inhomogeneous magnetic field. Sound and Alfvén speeds are constant. Under these conditions, it is possible to derive analytically the equations for gravito-MHD waves. The high values of the viscous and magnetic Reynolds numbers in the solar atmosphere imply that the dissipative terms in the MHD equations are negligible, except in layers around the positions where the frequency of the MHD wave equals the local Alfvén or slow wave frequency. Outside these layers the MHD waves are accurately described by the equations of ideal MHD. We consider waves that propagate energy upward in the atmosphere. For the plane boundary, z=0, between two isothermal plasma regions with horizontal but different magnetic fields, we discuss the boundary conditions and derive the equations for the reflection and transmission coefficients. In the simpler case of a gravitationally stratified plasma without magnetic field, these coefficients describe the reflection and transmission properties of gravito-acoustic waves.     

The generation of MHD shock waves during the impulsive phase of the February 27, 1992 flare

In this paper a unique 2.3–4.2 GHz radio spectrum of the flare impulsive phase, showing fast positively drifting bursts superimposed on a slowly negatively drifting burst, is presented. Analyzing this radio spectrum it was found that the flare started somewhere near the transition region, where upward propagating MHD waves were generated during the whole impulsive phase. Moreover, it was found that behind a front of these ascending MHD waves the downward propagating electron beams, which bombarded dense layers of the solar atmosphere, were accelerated. It seems that, simultaneously with the increase of beam bombardment intensity, the intensity of MHD waves was increasing and thus the MHD shock wave generation and the electron beam acceleration and bombardment formed a self-consistently amplifying flare process. At higher coronal heights this process was followed by a type II radio burst, i.e. by the MHD flare shock. To verify this concept, the numerical modeling of the shock-wave generation and propagation in space from a flare site near the transition region up to 3 solar radii was made. Comparing the thermal and magnetic field disturbances, it was found that those of magnetic origin are more relevant in this case. Combining the results of interpretation and numerical simulation, a model of the February 27, 1992 flare is suggested and new aspects of this model are discussed.     

The February 1986 solar activity: A comparison of GIOTTO,VEGA-1, and IMP-8 solar wind measurements with MHD simulations

Large disturbances in the interplanetary medium were observed by several spacecraft during a period of enhanced solar activity in early February 1986. The locations of six solar flares and the spacecraft considered here encompassed more than 100° of heliolongitude. These flares during the minimum of cycle 21 set the stage for an extensive multi-spacecraft comparison performed with a two-dimensional, magnetohydrodynamic (MHD) numerical experiment. The plasma instruments on the European Space Agency (ESA)'s GIOTTO spacecraft, on its way to encounter Comet Halley in March 1986, made measurements of the solar wind for up to 8 hours per day during February. We compare solar wind measurements from the Johnstone Plasma Analyzer (JPA) experiment on GIOTTO with the MHD simulation of the interplanetary medium throughout these events. Using plasma data obtained by the IMP-8 satellite in addition, it appears that an extended period of high solar wind speed is required as well as the simulated flares to represent the interplanetary medium in this case. We also compare the plasma and magnetometer data from VEGA-1 with the MHD simulation. This comparison tends to support an interpretation that the major solar wind changes at both GIOTTO and VEGA-1 on 8 February, 1986 were due to a shock from a W05° solar flare on 6 February, 1986 (06:25 UT). The numerical experiment is considered, qualitatively, to resemble the observations at the former spacecraft, but it has less success at the latter one.     

Generation of Intermediate Drift Bursts by Magnetohydrodynamic Waves in the Solar Corona

The plasma mechanism of radio emission generation in an inhomogeneous medium is investigated. In the model under study, the electron beam with loss-cone distribution generates upper-hybrid waves that, in turn, are transformed into radio emission. It is shown that the influence of the plasma density inhomogeneity limits the plasma waves’ intensity considerably due to variation in their wave vector. The results are used to interpret the intermediate drift (IMD) bursts. A model is proposed in which these bursts are reflections of propagating small-scale (with amplitudes of about 1% and sizes of hundreds of kilometers) magnetohydrodynamic (MHD) disturbances of magnetic tubes. It is shown that this model allows us to explain the spectral parameters of the bursts in question. At present, the lack of precise and independent data about the magnetic field does not allow us to decide definitively between the existing models (whistler or MHD waves) of the IMD bursts; nevertheless, if the proposed model is correct, it can be used to determine the characteristics of the coronal MHD waves.     

Interfacing MHD Single Fluid and Kinetic Exospheric Solar Wind Models and Comparing Their Energetics

An exospheric kinetic solar wind model is interfaced with an observation-driven single-fluid magnetohydrodynamic (MHD) model. Initially, a photospheric magnetogram serves as observational input in the fluid approach to extrapolate the heliospheric magnetic field. Then semi-empirical coronal models are used for estimating the plasma characteristics up to a heliocentric distance of 0.1 AU. From there on, a full MHD model that computes the three-dimensional time-dependent evolution of the solar wind macroscopic variables up to the orbit of Earth is used. After interfacing the density and velocity at the inner MHD boundary, we compare our results with those of a kinetic exospheric solar wind model based on the assumption of Maxwell and Kappa velocity distribution functions for protons and electrons, respectively, as well as with in situ observations at 1 AU. This provides insight into more physically detailed processes, such as coronal heating and solar wind acceleration, which naturally arise from including suprathermal electrons in the model. We are interested in the profile of the solar wind speed and density at 1 AU, in characterizing the slow and fast source regions of the wind, and in comparing MHD with exospheric models in similar conditions. We calculate the energetics of both models from low to high heliocentric distances.     

An MHD gadget for cosmological simulations

Various radio observations have shown that the hot atmospheres of galaxy clusters are magnetized. However, our understanding of the origin of these magnetic fields, their implications on structure formation and their interplay with the dynamics of the cluster atmosphere, especially in the centres of galaxy clusters, is still very limited. In preparation for the upcoming new generation of radio telescopes (like Expanded Very Large Array, Low Wavelength Array, Low Frequency Array and Square Kilometer Array), a huge effort is being made to learn more about cosmological magnetic fields from the observational perspective. Here we present the implementation of magnetohydrodynamics (MHD) in the cosmological smoothed particle hydrodynamics (SPH) code gadget . We discuss the details of the implementation and various schemes to suppress numerical instabilities as well as regularization schemes, in the context of cosmological simulations. The performance of the SPH–MHD code is demonstrated in various one- and two-dimensional test problems, which we performed with a fully, three-dimensional set-up to test the code under realistic circumstances. Comparing solutions obtained using athena , we find excellent agreement with our SPH–MHD implementation. Finally, we apply our SPH–MHD implementation to galaxy cluster formation within a large, cosmological box. Performing a resolution study we demonstrate the robustness of the predicted shape of the magnetic field profiles in galaxy clusters, which is in good agreement with previous studies.     

Emulating the OPAL equation of state

The equation of state for the structure of the Sun and stars has to be precise to allow comparisons with observations, i.e., helioseismic inversions of thermodynamic quantities. Among the two of the most popular formalisms are (1) the OPAL equation of state developed at Livermore and (2) the Mihalas-Hummer-Däppen (MHD) equation of state. While OPAL has a solid theoretical foundation, and matches the observational data better, the MHD formalism is more intuitive, easy to realize, and has the possibility of adjustable parameters. Furthermore, it an open-source product in contrast to the proprietary OPAL. Recently a version of MHD has been obtained by including the so-called “Plank-Larkin partition function” and by adding scattering-state terms. The resulting formalism matches OPAL rather well. Here, we report on the next logical step, the implementation of this MHD upgrade into the simple and popular CEFF equation of state. Such an implementation will make it a flexible and convenient tool, allowing an approximative on-line implementation of OPAL in solar and stellar models.     

An unstable arch model of a solar flare

The theoretical consequences of assuming that a current flows along flaring arches consistent with a twist in the field lines of these arches are examined. It is found that a sequence of magneto-hydrodynamic (MHD) and resistive MHD instabilities driven by the assumed current (which we refer to as the toroidal current) can naturally explain most manifestations of a solar flare.The principal flare instability in the proposed model is the resistive kink (or tearing mode in arch geometry) which plays the role of thermalizing some of the field energy in the arch and generating X-configured neutral points needed for particle acceleration. The difference between thermal and nonthermal flares is elucidated and explained, in part, by amplitude-dependent instabilities, generally referred to as overlapping resonances. We show that the criteria for the generation of flare shocks strongly depend on the magnitude and gradient steepness of the toroidal current, which also are found to determine the volume and rate of energy release. The resulting model is in excellent agreement with present observations and has successfully predicted several flare phenomena.     

High precision symplectic integrators for the Solar System

Using a Newtonian model of the Solar System with all 8 planets, we perform extensive tests on various symplectic integrators of high orders, searching for the best splitting scheme for long term studies in the Solar System. These comparisons are made in Jacobi and heliocentric coordinates and the implementation of the algorithms is fully detailed for practical use. We conclude that high order integrators should be privileged, with a preference for the new $(10,6,4)$ method of Blanes et al. (2013).     

Jet simulations extending radially self-similar magnetohydrodynamics models

We perform a numerical simulation of magnetohydrodynamics (MHD) radially self-similar jets, whose prototype is the Blandford & Payne analytical example. The final steady state that is reached is valid close to the rotation axis and also at large distances above the disc where the classical analytical model fails to provide physically acceptable solutions. The outflow starts with a subslow magnetosonic speed, which subsequently crosses all relevant MHD critical points and corresponding magnetosonic separatrix surfaces. The characteristics are plotted together with the Mach cones and the superfast magnetosonic outflow satisfies MHD causality. The final solution remains close enough to the analytical one, which is thus shown to be topologically stable and robust for various boundary conditions.     

Compressible magnetohydrodynamic turbulence: mode coupling, scaling relations, anisotropy, viscosity-damped regime and astrophysical implications

We present numerical simulations and explore scalings and anisotropy of compressible magnetohydrodynamic (MHD) turbulence. Our study covers both gas-pressure-dominated (high β) and magnetic-pressure-dominated (low β) plasmas at different Mach numbers. In addition, we present results for super-Alfvénic turbulence and discuss in what way it is similar to sub-Alfvénic turbulence. We describe a technique of separating different magnetohydrodynamic modes (slow, fast and Alfvén) and apply it to our simulations. We show that, for both high- and low-β cases, Alfvén and slow modes reveal a Kolmogorov   k ^−5/3  spectrum and scale-dependent Goldreich–Sridhar anisotropy, while fast modes exhibit a   k ^−3/2  spectrum and isotropy. We discuss the statistics of density fluctuations arising from MHD turbulence in different regimes. Our findings entail numerous astrophysical implications ranging from cosmic ray propagation to gamma ray bursts and star formation. In particular, we show that the rapid decay of turbulence reported by earlier researchers is not related to compressibility and mode coupling in MHD turbulence. In addition, we show that magnetic field enhancements and density enhancements are marginally correlated. Addressing the density structure of partially ionized interstellar gas on astronomical-unit scales, we show that the viscosity-damped regime of MHD turbulence that we reported earlier for incompressible flows persists for compressible turbulence and therefore may provide an explanation for these mysterious structures.     

Steady and Time-Dependent MHD Modelling of Jets

A brief review is given of some results of our work on the construction of (I) steady and (II) time-dependent MHD models for nonrelativistic and relativistic astrophysical outflows and jets, analytically and numerically. The only available exact solutions for MHD outflows are those in separable coordinates, i.e., with the symmetry of radial or meridional self-similarity. Physically accepted solutions pass from the fast magnetosonic separatrix surface in order to satisfy MHD causality. An energetic criterion is outlined for selecting radially expanding winds from cylindrically expanding jets. Numerical simulations of magnetic self-collimation verify the conclusions of analytical steady solutions. We also propose a two-component model consisting of a wind outflow from a central object and a faster rotating outflow launched from a surrounding accretion disk which plays the role of the flow collimator. We also discuss the problem of shock formation during the magnetic collimation of wind-type outflows into jets.