* Your assessment is very important for improving the workof artificial intelligence, which forms the content of this project

# Download Numerical modeling of plasma structures and turbulence transport in

Survey

Document related concepts

Van Allen radiation belt wikipedia , lookup

Magnetohydrodynamic generator wikipedia , lookup

Standard solar model wikipedia , lookup

Solar observation wikipedia , lookup

Energetic neutral atom wikipedia , lookup

Geomagnetic storm wikipedia , lookup

Ionospheric dynamo region wikipedia , lookup

Solar phenomena wikipedia , lookup

Advanced Composition Explorer wikipedia , lookup

Transcript

Numerical modeling of plasma structures and turbulence transport in the solar wind DISSERTATION zur Erlangung des Grades Doktor der Naturwissenschaften an der Fakultät für Physik und Astronomie der Ruhr-Universität Bochum vorgelegt von Tobias Wiengarten aus Soest Bochum 2015 First advisor: PD Dr. Horst Fichtner Second advisor: Prof. Dr. Julia Tjus Abstract The thesis focuses on the modeling of the inner-heliospheric solar wind with its embedded magnetic field by numerically solving the equations of magnetohydrodynamics (MHD). The main numerical tool used throughout this work is the state-ofthe-art MHD code Cronos, which is extended from its basic setup to adequately describe the supersonic solar wind. First, solar rotation effects are implemented in both an inertial and a co-rotating reference frame, as well as in a conservative formulation for the latter. This setup is used to investigate the formation and evolution of co-rotating interaction regions (CIRs), which can be a dominant agent shaping the inner-heliospheric environment. For a direct comparison with spacecraft data, it is necessary to employ observation-based inner boundary conditions, for which a potential field model for the solar corona in conjunction with the empirically Wang-Sheeley-Arge model is employed. It is shown that this setup can reproduce measured solar wind conditions as being present during Carrington Rotations 2059-2061, when stable and long-lasting CIRs were observed. These results are utilized to investigate these structures’ influence on energetic particle propagation by means of a stochastic differential equation (SDE) approach (not part of this thesis). As the transport coefficients needed in these latter models also depend on the turbulent component of the heliospheric magnetic field, these are provided by solving a set of equations for the transport of low-frequency turbulence alongside the Reynolds-averaged MHD equations. This particular setup is applied to the propagation of a coronal mass ejection (CME) scenario, for which it is found that, on the one hand, turbulence does not act back strongly on the large-scale quantities, but, on the other hand, that CMEs are strong drivers of turbulence. Therefore, such a model is not required for studies focusing on large-scale CME quantities only, but it does allow for self-consistently computed transport coefficients for future energetic particle propagation studies. The implementation of all these extensions has been respectively validated by comparing with analytical models or previous numerical studies, which is an important part of working with numerical codes and, thus, of this thesis. Another line of work is the application of the solar findings to other stars, their stellar winds and resulting astrospheres. For hot stars, these huge cavities are found to modulate the Galactic Cosmic Ray flux, which is an effect that could be responsible for tiny-scale anisotropies in corresponding all-sky maps. The results presented in this work can be utilized for subsequent studies on numerous topics, such as energetic particle propagation, space-weather effects, or a coupling to outer-heliospheric models. iii Zusammenfassung Der Schwerpunkt dieser Arbeit liegt in der Modellierung des Sonnenwindes in der inneren Heliosphäre. Dazu werden die Gleichungen der Magnetohydrodynamik (MHD) mit Hilfe des Cronos-Codes numerisch gelöst, welcher einer entsprechenden Erweiterung zur Beschreibung des Sonnenwindes bedarf: Zunächst wird die Rotation der Sonne berücksichtigt, und zwar sowohl in einem inertialen als auch einem korotierenden Bezugssystem, wobei die letztere Beschreibung zudem konservativ formuliert wurde. Auf diese Weise wird die Untersuchung von korotierenden Wechselwirkungsregionen (CIRs) ermöglicht, welche einen großen Einfluss auf die Strukturierung der inneren Heliosphäre haben. Um anschließend detaillierte Vergleiche mit Messdaten von Raumsonden vornehmen zu können, ist es notwendig Beobachtungsdaten für die Formulierung der inneren Randbedingungen zu verwenden. Dies wird erreicht durch ein Potentialfeld-Modell für die solare Korona, dessen Topologie durch Anwendung des Wang-Sheeley-Arge Modells Aufschluss über die resultierenden Sonnenwindgeschwindigkeiten gibt, welche im anschließenden Vergleich mit Messdaten gute Übereinstimmung zeigen. Insbesondere Carrington Rotationen 2059-2061 werden untersucht, während derer stabile CIRs auftraten. Die Ergebnisse der MHD Simulationen können zur Untersuchung des Einflusses dieser Strukturen auf den Transport von energetischen Teilchen genutzt werden. Die diesen Transport beschreibenden Parameter hängen allerdings auch von der Turbulenz im Sonnenwind ab, sodass der Cronos-Code dahingehend erweitert wurde, die Gleichungen für den Transport von Turbulenz selbst-konsistent mit den MHD Gleichungen zu koppeln. Dieses Modell wird angewendet auf die Evolution eines koronalen Massenauswurfs (CME) mit dem Ergebniss, dass Turbulenz kaum Auswirkungen auf die großskalige Ausbreitung eines CME hat, aber ein CME die Turbulenzentwicklung maßgeblich antreibt. Somit stellen die Resultate eine Möglichkeit dar, selbst-konstitent berechnete Transportkoeffizienten für anschließende Studien zur Propagation energetischer Teilchen in CMEs zu nutzen. Die Implementierung all dieser Erweiterungen wird jeweils validiert durch detaillierte Vergleiche mit analytischen Modellen oder bestehenden numerischen Ergebnissen, was einen wichtigen Teil von Studien mittels numerischer Methoden, und somit dieser Arbeit, darstellt. Ein weiteres Arbeitsfeld ist die Anwendung der bezüglich unserer Heliosphäre gewonnen Erkenntnisse auf andere Sterne und deren Astrosphären. Es wird gezeigt, dass letztere bei heißen Sternen so groß sind, dass sie galaktische kosmische Strahlung modulieren und somit für kleinskalige Anisotropien in entsprechenden Messdaten verantwortlich sein könnten. Die Ergebnisse dieser Arbeit können vielseitig weiter verwendet werden, unter anderem für Studien zum Transport energetischer Teilchen, Weltraumwetter oder als Startpunkt für Simulationen der äußeren Heliosphäre. v Contents 1 Introduction 1.1 General overview . . . . . . . . . . . . 1.1.1 The Sun . . . . . . . . . . . . . 1.1.2 The solar wind . . . . . . . . . 1.1.3 The heliosphere . . . . . . . . . 1.1.4 Stellar winds and astrospheres . 1.2 Model formulation and numerical setup 1.2.1 Principal simplifications . . . . 1.2.2 The Cronos code . . . . . . . 1.3 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 1 6 12 15 17 17 18 22 2 MHD simulation of the inner-heliospheric magnetic field 27 2.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2 Wiengarten et al. (2013) – MHD simulation of the inner-heliospheric magnetic field Journal of Geophysical Research . . . . . . . . . . . . . . . . . . . . . 32 2.3 Errata and further development . . . . . . . . . . . . . . . . . . . . . 49 2.4 Addenda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.4.1 Weber-Davis model in the inertial and co-rotating reference frame . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.4.2 Divergence-free initialization for arbitrary input data . . . . . 54 3 Modeling background solar wind 57 3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.2 Wiengarten et al. (2014) – Cosmic Ray transport in heliospheric magnetic structures. I. Modeling background solar wind using the Cronos magnetohydrodynamic code The Astrophysical Journal . . . . . . . . . . . . . . . . . . . . . . . . 61 3.3 Further development and future directions . . . . . . . . . . . . . . . 78 4 Implementing turbulence transport 81 4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 vii CONTENTS 4.2 Wiengarten et al. (2015) – Implementing turbulence transport in the Cronos framework and application to the propagation of CMEs The Astrophysical Journal . . . . . . . . . . . . . . . . . . . . . . . . 85 4.3 Further development . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.4 Addendum - Conservative treatment of the co-rotating frame of reference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5 Extended turbulence transport 5.1 Model formulation . . . . . . 5.2 Model validation . . . . . . . 5.3 Model extensions . . . . . . . model 105 . . . . . . . . . . . . . . . . . . . . . . 105 . . . . . . . . . . . . . . . . . . . . . . 108 . . . . . . . . . . . . . . . . . . . . . . 111 6 Modeling of astrospheres 117 6.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 5.2 Scherer et al. (2015) – Cosmic rays in astrospheres Astronomy & Astrophysics . . . . . . . . . . . . . . . . . . . . . . . . 120 7 Summary and outlook 129 A Non-equidistantly spaced grids A.1 General overview . . . . . . . . . . . A.2 Application . . . . . . . . . . . . . . A.2.1 NESG for the radial direction A.2.2 NESG for the polar angle . . A.3 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 . 135 . 137 . 138 . 138 . 140 Bibliography 141 Acknowledgments 157 viii Chapter 1 Introduction This chapter starts with a general overview of the heliospheric environment (Section 1.1), before outlining principal simplifications for modeling purposes and introducing the main numerical code used in this work in Section 1.2. The further outline of the thesis – being based on four publications, and putting them into perspective regarding their significance for solar and heliospheric research – is given in Section 1.3. 1.1 General overview This section gives a general overview of the protagonists in this thesis: the Sun with its magnetic field (Section 1.1.1), emanating and shaping the turbulent solar wind (Section 1.1.2) that blows a bubble into the interstellar medium to form the heliosphere shielding the Earth from Galactic Cosmic Rays (Section 1.1.3). The hierarchy of the overview is to start from the center of the Sun and advance further and further outwards to the outer boundary of the heliosphere, while eventually the application of the solar findings to other – and thus even more distant – stars and their astrospheres (Section 1.1.4) is outlined. 1.1.1 The Sun Before embarking on describing the different distinguishing parts of the Sun, a summary of its large-scale features is given (see, e.g., Stix, 2004; Priest, 2014). The Sun is a main-sequence star of spectral type G2 V and is about 4.9 billion years of age. Although it is a rather ordinary star, Earth’s proximity with fairly intelligent lifeforms inhabiting it makes it possibly a unique one. Studying the Sun is of great importance as it has profound influence on our climate, and since it is the main driver for space weather (e.g., Baker, 1998) it is crucial to understand the Sun and the solar wind to circumvent dangers to space-borne and earthbound technology, on which mankind increasingly depends. Furthermore, the insights into stellar physics 1 CHAPTER 1. INTRODUCTION gained from studying our close-by star can be transferred to other, more distant stars. The Sun consists mainly of hydrogen (92%) and helium (8%) with low abundances of heavier elements (0.1%) so that the mean mass density ρ = 1.4 × 103 kg/m3 , while the total mass of M = 1.99 × 1030 kg constitutes 99,86% of the mass in the solar system. On a grand scale the Sun is in hydrostatic equilibrium, where the gravitational pull is balanced by the radiation pressure amounting to a luminosity of L = 3.86 × 1026 W. The radius of the Sun is defined as R = 6.955 × 108 m, so that the mean heliocentric distance of the Earth is one Astronomical Unit (1 AU = 1.5 × 1011 m ≈ 215R ). The outer layers of the Sun exhibit differential rotation with respective synodic (as seen from Earth) rotation periods ranging from 26.24 days at the equator up to about 36 days toward the poles (Snodgrass and Ulrich, 1990). The mass-loss rate due to the emitted solar wind is |Ṁ | = 109 kg/s. The solar interior The interior of the Sun and especially its inner core are inaccessible to direct observations due to the extremely high opacity, although the relatively new field of helioseismology (see, e.g., the review of Kosovichev, 2011) allows for indirect conclusions about the inner structure of the Sun by measuring oscillations on the Sun’s visible surface, the photosphere. However, as the Sun constantly emits energy in the form of radiation, the question as to how the energy is generated inside the Sun has puzzled mankind for a long time. It was not until the early 20th century – after establishing quantum mechanics – that it was recognized by Eddington that thermonuclear fusion held the answer, whereas earlier ideas like gravitational contraction, which put the Sun’s lifetime at a few million years only, could be discarded. The core’s extreme density (ρc ≈ 1.5 × 105 kg/m3 ) and temperature (Tc ≈ 1.5 × 107 K) are sufficient to drive fusion processes like the proton-proton chain or the CNO cycle, in which four protons are fused into one helium nucleus, and the mass difference between input and output is converted into energy according to Einstein’s E = mc2 . A byproduct of these processes are extremely weakly interacting neutrinos, whose detection on Earth (Davis et al., 1968) is another confirmation of the theory that the Sun’s energy originates from fusion. The transport of this energy within the solar interior is two-fold (see Figure 1.1). First, above the core from about 0.25R outwards, radiation is the dominant transport process: photons, initially with energies corresponding to gamma rays, are constantly absorbed and re-emitted, leading to a random walk inside the Sun that takes about two hundred thousand years (Mitalas and Sills, 1992) until they leave the Sun mainly as visible light. Second, above some 0.7R convective instability sets in due to a very large temperature gradient (see Schwarzschild, 1992). Consequently, blobs of hot plasma rise, expand and cool so that they eventually descend 2 1.1. GENERAL OVERVIEW Figure 1.1: Schematic illustration of the solar interior and the lower atmospheric layers (taken from Priest, 2014). again. The convection zone exhibits a number of complex motions such as differential rotation, meridional flow and 5-minute oscillations, while at the photosphere, convection is manifested in the so-called granulation. Granules occur at different scales (meso- and supergranulation), but a common feature are the bright centers containing upwelling hot plasma and dim edges of descending cool plasma (see left panel of Figure 1.2). The transition region between the radiative interior and the convection zone is called tachocline (Spiegel and Zahn, 1992), which is a region of strong shear that is thought to be largely responsible for the generation of the solar magnetic field by means of a dynamo process. Magnetic flux emerges through the photosphere on a variety of scales, which can be pictured as a carpet of magnetic loops closing at different heights (see right panel of Figure 1.2). The most prominent feature caused by the magnetic field in the photosphere are sunspots, which appear darker than their surroundings as the magnetic field inhibits convection so that cooled plasma is prevented from descending quickly (Figure 1.2 (left)). Sunspots have radii ranging from a few kilometers up to 15000 km, which make the largest ones visible to the 3 CHAPTER 1. INTRODUCTION Figure 1.2: Left: Granules and sunspots in the photosphere seen in the H-alpha wavelength on August 4, 2003 with the Swedish 1-m Solar Telescope (SST). [http://www.boston.com/bigpicture/2008/10/the sun.html]. Right: Solar magnetic loops captured by the Transition Region And Coronal Explorer. [http://sunearthday.gsfc.nasa.gov/2008/images/gal 023a art.jpg]. naked eye from the Earth under certain favoring conditions. A central dark umbra and a less dark, highly filamentary penumbra can be distinguished. Sunspots are found as pairs of opposite polarities within active regions located within ±30◦ of the equator, and their number is a measure for solar activity that roughly follows an 11-year cycle, during which the Sun exhibits many (solar maximum) to almost no (solar minimum) sunspots. The Sun’s lower atmosphere In contrast to the solar interior, the solar atmosphere is the part of the Sun that photons can directly escape from. Most of the solar radiation originates from the photosphere, which radiates like a black body at a temperature of about 6000K. Contrary to the expectation that the temperature should further decline with radial distance, there is actually an increase in temperature in the chromosphere with a sharp peak in the temperature gradient at the so-called transition region giving rise to coronal temperatures of several million degrees (Figure 1.3, with temperature still increasing beyond the heights shown here). The respectively labeled coronal heating problem remains unsolved, although the most likely candidates for the heating mechanisms are either damping and dissipation of waves originating from the turbulent motions in the convection zone (e.g., Carlsson et al., 2007), or magnetic reconnection (e.g., Parker, 1972; Cargill and Klimchuk, 2004). There is a vast number of magnetic structures in the solar atmosphere, which are often organized in networks. A detailed description would go beyond this brief overview (but see, e.g., Priest, 2014). The most common features in the chromosphere are spiculae of heated, accelerated and uprising, as well as cooled, descending plasma. The extension of the chromosphere is roughly 2.5 Mm, the transition region 4 1.1. GENERAL OVERVIEW Figure 1.3: Temperature and density profiles of the solar atmosphere (taken from Priest, 2014). is rather narrow (a few 100 km), and the corona’s outer boundary is usually defined by the Alfvénic critical surface at about rA ≈ 15 − 20R . At this surface the solar wind speed equals the Alfvén speed (vsw (rA ) = vA (rA )), so that no information can travel back to the corona. Two distinct topologies are typical for the coronal magnetic field, namely regions of closed loops on top of active regions, and open field lines in so-called coronal holes above unipolar regions. During solar minimum conditions the global magnetic field is approximately that of a dipole due to the strong unipolar fields of respectively opposing polarity in both polar regions, and large coronal holes occupy the higher latitudes, while closed loops are restricted to the equatorial regions (left panel of Figure 1.4). The influence of the solar wind (see Section 1.1.2)) on the magnetic field becomes more dominant with increasing radial distance and effects an opening of the field at mid-latitudes as well, usually at about 2.5R (Altschuler and Newkirk, 1969), and a radial magnetic field at all latitudes is present beyond the Alfvénic critical surface. During the solar cycle the large polar coronal holes become ever less pronounced towards solar maximum and eventually the predominant polarity of the northern and southern hemisphere changes sign. A magnetic- or Hale-cycle therefore takes about 22 years. In the more active phases around solar maximum the coronal field exhibits a more complex structure, where one can find coronal holes and closed loops at all latitudes (right panel of Figure 1.4). 5 CHAPTER 1. INTRODUCTION Figure 1.4: Corona near solar minimum (left) and maximum (right) as seen during solar eclipses (taken from Lang, 2006). 1.1.2 The solar wind Large-scale structure The history of the solar wind is a rather modern one, having found its way into science only in the second half of the 20th century, although it had been proposed already in the 19th century by Sun-observing pioneers such as Carrington and Birkeland that solar activity was responsible for geomagnetic effects such as the aurora. In the 1950s, Biermann (1951) proposed a corpuscular emission from the Sun to explain the gaseous straight tails of comets. Meanwhile, Chapman and Zirin (1957) argued that – due to the high coronal temperatures already known from X-ray observations – the coronal plasma should expand at least beyond Earth’s orbit. The latter authors, however, assumed a static configuration. By first demonstrating that such a static solution results in finite – and thus unrealistic – pressure at infinite distance, the modern basic idea of the solar wind as a dynamic, supersonic flow was postulated by Parker (1958). The driving force accelerating the solar wind in his isothermal and hydrodynamic model is the thermal pressure of the hot corona, resulting in terminal solar wind speeds in a range from about 400 to 800 km/s (depending on temperature, see left panel of Figure 1.5). The magnetic field is ”frozen” into (i.e., it is passively advected with) the solar wind, while the footpoints of the field lines remain anchored at the rotating Sun resulting in an Archimedean spiral pattern (see right panel of Figure 1.5). Parker’s idea met some resistance at first, but when the first measurements (e.g. with the Mariner 2 spacecraft, Snyder and Neugebauer, 1963) confirmed a supersonic flow it became broadly accepted and solar wind measurements made during the space age have confirmed the general picture. The solar wind consists mainly of protons and electrons, with low abundances of alpha particles. Typical values of solar wind 6 1.1. GENERAL OVERVIEW Figure 1.5: Parker’s solar wind model for radial velocity (left) and embedded magnetic field lines bent into Archimedean spirals by solar rotation (right) (taken from Parker, 1958). Speed Number density Proton temperature Magnetic field Average value 400 km/s 6.5 cm−3 50000 K 6 nT Table 1.1: Average observed solar wind properties at 1 AU in the ecliptic plane (Priest, 2014). properties at 1 AU in the ecliptic plane are summarized in Table 1.1. At the average speed of about 400 km/s, a solar wind parcel takes about 4-5 days to travel the distance from the Sun to the Earth. For a long time in-situ observations were restricted to the ecliptic plane. It was only after 1990 when the Ulysses spacecraft was brought into a high-latitude orbit after a swing-by maneuver at Jupiter that the three-dimensional picture of the solar wind topology was revealed. During Ulysses’ 19-year mission it covered almost two solar cycles, and the data taken by Ulysses are nicely summarized in plots such as the one presented in Figure 1.6, which show the following pattern: near activity minimum the solar wind is highly structured, with remarkably steady, fast, hot and tenuous streams (about 750 km/s) towards the poles, while the region within about ±20◦ of the solar equatorial plane is occupied by variable, slow, dense and cooler streams (about 400 km/s). This pattern was measured during the first (1992-1998; left panel of Figure 1.6) and third (2004-2009; not shown) orbits of the mission. In contrast, there is no discernible global structure during the maximum of solar activity, when fast as well as slow streams can be found at any given latitude. Naturally, this raises the question as to where the different types of solar wind originate from 7 CHAPTER 1. INTRODUCTION in the solar corona and how they are accelerated. As the plasma in the corona has to follow the magnetic field lines, only open field lines allow for an escape. The fast solar wind is clearly associated with the large coronal holes, but the origin of the slow solar wind is still unclear and under debate (see, e.g., Shergelashvili and Fichtner, 2012). Candidates are high-latitude coronal-hole/streamer boundaries for solar minimum, and small coronal holes and active regions for solar maximum. Being related to the coronal heating problem, the acceleration processes for the different types of solar wind are also still subject to ongoing research. One of the prevailing ideas is that the dissipation of Alfvén waves in the corona is the main driver to accelerate the solar wind, where their energy is mainly dumped below/above the Alfvénic critical surface for the slow/fast wind (Cranmer and van Ballegooijen, 2005; Cranmer et al., 2007). Figure 1.6: Solar wind velocity profiles near solar minimum (left) and maximum (right) as measured by the Ulysses spacecraft (taken from Meyer-Vernet, 2007). The white arrows indicate time progressing in an anti-clockwise fashion. During solar minimum, the coronal magnetic field consists of two oppositely directed large polar coronal holes, which are separated by a so-called helmet streamer, which is a closed magnetic arcade located below the heliospheric current sheet (HCS) that extends into the heliosphere. Because of slight asymmetries and the tilt of the rotation axis of the Sun with respect to the ecliptic plane, in which the planets perform their orbits, there is a typical structure of alternating heliospheric magnetic field (HMF) directions measured at Earth during one rotation around the Sun (see ,e.g., Schwenn and Marsch, 1990). A crossing of the HCS or neutral line is called a sector boundary, while the recurring phases of alternating polarities are known as sector structure. The situation can be visualized as a flying ballerina skirt or like the brim of a sombrero hat as depicted in Figure 1.7 (left). As the band of slow solar 8 1.1. GENERAL OVERVIEW wind aligns roughly with this wavy current sheet, there are also respective recurrent streams of fast wind extending to ecliptical regions, which are azimuthally embedded into otherwise slow solar wind. Due to solar rotation, these streams become radially aligned so that the fast wind runs into the slow one ahead – leading to a compression region – and outruns the trailing one, resulting in a rarefaction region (see right panel of Figure 1.7). The embedded frozen-in magnetic field prevents the fast wind from overtaking and strong enhancements in plasma and magnetic field density at the compression region develop, which generate pairs of waves traveling forward and backward with respect to the stream interface that separates fast and slow wind. As the critical/characteristic speeds decrease with radial distance, these waves steepen into a forward-reverse shock pair, which usually occurs at about 2 AU (Balogh et al., 1999). These kind of structures as a whole are known as corotating interaction regions (CIRs), which have important effects on the propagation of energetic particles because of their associated shocks and turbulence properties. Particularly during the descending and ascending phases in the solar cycle, CIRs can be long-lasting recurrent structures over several solar rotations, constituting the main driver shaping the heliospheric environment. Figure 1.7: Left: Sketch of the wavy heliospheric current sheet (taken from Sakurai, 1985). Right: Sketch of a co-rotating interaction region (taken from Hundhausen, 1972). Towards solar maximum the HCS becomes increasingly tilted until a polarity reversal occurs and it becomes flatter again with decreasing solar activity. The topology of the coronal magnetic field and the associated solar wind is not only rather complex during high solar activity with fast and slow streams at arbitrary latitudes, but there are also a number of impulsive events such as flares that are often associated with coronal mass ejections (CMEs) and distort the above outlined picture of otherwise rather steady configurations during a solar rotation. Both flares and CMEs are 9 CHAPTER 1. INTRODUCTION caused by reconnection events in the solar corona, but there are many more flares than CMEs (e.g., Priest, 2014). Although there are still many open questions, it is believed that if there is an overlying magnetic arch above the reconnection site, there is usually just a flare, whereas a CME occurs when the reconnection site detaches and releases huge amounts of matter. In a flare the released magnetic energy leads to both the emission of short-wavelength radiation and the acceleration of solar energetic particles (SEPs) up to a significant fraction of the speed of the light (∼0.8c), so that there is only little time to make arrangements for the protection of e.g. astronauts from this radiation. In case of a CME, magnetic energy is also channeled into the acceleration of the ejected mass, which can reach velocities of a few thousand km/s and the propagation time to Earth can be as short as a day or two. CMEs often drive a shock that can also accelerate SEPs. Figure 1.8: Artistic sketch of a CME set loose at the Sun and propagating towards Earth with its geomagnetic field. [www.land-of-kain.de/docs/spaceweather/nasa spaceweather small.jpg]. The geomagnetic field of the Earth is usually able to deflect the solar wind around the Earth’s magnetosphere, but CMEs can strongly distort this protective shield and particles can reach Earth especially from above the poles (see Figure 1.8). In the atmosphere these particles ionize oxygen and nitrogen atoms, leading to the beautiful phenomenon of the aurora, but in some cases particles can also reach the Earth’s 10 1.1. GENERAL OVERVIEW surface (ground level events) and cause shortcuts and other damage in electrical systems. As an ever increasingly technology-dependent society, it is of vital importance to us to be able to understand and predict these kind of events. Respective efforts are made by scientists over the world and are channeled in institutions such as the Space Weather Prediction Center (SWPC, http://www.swpc.noaa.gov/) and the Coordinated Community Modeling Center (CCMC, http://ccmc.gsfc.nasa.gov/). Turbulence in the solar wind Fluctuations that are super-imposed onto the large-scale magnetic field, solar wind velocity and density, have already been seen in the very first solar wind measurements, and respective studies (e.g., Coleman, 1968; Belcher and Davis, 1971) showed that the solar wind has to be considered as an inherently turbulent plasma. Turbulence is not only generated in the corona from where fluctuations are convected into interplanetary space, but also in-situ by the plasma itself at density gradients throughout the heliosphere, as well as on transients with their associated shocks such as CIRs and CMEs and by the isotropization of newly born pickup ions (see next section). Turbulence has several important effects (see the review by Miesch et al., 2015). For example, the associated turbulent energy decays with time and radial distance and heats the solar wind plasma particularly in the outer heliosphere (as measured by Voyager 2, Gazis et al., 1994)) due to the increased pickup ion production. Furthermore, turbulence is a catalyst for reconnection (e.g. Yokoi et al., 2013). Another important factor is the influence of turbulence on the propagation of energetic particles: the turbulent magnetic field guides their transport, where the particles not only gyrate about and follow the mean magnetic field lines but are also able to perform a perpendicular diffusion towards neighboring field lines due to the turbulent fluctuations. Respective models for the transport of energetic particles in the heliosphere have – for a long time – been coupled to the solar wind dynamics by parameterizing all relevant processes in terms of the large-scale background quantities, in particular for the treatment of spatial diffusion (see the review of Potgieter, 2013). Only recently have models been developed that formulate the diffusion coefficients explicitly based on the small-scale fluctuations (ab-initio models, e.g., Engelbrecht and Burger, 2013), which has been made possible by new improvements in analyzing turbulence in solar wind measurements (e.g. Horbury and Osman, 2008; Horbury et al., 2012) on the one hand, and respective modeling on the other hand, as described next. After the pioneering papers by Tu et al. (1984) and Zhou and Matthaeus (1990) and first systematic studies of the radial evolution turbulence quantities (Zank et al., 1996) major progress has only been achieved during the recent decade. First, Breech et al. (2008) improved the previous modelling by considering off-ecliptic latitudes. 11 CHAPTER 1. INTRODUCTION Second, in Usmanov et al. (2011) this turbulence model was extended to full time dependence and three spatial dimensions. In the same paper the turbulence evolution equations were simultaneously solved along with a large-scale MHD model of the supersonic solar wind. While Usmanov et al. (2012, 2014) extended this study further with the incorporation of pick-up ions and electrons as separate fluids, Oughton et al. (2011) undertook another extension of the modelling by considering not only one but two, mutually interacting, incompressible components, namely quasi-twodimensional turbulent and wave-like fluctuations. All of the mentioned models have one limitation in common, namely the condition that the Alfvén speed is significantly lower than that of the solar wind, which precludes their application to the heliocentric distance range less than about 0.3 AU. The turbulence in this region and its consequences for the dynamics of the solar wind has been studied on the basis of the somewhat simplified transport equations for the wave power spectrum and the wave pressure (e.g., Tu and Marsch, 1995; Hu et al., 1999; Vainio et al., 2003; Shergelashvili and Fichtner, 2012). The Alfvén speed limitation in the more detailed models was very recently removed with a new, comprehensive approach to the general problem of the transport of low-frequency turbulence in astrophysical flows by Zank et al. (2012), see also Dosch et al. (2013) and Adhikari et al. (2015). Not only are the derived equations formally valid close to the Sun in the sub-Alfvénic regime of the solar wind, but they also represent an extension of the treatment of turbulence by non-parametrically and quantitatively considering the evolution of the so-called energy difference in velocity and magnetic field fluctuations, and by explicitly describing correlation lengths for sunward and anti-sunward propagating fluctuations as well as for the energy difference (sometimes also referred to as the residual energy). 1.1.3 The heliosphere The heliosphere is formed as a consequence of the interaction between the solar wind and the interstellar medium (ISM), through which the Sun moves on its 200 million year lasting orbit around the Galactic center with a speed of 240 km/s, while the relative speed between the Sun and the ISM is only about 25 km/s. From the point of view of the Sun this looks like an inflow of interstellar plasma from the so-called nose region, which leads to a number of interesting structures in the heliosphere (see Figure 1.9): Both the ISM and the solar wind plasma have embedded magnetic fields and are unable to permeate the other under idealized conditions. The boundary between the two plasmas is the bullet-shaped heliopause (HP), a tangential discontinuity at which the solar wind pressure equals the one of the ISM. The radially out-flowing solar wind adapts to this obstacle by first transiting to a heated, sub-sonic state at the termination shock (TS), and in the region between the TS and the HP – the inner heliosheath – the flow is bent towards the heliospheric tail. Similarly, the 12 1.1. GENERAL OVERVIEW Figure 1.9: Overview of the heliospheric environment: The solar wind emanating from the Sun (at center) blows a ”bubble”, the heliosphere, into the interstellar medium. Typical structures are the termination shock and the heliopause (see text). The trajectories of mankind´s most distant spacecraft, Voyager 1 and 2, are also indicated. [www.nasa.gov/vision/universe/solarsystem/voyager agu.html]. ISM flow is also decelerated, resulting in a bow wave or bow shock, depending on whether the ISM speed is supersonic or not, which is still under debate (see Scherer and Fichtner, 2014). While the solar wind plasma can be considered as fully ionized due to the high coronal temperatures, the ISM plasma is a mixture of ions and electrons, as well as neutral atoms (i.e. gas, 99%), but also dust (1%, Boulanger et al., 2000)). As the interstellar plasma slows down in the outer heliosheath’s nose region and thereby becomes denser, the cross section for charge exchange processes, in which an electron from a neutral atom ”jumps” to an ion, increases strongly. This leads to the formation of the so-called hydrogen wall, which is an accumulation of neutral hydrogen that can be observed in the Lyman-α absorption line (Wood et al., 2000). The majority of neutral atoms, however, can penetrate the heliopause and will eventually (at least beyond about 8 AU) be ionized by photo-ionization, electron impact or charge-exchange with solar wind particles (see, e.g., Scherer et al., 2014), so that they become so-called pickup ions, because they are then assimilated in and convected with the solar wind flow. During this process, in which the initially 13 CHAPTER 1. INTRODUCTION ring-distributed pickup ions isotropize, Alfvén waves are generated that contribute to the turbulence spectrum in the outer heliosphere, eventually heating the latter. Furthermore, pickup ions are considered as seed population for the generation of anomalous cosmic rays accelerated at the TS (the term ”anomalous” stems from the different composition – higher abundances of heavier elements – of these cosmic rays as compared to Galactic Cosmic Rays (GCRs) or SEPs, which are mainly protons, Cummings et al., 1995). If pickup ions are generated in a charge exchange process the former solar wind ion becomes an energetic neutral atom (ENA) and can flow unimpeded by the magnetic field. These atoms are measured by the IBEX spacecraft orbiting Earth and mapping the whole sky (McComas et al., 2014). An unpredicted phenomenon that was found is the so-called IBEX ribbon of enhanced ENA flux that aligns with the plane where the interstellar magnetic field is perpendicular to the heliocentric radius vector. While an undisputed theory has not yet been found to explain the ribbon (see the review of McComas et al., 2011), the IBEX mission provides insights into outer-heliospheric processes and structures by global remote observations, as opposed to the few in-situ spacecraft measurements described next. Mankind´s most distant spacecraft – Voyager 1 and 2 (V1, V2) – are both headed towards the upwind region of the heliosphere, but offset from the ecliptic by about 30◦ towards north and south, respectively. V1 crossed the TS in 2004 at 94 AU, whereas V2 crossed it at 84 AU three years later. V1 is the first and of yet only spacecraft that has also left the heliosphere in crossing the heliopause in 2012 at about 122 AU, when particles of solar origin became absent while the GCR intensity went up. The crossing came with some unexpected findings: first, the inner-heliosheath thickness was expected to be much larger – hinting at energy losses there reducing the pressure (Izmodenov et al., 2014), while it was argued by Fisk and Gloeckler (2014) that the heliopause has not been reached yet after all. Second, the magnetic field did not change much after V1 crossed the heliopause, where it was expected to change from a Parker spiral field of solar origin to the one of the ISM that drapes around the heliopause (see Röken et al., 2015, for a recent model, visualizations and further references). This mystery is still under investigation, but there are ideas that, for example, involve heliopause instabilities or reconnection, yielding a rather turbulent transition (Grygorczuk et al., 2014; Opher and Drake, 2013). V2 is expected to cross the heliopause soon, which will hopefully shed some light on the unresolved issues as this spacecraft – unlike V1 – still benefits from a working plasma instrument. Both V1 and V2 have enough power left to operate until about 2020, which will most likely, however, not be enough time to reach the bow wave/shock expected at about 230 AU. The heliospheric tail is largely uncharted and no current spaceprobe is scheduled to investigate. Its extent is probably a few thousand AU, and most models predict a single tail that can ”store” several solar cycles with their respective polarity changes (e.g., Pogorelov et al., 2014). However, it was recently proposed by Opher et al. (2015) (reviving an idea by Yu, 1974) that there is a split tail caused by the twisted HMF that channels two separate jets. 14 1.1. GENERAL OVERVIEW It should be stressed though that the distances to the heliospheric structures as given above are very likely to be time- and direction-dependent: on the one hand the solar wind conditions change during the solar cycle and transient structures such as CMEs and CIRs merge and interact with increasing radial distance to form global merged interaction regions (GMIRs) with enhanced pressure, which probably make the TS constantly move back and forth around a mean heliospheric distance of about 90 AU (le Roux and Fichtner, 1999). On the other hand the ISM is also not homogeneous (e.g. de Avillez et al., 2015) and changes the overall size of the heliosphere, albeit on much longer time-scales. Under extreme conditions it is possible that during the Sun’s passage through a hot and dense interstellar cloud the TS distance can shrink down to be within the planetary orbits (Scherer et al., 2008). The heliosphere also acts as a shield from highly energetic GCRs, whose influence on the Earth’s climate is still under debate and probably not negligible, so that during the Sun’s orbit around the galactic center there have been phases of more intense impacts of GCRs on Earth (Scherer et al., 2006; Kopp et al., 2014). The sources of GCRs are believed to be energetic phenomena such as jets driven by black holes or supernovae remnants (e.g., Büsching et al., 2005), which can accelerate particles to very high energies. The exact origin of such a particle that might eventually end up in a detector near Earth cannot be pinpointed because the transport of GCRs is mainly along the poorly known (inter-)galactic magnetic field. Meanwhile, the modulation of the cosmic ray spectrum in the heliosphere (and possible other astrospheres) and the resulting cosmic ray flux at Earth are more directly accessible due to an increasingly better understanding of the heliospheric structures. 1.1.4 Stellar winds and astrospheres It was realized already in the late middle ages by Tycho Brahe that stars are not perfect, ever constant celestial objects, but suffer from mass loss (although at that time it was not the stellar wind that was observed, but most likely the more spectacular novae), and it was much later also established that respective line profiles in the continuum of stars could be explained by expanding stellar atmospheres (Beals, 1929) long before the confirmation of the existence of the solar wind. The term stellar wind, however, was again coined by Parker (1960). Recent observational advancements (e.g. Hubble space telescope, or the Wide-field Infrared Survey Explorer (WISE) mission) allowed for all-sky surveys measuring dust emission (mostly infrared, e.g., Groenewegen et al., 2011, see Figure 1.10) or the Lyman-α line (Linsky and Wood, 2014) associated with bow shocks around other stars, thus indicating the existence of respective astrospheres. Therefore, many findings about the solar wind and our heliosphere can be transferred to other stars and their astrospheres. However, the interstellar environment around those stars is most likely different from our local environment, and the processes driving a stellar wind can be very different from the coronal solar wind (see Lamers and Cassinelli (1999) for an overview). For 15 CHAPTER 1. INTRODUCTION Figure 1.10: Infrared images from the Herschel MESS key program (Groenewegen et al., 2011) of stars with bow shocks. [http://ster.kuleuven.be/~nick/research astrospheres.html]. example, hot (O/B-) stars are associated with radiation or line driven winds by the absorption in spectral lines. These stars posses a completely radiative interior, so that the strength or even the existence of their magnetic fields is highly debated (although some examples have recently been found (Hubrig et al., 2015)). Line driven winds can have terminal velocities of several thousand km/s leading to huge astrospheres with extents of several parsecs. Surprisingly, such high-velocity winds are also found to be generated by the much cooler M-dwarf stars, where the interior energy transfer is completely convective and are, therefore, magnetically very active. It can be speculated that a considerable amount of the mass loss is probably attributable to flares and CMEs. The astrospheres around hot stars are of interest, because they represent huge cavities in the ISM that modulate the cosmic ray flux through them and might therefore be responsible for tiny-scale anisotropies in respective all-sky maps of the GCR flux (Scherer et al., 2015). For cool stars that are relatively close-by, it is even possible to infer their photospheric magnetic field distribution (Donati et al., 2008), which allows to apply solar and heliospheric numerical models to such stars’ stellar winds as performed by, e.g.,Vidotto et al. (2011, 2015); Vidotto (2014). A particular focus lies on M-dwarfs with orbiting exoplanets, and the possible interaction between these stellar winds with planetary magnetospheres (e.g., Ip et al., 2004; Preusse et al., 2007), which should be considered as one of the important factors for habitability. As discussed for the heliosphere above, the size of a respective astrosphere and its 16 1.2. MODEL FORMULATION AND NUMERICAL SETUP capability to shield planets from cosmic rays may also be important, so that global astrospheric models for these cases are also planned. 1.2 1.2.1 Model formulation and numerical setup Principal simplifications From a modeling point of view, it is customary and necessary to sub-divide the heliospheric environment because of the very different physical processes and associated length and time scales involved in one region or the other. Roughly speaking, these regions can be thought of as spherical shells of increasing radial distance from the center of the Sun as outlined in the preceding sections, so that – for an outward flow of information – a particular model of one region usually relies on the region below and determines the following one. In particular, this thesis is primarily concerned with the formation of the solar wind in the Sun’s lower atmosphere and its continued dynamical behavior in the supersonic part of the heliosphere, i.e. within the TS. Following the above reasoning, such a model naturally depends on processes in the solar interior. It is, however, simpler to use the directly accessible observational data from the solar photosphere as a starting point. For large-scale models the knowledge of the photospheric magnetic field suffices to determine the coronal magnetic field, which shapes the emanating solar wind filling the heliosphere. Important processes in the super-Alfvénic heliosphere above the corona primarily involve the interaction between differently fast streams, which leads to interesting structures such as evolving shocks, which have important implications for, e.g., the transport of energetic particles in the heliosphere, or planetary magnetospheres. With increasing distance from the Sun, the initially rather inhomogeneous solar wind perturbed by CIRs and CMEs – which eventually smear out in GMIRs – becomes relatively homogeneous. Such (time-dependent) conditions can then be used as input to outer-heliospheric models (as done in, e.g., Pogorelov et al., 2014), which is particularly important for a better understanding of the Voyager 1 measurements. The above mentioned necessity for subdividing the heliosphere is also found in contemporary modeling of the formation of the solar wind: typically one distinguishes the sub-Alfvénic corona and the super-Alfvénic part beyond the so-called heliobase (Zhao and Hoeksema, 2010) at about 20R . There are two common approaches to compute the large-scale coronal magnetic field, namely (i) full three-dimensional (3D) MHD models and (ii) potential field models (see Chapters 2 and 3 for details). While the full MHD models are principally better suited to address all the involved physics self-consistently, they have the drawback of being quite complex and requiring large computer resources. Particularly this last issue has made potential field models a popular approach to describe the large-scale coronal magnetic field, which has also been shown to be remarkably similar in both MHD and potential field 17 CHAPTER 1. INTRODUCTION models (Riley et al., 2006). Both approaches can be applied to observational data of the photospheric magnetic field (i.e. magnetograms) for usage as inner boundary conditions, and both eventually provide the global coronal magnetic field. However, the drawback of the potential field models is that they do not self-consistently compute the resulting solar wind speeds. Instead, empirical models have to be used like the so-called Wang-Sheeley-Arge (WSA; Arge and Pizzo, 2000) model, which uses the coronal magnetic field topology to determine solar wind speeds, i.e. fast wind from coronal holes and slow wind from their boundaries as outlined in Section 1.1.2. Thus, both MHD and empirically based models of the corona provide the resulting magnetic field and solar wind speeds at the heliobase, which can be further utilized in subsequent MHD simulations to investigate the interaction of different solar wind streams. Compared to coronal MHD simulations, these inner-heliospheric MHD simulations are not requiring as much spatial and time resolution, and are more simple to formulate, because the solar wind is super-Alfvénic here and information cannot travel backwards – an issue that otherwise complicates boundary conditions considerably. Large-scale heliospheric models are usually based on a fluid treatment of the plasma, so that the latter is described by the MHD equations. These equations are only valid under certain conditions (see, e.g., Priest, 2014), which are not well fulfilled everywhere in the heliosphere (such as the assumptions of distribution functions being close to Maxwellian despite the non-collisionality in the solar wind). There are still several reasons to pursue modeling approaches based on MHD: first, alternative descriptions either based on non-Maxwellian distribution functions (meso-scale) or on particle trajectories themselves (micro-scale) are too fine-scaled for a numerical modeling of the large-scale heliospheric structures due to computational limitations. Such approaches (e.g. particle in cell (PIC) simulations) are instead applied to much smaller computational domains (e.g., Kunz et al., 2014). Second, the fluid picture is still a reasonable good approximation as many respective studies (see, e.g., the review by Zank, 1999) were successful in recreating solar wind properties as compared to spacecraft data. Consequently, this approach is also taken in this thesis by using the state-of-the-art MHD code Cronos, which is described in more detail in the following section. 1.2.2 The Cronos code The equations of MHD can be solved analytically for only a few simplified cases so that a numerical approach has to be taken. This section gives a descriptive overview of the Cronos code that is used in this work. The code was developed by R. Kissmann (based on Grauer et al., 1998; Kleimann et al., 2009), and was used in recent years primarily by his current working group at the University of Innsbruck as well as by his former group in Bochum. Studies focused on astrophysical [ISM turbulence (Kissmann et al., 2008), accretion disks (Kissmann et al., 2011), colliding 18 1.2. MODEL FORMULATION AND NUMERICAL SETUP wind binaries (Reitberger et al., 2014b,a), and astrospheres (Scherer et al., 2015)] as well as heliospherical scenarios [CMEs (Kleimann et al., 2009), magnetic clouds (Dalakishvili et al., 2011), background solar wind (Wiengarten et al., 2013, 2014), and turbulence in the solar wind (Wiengarten et al., 2015)]. The following text summarizes the main features of Cronos without going into details on numerical issues, but see Appendix C of Wiengarten et al. (2015) and references therein. In its basic setup the code solves the equations of ideal MHD, given here in normalized form: ∂t ρ + ∇ · (ρu) ∂t (ρu) + ∇ · ρuu + (p + kBk2 /2) 1 − BB ∂t e + ∇ · (e + p + kBk2 /2) u − (u · B)B ∂t B + ∇ × E = = = = 0 0 0 0 (1.1) (1.2) (1.3) (1.4) where ρ is the mass density, u is the velocity of a fluid element, B and E describe the electromagnetic fields, e is the total energy density, and p is the scalar gas pressure. Furthermore, 1 denotes the unit tensor, and the dyadic product is used. So far the system is under-determined, and requires the closure relations ρ kuk2 kBk2 p/(γ − 1) : γ 6= 1 + + (1.5) e = 0 : γ=1 2 2 E+u×B = 0 (1.6) ∇·B = 0 (1.7) to complete the set of equations. Equation (1.5) for the total energy density e distinguishes between adiabatic (γ 6= 1) and isothermal (γ = 1), where in the latter case Equation (1.3) becomes redundant and does not enter the computations. The Cronos code is written in C++. It has been designed in a user-friendly fashion such that the part in which the user adopts an own specific physical model is strictly separated from the part implementing the algorithms of the underlying numerical schemes, i.e. the code is modularized. This also allows for extensions of the basic equations of ideal MHD, such as additional source or flux terms, as well as the incorporation of equations for additional quantities to be solved along with the original ones. Cronos employs the semi-discrete finite-volume scheme CWENO (Kurganov and Levy, 2000) and integrates the equations forward in time by means of the RungeKutta method of adaptable accuracy. The initial time step has to be chosen small enough, so that the Courant-Friedrichs-Lewy condition (CFL; Courant et al., 1928) is fulfilled. Afterwards, adaptive time-stepping is applied such that a subsequent time step is chosen to be the largest possible one fulfilling CFL. The solenoidality condition of the magnetic field (1.7) is ensured by a constrained transport scheme. 19 CHAPTER 1. INTRODUCTION The coordinate systems supported by Cronos are currently Cartesian, cylindrical, and spherical, where the latter have recently been extended to cover the coordinate singularities (e.g. the polar axis in spherical coordinates) as well. As part of this thesis, non-linearly spaced grids have also been implemented (see Appendix A), but the general description here concentrates on linearly spaced grids for simplicity. The simulation is carried out on a three-dimensional grid, which in C++ is saved in an array whose (integer) indices i,j,k run from 0 to (Nx,y,z − 1). Thus, the grid comprises Nx · Ny · Nz cells where Nx,y,z stands for the number of cells for the corresponding dimension. To construct a grid, the user first specifies the desired range for each dimension, either by fixing the first and last cell center locations or the first cells left and last cells right face coordinate. Here, the first possibility is described, so that, considering only the x direction, xs and xe denote the start and end cells’ centers, and the length of the interval xl becomes xl = xe − xs . For a linearly spaced grid the cell size ∆x is then given by xl (1.8) ∆x = Nx − 1 and, accordingly, the coordinate of a cell center can be calculated via xs + i∆x ri,j,k = ys + j∆y . zs + k∆z (1.9) The grid can be thought of as being composed of rectangular cells with the grid points at the center. In Cronos the fluid variables (ρ, vx , vy , vz , e, p) are computed at the cells’ centers so that the corresponding arrays, in which they are stored, have the same indices as the array storing the grid points. For the magnetic field a constrained transport scheme (Brackbill and Barnes, 1980) is used to ensure conservation of magnetic divergence to machine precision. To compute the flux going in and out of the cells, the magnetic field components at the corresponding orthogonal surfaces of each cell are needed. These would need semi-integer indices, which is not possible in C++. Therefore, the arrays storing magnetic field components are offset from those storing the cell center coordinates by one half, which has to be kept in mind. For example, when computing the absolute magnetic field strength |B| at a cell center, averaging the components of B from the corresponding adjacent cell surfaces is required, as can be seen from the layout of the grid shown in Figure 1.11. Cronos makes use of so-called ghost cells (also depicted in Figure 1.11), which are an extension of the actual grid at its boundaries. They are used to implement boundary conditions that can be defined for all six boundaries (two boundaries for each of the three dimensions). Here, variables are not calculated according to the set of equations and instead are given values according to the applied boundary condition. Pre-implemented options for boundary conditions are Periodic, Extrapolation, Outflow, and Reflecting, which are self-explanatory. Additionally, the user 20 1.2. MODEL FORMULATION AND NUMERICAL SETUP Figure 1.11: Sketch of the grid layout: grid points are cell-centered (integer indices) and store fluid variables (abbreviated as Q), while magnetic field components are stored at the corresponding adjacent cell surfaces (semi-integer indices in the corresponding direction, integer indices for the others). The concept of ghost cells is also demonstrated as they can be considered as an extension at the simulation box’ faces with negative indices (or indices beyond N at the outer boundary); y indices were chosen arbitrarily. has the opportunity to implement own boundary conditions where necessary. If coordinate singularities are detected, the preset boundary conditions are ignored and a respective scheme is applied automatically (see Appendix A). At the beginning of a simulation the MHD variables have to be initialized (i.e. they are given values according to the model being simulated). This also covers the ghost cells, so that the initialization process goes hand in hand with the implementation of boundary conditions, especially if a boundary condition is such that it holds on to the initialized values. Special care must be taken when initializing the magnetic field, which has to be done in a divergence-free manner as otherwise the constraint transport algorithm maintains non-zero divergences. A comfortable way of achieving this is to use the vector potential A, which by definition gives ∇ · B = ∇ · (∇ × A) = 0 . (1.10) In case a vector potential is unavailable for complicated magnetic fields, other methods must be devised, such as the one developed in this thesis (see Section 2.4.2). The integration proceeds up to a predefined point in time at which the simulation ends. Setting this point requires some experience and depends on the physical model: if a steady-state is sought from time-independent inner boundary-conditions and outward flow of information, the simulation should end by the time a converged solution has been found. This has to be done manually, but can be estimated by the propagation time of the slowest plasma parcels to the outer boundary. In other cases – investigating time-dependent phenomena – this consideration can be more complex. 21 CHAPTER 1. INTRODUCTION It is possible to let the code run in parallel using MPI (Message Parsing Interpreter) on stand-alone computers using multiple cores, but also on several computers linked by a network connection. Output is saved in the HDF5 (Hierarchical Data Format; http://www.hdfgroup.org/HDF5/) format, which can be written to individual files by each core (and have to be reassembled into a single file afterward), but also in a parallel fashion into a single file. Output can be set to be done at fixed intervals in time and should be limited to actual needs (output can amount to several Gigabytes): in simulations converging towards steady-state the initial and final state may suffice, but output is also important for debugging when errors are encountered. It is also possible to create movies of a simulation, which requires a lot of time steps with output. For analysis and visualization, the files can later be read using plotting tools such as SectorPlot written by J. Kleimann in IDL (interactive data language1 ). It allows the user to choose quantities to be visualized from different output time steps of a selected project. Besides several other options, the user has to select a (2D) cross section to be plotted since 3D visualizations are difficult to handle and are not implemented here. Alternative tools for visualizing results also in (pseudo-) 3D are Visit2 and Paraview3 . 1.3 Outline of the thesis This section puts the work described in this thesis in perspective with respect to the previously outlined interplay of the protagonists in the heliosphere. Details are postponed to the subsequent chapters. The building blocks of the present thesis are four accepted publications (three as first author) in peer-reviewed astrophysical journals, as well as some recent work, which will be published in the near future. For each of these there is a separate chapter, respectively starting with a summarizing introduction, followed by a copy of the actual publication – sometimes appended with additional details – and concluded with follow-up thoughts and developments leading to the subsequent chapters. Chapter 2 is centered around my initial first-author publication MHD simulation of the inner-heliospheric magnetic field (Wiengarten et al., 2013), where I pursuit the goal to implement observations-based inner boundary conditions to drive 3D simulations of the inner heliosphere (0.1 to 1 AU) in the Cronos framework. This required several extensions to the code such as introducing the effects of solar rotation and solving the full energy equation (instead of previously adopted isothermal 1 http://www.ittvis.com/language/en-us/productsservices/idl.aspx https://wci.llnl.gov/simulation/computer-codes/visit/ 3 http://www.paraview.org 2 22 1.3. OUTLINE OF THE THESIS models (Kleimann et al., 2009; Wiengarten, 2011)). The respective implementations were validated by reproducing the analytic models of Parker (1958) and Weber and Davis (1967). These latter models have been quite successful in describing general features of the solar wind and the HMF, but are too simple for reproducing more complex conditions (e.g. wavy current sheets, CIRs) as measured by spacecraft. To that end, in a collaboration with colleagues from the Max Planck Institute for Solar System Research (MPS) in Katlenburg-Lindau (now in Göttingen), maps of the radial magnetic field at the heliobase for solar minimum and maximum conditions served as inner boundary conditions for my subsequent MHD simulations. These maps are the result of a coupled solar surface flux transport (SFT) and a coronal potential field model as described in Jiang et al. (2010). The common approach in potential field models is to use magnetogram data as inner boundary conditions, which has the drawback that the data is composed of measurements taken during one solar rotation, so that especially the far side of the Sun is described by somewhat outdated data. SFT models aim to model the photospheric motions and evolution of sunspots to allow for a better treatment of the far side of the Sun. Moreover, by using sunspot group data, which is available for longer times than modern magnetograms, the model can recreate conditions from before the space age. The subsequent utilization of these results as input to MHD simulations therefore provided a novel possibility to model the accompanying solar wind conditions. As the Cronos code is usually operated with analytic prescriptions for initial and inner boundary conditions, it was necessary to extend the code by devising routines to read-in and interpolate the input data on arbitrary grids, as well as to find a method for a divergence-free initialization of the magnetic field. Furthermore, inner boundary conditions for the remaining MHD quantities had to be set, which follow from the coronal magnetic field topology (fast wind from coronal holes, slow wind from their boundaries (Wang and Sheeley, 1990)), where I adopted an empirical model of Detman et al. (2011), based on the WSA model (Arge and Pizzo, 2000; Arge et al., 2003). MHD simulations were performed for typical solar minimum and maximum conditions, and the results I obtained appeared reasonable at the time as they were roughly compared to spacecraft data and agreed with the general topologies expected for the respective phases in the solar cycle. A subsequent detailed comparison with spacecraft data after publication revealed that the input data used so far were insufficient, so that a different approach for the coronal modeling became necessary and was subsequently taken as described next. The work described in Chapter 3 was done in the framework of a joint project with the University of Kiel and the North-West University, Campus Potchefstroom, South Africa, on the influence of CIRs on the propagation of energetic particles, where the respective CIR disturbed background solar wind resulting from my MHD simulations is used as input to a respective stochastic differential equation (SDE) model for the transport of energetic particles in the heliosphere. My MHD modeling 23 CHAPTER 1. INTRODUCTION is described in Cosmic Ray transport in heliospheric magnetic structures. I. Modeling background solar wind using the Cronos magnetohydrodynamic code (Wiengarten et al., 2014), while the subsequent SDE modeling paper is in preparation, but see Kopp et al. (2012) for a detailed description of the idea behind the SDE modeling as well as the respectively developed code. Once again, a test case for validating the ability of the Cronos code to correctly compute CIR structures was performed by reproducing the numerical results of Pizzo (1982). After successful validation, a time period in the ascending phase of solar cycle 23 (August 2007) – during which a stable CIR configuration was present – was chosen for a subsequent detailed comparison with spacecraft measurements in the ecliptic at 1 AU (Stereo-A/B; see Kaiser, 2005) and out-of-ecliptic (Ulysses). The required potential field modeling of the solar corona providing the boundary conditions for my simulations had to be performed in a different fashion as compared to my previous work, because the input data used there were found to be incorrect. The employed WSA model for deriving the resulting solar wind speeds is tuned to the usage of GONG (Global Oscillation Network Group) magnetograms in coupled potential field source surface (PFSS) (Altschuler and Newkirk, 1969) and Schatten current sheet models (Schatten, 1971), so that these were adopted. Instead of using the common spherical harmonics approach to solve the Laplace equation in these models, a finite difference approach was taken and for the first time coupled to the WSA formalism. For the analysis of the resulting coronal field topology a field line tracing algorithm was devised, which is also required to derive the resulting solar wind speeds. The results showed good agreement with spacecraft data allowing to expect insights for energetic particle propagation from the ongoing SDE modeling. For the latter, transport coefficients are employed that actually depend not only on the large-scale MHD quantities, but also on the fluctuations of the magnetic field. These cannot be resolved in large-scale MHD simulations, but there are models describing the evolution of integral turbulence quantities such as the total turbulent energy density and the ratio between inward and outward propagating modes. In order to eventually be able to provide these quantities for an improved SDE modeling, a turbulence transport model was implemented alongside the large-scale MHD equations in the Cronos framework, which is the topic of Chapter 4 based on my publication Implementing turbulence transport in the CRONOS framework and application to the propagation of CMEs (Wiengarten et al., 2015). I started again by validating the implementation by comparing with previous work by Usmanov et al. (2011), during which some discrepancies were found that could be traced back to errors made by these authors (A. Usmanov, priv. comm). A corrected model was used, which was further extended by removing the constraint of being only applicable in regions of highly super-Alfvénic solar wind (U VA ), while I had the goal to have a model that can eventually be coupled with my previous WSA driven MHD simulations that start in regions where U ≥ VA . This was achieved by simplifying a more general model of turbulence transport by Zank et al. (2012). The resulting model was then applied for the first time to the propagation of CMEs, 24 1.3. OUTLINE OF THE THESIS and although there is little observational data the obtained results for the turbulence levels are in general agreement with estimates by Subramanian et al. (2009). Furthermore, it is found that, on the one hand, turbulence does not act back strongly on the large-scale quantities, but, on the other hand, that the CME is a strong driver of turbulence. Therefore, such a model is not required for studies focusing on large-scale CME quantities only, but it does provide another interesting scenario for a related energetic particle propagation study. Chapter 5 describes the subsequent implementation of the complete model of Zank et al. (2012), which thus far had been solved neither in 3D nor self-consistently coupled with the MHD equations. A simplified implementation in the ecliptic plane only and neglecting some of the involved terms was done by Adhikari et al. (2015), whose results were used for validation purposes. Afterward, the full model is solved, which gives some yet not fully understood results. Consequently, this work is unpublished as of yet, but ongoing. Another line of work in our team is the modeling of the large-scale heliosphere and astrospheres, i.e. the interaction between the solar/stellar wind and the interstellar medium. An analytic model for the interstellar magnetic field in the vicinity of the heliopause was recently presented by Röken et al. (2015), in which it was also compared to numerical results obtained with the Cronos code by J. Kleimann. This numerical setup was altered by K. Scherer and myself for the application to the astrosphere of the O-star λ Cephei. Due to the large terminal velocities of this star’s stellar wind, its astrosphere is a cavity of several parsecs in diameter and is thus able to modulate high-energy cosmic rays that pass through it, which could be responsible for tiny-scale anisotropies in all-sky maps of the cosmic ray flux as measured, e.g., by IceCube. This is the topic of Chapter 6 based on the paper Cosmic rays in astrospheres (Scherer et al., 2015). 25 CHAPTER 1. INTRODUCTION 26 Chapter 2 MHD simulation of the inner-heliospheric magnetic field 2.1 Overview On a large scale the unperturbed HMF inside the termination shock is remarkably well described by the Archimedean spiral pattern frozen-into the solar wind as proposed first by Parker (1958). Close to the Sun, where the magnetic field is energetically dominant as compared to the kinetic energy density of the solar wind, the latter follows the magnetic field anchored on the rotating Sun, leading to a transfer of angular momentum to the solar wind, which has slowed the Sun’s solar rotation over its lifetime (Weber and Davis, 1967). The photospheric motion of these anchored footpoints motivated yet another analytic prescription of the resulting HMF (Fisk, 1996), which in this model also has a latitudinal component (see Figure 2.1). Several adaptions have been proposed by e.g. Jokipii and Kota (1989); Smith and Bieber (1991); Schwadron (2002); Burger and Hitge (2004); Burger et al. (2008); Schwadron et al. (2008) that have been compared by Wiengarten (2009) and Scherer et al. (2010). These analytic models can account for the large-scale behavior of the HMF particularly during solar quiet times. However, for a more quantitative and realistic modeling MHD simulations have to be performed. It is customary to distinguish the sub-Alfvénic corona from the super-Alfvénic part beyond the so-called heliobase (Zhao and Hoeksema, 2010) at about 20R , with two common approaches to compute the coronal magnetic field: on the one hand, modern computers have made it possible to use fully 3D MHD codes to model the corona, thereby directly taking into account the effect of the evolving solar wind (e.g., Usmanov and Goldstein, 2003; Cohen et al., 2007; Feng et al., 2010; Riley et al., 2011; van der Holst et al., 2014, and many more). These models account for more realistic physics, but require enormous amounts of computing power and are of considerable complexity as they have to address the coronal heating problem and respective acceleration of the solar wind, as well as the tran27 CHAPTER 2. MHD SIMULATION OF THE INNER-HELIOSPHERIC MAGNETIC FIELD Figure 2.1: Field lines in the Parker (left) and Fisk model (right) ( taken from Sternal et al. (2007)). sition from sub- to super-Alfvénic conditions, which leads to a severe complication of boundary conditions as information can travel back towards the inner boundary (Nakagawa, 1981). On the other hand, potential field models for the coronal magnetic field (Altschuler and Newkirk, 1969; Schatten et al., 1969) have been used for decades already, as they are computationally cheap and were shown to yield remarkably similar results to respective MHD models (Riley et al., 2006). Such potential field models assume the corona to be current-free (∇ × B = J = 0) that together with ∇ · B = 0 allows to compute the coronal magnetic field by solving the Laplace equation for a scalar potential ∆Ψ = 0. While the imposed inner boundary condition is a magnetogram, the outer boundary condition is to force the field to become radial at the so-called source surface (typically at 2.5R ), which is to mimic the effect of the solar wind eventually dragging the magnetic field. An example for a PFSS solution within the source surface radius as presented in Altschuler and Newkirk (1969) is shown in Figure 2.2 for solar minimum conditions in 1966. More elaborate potential field models have also been proposed, e.g. to realize thin and sharp current sheets (Schatten, 1971) or for a better treatment of the opening of field lines by introducing an intermediate cusp surface (Bogdan and Low, 1986). Another more elaborate approach can also be taken via force-free models assuming B k J (see, e.g., the review by Wiegelmann et al., 2014). The resulting solar wind speeds that are not inherently taken into account in potential field models can, however, be deduced from the coronal magnetic field topology: as found by Wang and Sheeley (1990), there is an inverse relationship between the topology-characterizing flux tube expansion factor and the respective solar wind 28 2.1. OVERVIEW Figure 2.2: Example of a PFSS solution for the coronal magnetic field within the source surface at 2.5R (taken from Altschuler and Newkirk (1969)). speed at the heliobase – which basically quantified the association of fast wind with coronal holes and slow wind with coronal hole boundaries. This finding motivated the WSA model (Arge and Pizzo, 2000; Arge et al., 2003) as a space weather forecasting tool, estimating the solar wind speed at Earth either by means of simple kinematic propagation (e.g., the HAF model Akasofu et al., 1983) or by MHD simulations (e.g. Odstrčil et al., 2004). Respective frameworks have transitioned to operation at the space weather prediction center (SWPC) and the coordinated community modeling center (CCMC). In my following initial first-author publication MHD simulation of the inner-heliospheric magnetic field (Wiengarten et al., 2013, in the following abbreviated as W13) I extended the work done in my Master’s thesis (Wiengarten, 2011). A first step was the numerical reproduction of known analytic formulations for the HMF by Parker and Weber & Davis, for which the Cronos code had to be extended to account for solar rotation effects (Section 2 in W13), which was realized in two frames of reference: co-rotating with the Sun, and an inertial observer frame. The Weber & Davis test case is illustrated in Section 3 of W13, and a more thorough description 29 CHAPTER 2. MHD SIMULATION OF THE INNER-HELIOSPHERIC MAGNETIC FIELD of its analytic formulation in its original form and the previously unpublished form in the co-rotating frame of reference are discussed in more detail in an addendum found in Section 2.4.1. For a more realistic modeling of the solar wind as compared to these analytic models, the computationally easier way of a potential field description of the coronal magnetic field was pursued, the results of which serve as input to my MHD simulations of the inner heliosphere. The potential field modeling was performed by our collaborators from the MPS, and the respective work is described in Jiang et al. (2010) and Section 4.1 of W13. The common approach in potential field models is to use magnetograms of the observed photospheric magnetic field as inner boundary conditions. This has the disadvantage that the solar disk cannot be completely observed from Earth at any given time. Instead, data is accumulated during one solar rotation, so that changes especially at the far side of the Sun during a rotation remain neglected. This can be circumvented by surface flux transport (SFT) models, which take into account the photospheric motions and evolution of magnetic flux. In Jiang et al. (2010) such a model was fed with sunspot group data, which are also available for longer times in the past in contrast to magnetograms, which allows for estimates of the heliospheric open magnetic flux before the space age. The resulting modeled ”magnetogram” was then used in a potential field approach (current sheet source surface (CSSS) model) yielding the radial magnetic field at a source surface at 10 solar radii, and were provided for typical solar minimum and maximum conditions for utilization as inner boundary conditions in my MHD simulations. The code also requires boundary conditions for the remaining plasma quantities, for which I used a set of empirical formulas as in Detman et al. (2006, 2011) based on the WSA approach. The radial velocity formula relies on the finding by Wang and Sheeley (1990) that it is inversely proportional to the flux-tube expansion factor (EF), which relates field strengths at the source surface with the one at the photosphere as connected by a field line. A large expansion is found for flux tubes originating at coronal hole boundaries (resulting in slow solar wind speed), while small expansion is found within coronal holes. Therefore, another quantity that can be used to set the resulting solar wind speed is the so-called footpoint distance (FPD) to the nearest coronal hole boundary. Both EF and FPD were found to give better results if used in conjunction than used alone (Arge et al., 2003). The required quantities needed to compute them, which are the footpoint locations of open field lines in the photosphere and their mapping via these field lines to the source surface, were also provided along with the radial magnetic field. I developed an algorithm to identify the position of the coronal hole boundaries and to compute the DFP, so that all the inner boundary conditions could be set. Another complication, however, arises from the necessity to initialize the magnetic field in the whole computational domain in a divergence-free manner. This cannot be done with analytic prescriptions for this arbitrary read-in magnetic field, but instead another algorithm was constructed to achieve this (for details see Section 2.4.2). 30 2.1. OVERVIEW The obtained results at Earth’s orbit for the solar minimum case were crudely compared to spacecraft data from that time, proving the existence of the modeled high-speed streams in the ecliptic at that time. Meanwhile, the typical structure of high-speed polar flows and low speeds in the ecliptic was also recreated. Therefore, my results were considered to be suitable for further usage in outer heliospheric models that put the inner boundary at a few AU and investigate the interaction with the ISM. A more detailed analysis of the results and comparison with spacecraft data for several more sets of input data from our collaborators at the MPS revealed later on that these data used in this study were inaccurate (see Section 2.3). This forced me to model the input data (i.e. the coronal magnetic field) by myself, which was the starting point for the work described in Chapter 3. 31 2.2 Wiengarten et al. (2013) MHD simulation of the inner-heliospheric magnetic field Journal of Geophysical Research (Space Physics) 32 JOURNAL OF GEOPHYSICAL RESEARCH: SPACE PHYSICS, VOL. 118, 29–44, doi:10.1029/2012JA018089, 2013 MHD simulation of the inner-heliospheric magnetic ﬁeld T. Wiengarten,1 J. Kleimann,1 H. Fichtner,1 R. Cameron,2 J. Jiang,2 R. Kissmann,3 and K. Scherer1 Received 3 July 2012; revised 14 November 2012; accepted 21 November 2012; published 27 January 2013. [1] Maps of the radial magnetic ﬁeld at a heliocentric distance of 10 solar radii are used as boundary conditions in the MHD code CRONOS to simulate a three-dimensional inner-heliospheric solar wind emanating from the rotating Sun out to 1 AU. The input data for the magnetic ﬁeld are the result of solar surface ﬂux transport modeling using observational data of sunspot groups coupled with a current-sheet source surface model. Among several advancements, this allows for higher angular resolution than that of comparable observational data from synoptic magnetograms. The required initial conditions for the other MHD quantities are obtained following an empirical approach using an inverse relation between ﬂux tube expansion and radial solar wind speed. The computations are performed for representative solar minimum and maximum conditions, and the corresponding state of the solar wind up to the Earth’s orbit is obtained. After a successful comparison of the latter with observational data, they can be used to drive outer-heliospheric models. Citation: Wiengarten, T., J. Kleimann, H. Fichtner, R. Cameron, J. Jiang, R. Kissmann, and K. Scherer (2013), MHD simulation of the inner-heliospheric magnetic field, J. Geophys. Res. Space Physics, 118, 29–44, doi: 10.1029/2012JA018089. and its deviations from the Parker ﬁeld, and numerical simulations [Lionello et al., 2006] have revealed that the original choice of parameters by Fisk [1996] represented an extreme, which has recently been conﬁrmed via a study of the consequences of the Fisk-ﬁeld for cosmic ray modulation [Sternal et al., 2011]. [4] Several authors have suggested modiﬁcations of both ﬁelds [Jokipii and Kota, 1989; Smith and Bieber, 1991; Schwadron, 2002; Burger and Hitge, 2004; Burger et al., 2008; Schwadron et al., 2008], which have been quantitatively compared in a study by Scherer et al. [2010]. Although very useful for many purposes, these analytical representations of the HMF cannot be used to reproduce actual measurements in sufﬁcient detail at any given time in the solar activity cycle. [5] For solar activity minima they are often too crude to catch small-scale HMF structures and for the HMF during maximum solar activity none of these analytical representations is valid. The ﬁeld expressions suggested by Zurbuchen et al. [2004] for use during solar maximum are a simple modiﬁcation of a representation for solar minimum [Zurbuchen et al., 1997] obtained by assuming a rather unlikely so-called Fisk angle [see Sternal et al., 2011]. [6] While, physically, the HMF originates in the Sun, it is conceptually customary to distinguish for modeling purposes between the coronal magnetic ﬁeld—ﬁlling the region from the solar surface out to a spherical “heliobase” [Zhao and Hoeksema, 2010] at several (tens of) solar radii—and the HMF beyond. This concept is used in many modeling attempts, which can be divided in two groups, potential ﬁeld reconstructions and MHD models [see Riley et al., 2006]. The latter approach is computationally more expensive but can account for more physics, direct time dependence and self-consistency. Examples for such 1. Introduction [2] Since Parker’s explanation of the expansion of the solar atmosphere as the solar wind [Parker, 1958], its basic, long-term average features have been validated by direct in situ as well as indirect remote measurements. At a radial distance of a few solar radii, i.e., at the so-called source surface, the heliospheric magnetic ﬁeld (HMF) is essentially directed radially and “frozen in” into the radially expanding solar wind plasma. Because the magnetic ﬁeld remains, due to an anchoring of the ﬁeld line footpoints in the solar atmosphere, inﬂuenced by the solar rotation, the ﬁeld lines form three-dimensional Archimedean or so-called Parker spirals [Parker, 1958] if, in a ﬁrst approximation, the Sun is considered as a rigid rotator. [3] For shorter periods, a more complex magnetic ﬁeld structure was suggested by Fisk [1996], who took into account both the motion of magnetic footpoints on the solar surface and the existence of a nonvanishing latitudinal ﬁeld component. He was able to derive analytical expressions that generalize the Parker ﬁeld correspondingly. Subsequent analytical investigations [Zurbuchen et al., 1997; Kobylinski, 2001; Schwadron, 2002; Schwadron et al., 2008] have quantiﬁed the complex structure of this so-called Fisk-ﬁeld 1 Institut für Theoretische Physik IV, Ruhr-Universität Bochum, Bochum, Germany. 2 Max-Planck-Institut für Sonnensystemforschung, Katlenburg-Lindau, Germany. 3 Institut für Astro-und Teilchenphysik, Universität Innsbruck, Innsbruck, Austria. Corresponding author: T. Wiengarten, Institut für Theoretische Physik IV, Ruhr-Universität Bochum, 44780 Bochum, Germany. ([email protected]) ©2012. American Geophysical Union. All Rights Reserved. 2169-9380/13/2012JA018089 29 CHAPTER 2. MHD SIMULATION OF THE INNER-HELIOSPHERIC MAGNETIC FIELD 2.3 Errata and further development Some minor mistakes were introduced unnoticed by the publisher during the production process: h • Equation (26) should read Bϕ = Br v vr ϕ i . • On page 36, paragraph [38]: reference to the boundary conditions for number density should be ”Panel (c) of Figures 5 and 6”. As mentioned in the introduction of this chapter, a more detailed subsequent evaluation of the results presented in W13 by comparing with multi-spacecraft observations revealed that they are inaccurate. The first possibility I considered was that the empirical interface used to derive the boundary conditions – in particular for the solar wind speeds – would need additional tuning since the approach of Detman et al. (2011) uses a different potential field model and source surface location. Although such a tuning would probably indeed have been necessary, a more severe problem turned out to be the location of the current-sheet in the radial magnetic field maps, because these maps were used as direct input in the boundary conditions and could therefore not be tuned. Although there are always slight differences in the location of the current-sheet depending on the magnetogram source and the potential field model used, the differences found for the input data used thus far were rather severe, as demonstrated in Figure 2.3. The Figure shows the resulting radial magnetic field from the SFT/CSSS model of Jiang et al. (2010) at ten solar radii around 2007.5 (a), 2007.87 (c), 2008.24 (e) in the left panels compared to GONG maps1 of footpoints of open field lines (red/green dots), as well as projections of the highest closed coronal field lines (blue), and the indicated neutral line position (black) at the respective times (right panels). Although these figures are in a different format (Hammer projections on the left) and show different quantities, they do allow for a comparison of the neutral line/current-sheet position. While the tilts in the current-sheet in the left panels vary only slightly over the course of more than half a year, this is not the case in the right panels, and there is only a resemblance found between the bottom two panels. A discussion of this finding with our colleagues from the MPS yielded that their model did not properly take into account newly emerging active regions near the equator. This issue was addressed in Jiang et al. (2014), but for my work I chose to take a different approach as described in the next chapter. 1 obtainable from http://gong.nso.edu/ 34 2.3. ERRATA AND FURTHER DEVELOPMENT (a) (b) (c) (d) (e) (f) Figure 2.3: Resulting radial magnetic field strength from the SFT/CSSS model of Jiang et al. (2010) at ten solar radii around 2007.5 (a), 2007.87 (c), 2008.24 (e) on the left. In the right panels GONG maps of footpoints of open field lines (red/green dots), as well as projections of the highest closed coronal field lines (blue), and the indicated neutral line position (black) at the respective times is shown for comparison. 35 CHAPTER 2. MHD SIMULATION OF THE INNER-HELIOSPHERIC MAGNETIC FIELD 2.4 2.4.1 Addenda Weber-Davis model in the inertial and co-rotating reference frame Parker’s assumption of magnetic field lines following the plasma flow is not valid close to the Sun since there the kinetic energy density is smaller than the magnetic energy density. Consequently, one can assume that the plasma flow follows the magnetic field lines well inside the Alfvénic radius (r rA ). This co-rotation should weaken when approaching rA , and for the limiting case r rA Parker’s model should be a good approximation. This was calculated by Weber and Davis (1967, hereafter WD) who assumed an azimuthal component of the flow velocity and solved a steady-state model with axial symmetry (no ϕ-dependence) in the equatorial plane. An extension for arbitrary latitudes was given by Sakurai (1985), but as it is solved numerically a comparison would go beyond the scope of a test case as intended here (but see, e.g., Preusse et al., 2005). Therefore the test case is restricted to the equatorial plane. In the WD model the components of the magnetic field and the flow velocity have azimuthal components so that B = Br (r)er + Bϕ (r)eϕ v = vr (r)er + vϕ (r)eϕ . The solenoidality of B (1.7) gives (R = b solar radius) 2 R r2 (2.1) r2 ρvr = const. (2.2) Br (r) = B0 and the equation of continuity (1.1) yields In the considered steady state the induction equation (1.4) reads ∂r (r(vr Bϕ − vϕ Br )) = 0 ⇒ r(vr Bϕ − vϕ Br ) = const. (2.3) WD determine this constant by noting that in a perfectly conducting fluid v is parallel to B in a frame that rotates with the Sun, yielding 2 vr Bϕ − vϕ Br = −ΩW D R B0 /r (2.4) where ΩW D is the angular velocity of the roots of the lines of force in the Sun (a detailed derivation is given in MacGregor and Pizzo (1983)). In their paper WD omitted any subscript on ΩW D so that one is tempted to consider it to be the photospheric solar rotation frequency. This in turn makes the original results of WD 36 2.4. ADDENDA seemingly unphysical and, consequently, Barker and Marlborough (1982, hereafter BM) took a different approach to determine the constant in Equation (2.3) by taking values at the solar surface Bϕ (R ) = B0,ϕ , Br (R ) = B0,r , vϕ (R ) = v0,ϕ = ΩR and vr (R ) = v0,r , yielding 2 B0 /r · (1 + f ) vr Bϕ − vϕ Br = −ΩR where f= v0,r |B0,ϕ | v0,ϕ B0,r (2.5) (2.6) and Ω (without subscript) denotes the solar rotation frequency. Later, MacGregor and Pizzo (1983) compared the approaches of WD and BM concluding that they just differ by the definition of the angular velocity Ω (as is obvious here from comparing Equations (2.4) and (2.5)), yielding ΩW D = (1 + f )Ω (2.7) which does not change the involved physics. Nevertheless, for comparing WD with the numerical results it is better to deal with the known value of the photospheric solar rotation frequency as was used by BM instead of the poorly known value for the angular velocity of the roots of the lines of force in the Sun, which are located within the Sun. Additionally, the values at the solar radius used by BM are also easily accessible for the numerical results. The following steps are again taken from WD who continue to solve (2.4) for Bϕ : Bϕ (r) = vϕ (r) − rΩW D Br (r) . vr (r) (2.8) To derive the azimuthal component of the flow velocity, the ϕ-component of the equation of motion (1.2) ρ vr Br ∂r (rvϕ ) = ∂r (rBϕ ) r µ0 r is considered, and after multiplication by r3 Br r2 2 ρvr r ∂r (rvϕ ) = ∂r (rBϕ ) µ0 (2.9) (2.10) is obtained, which allows to identify the terms in square brackets as constants (see Equations (2.1) and (2.2)). Integration yields rvϕ − Br rBϕ = const. ≡ L , µ0 ρvr (2.11) which corresponds to the total angular momentum carried away from the Sun per unit mass loss. Furthermore, the radial Alfvén Mach number vr MA ≡ (2.12) Br /(µ0 ρ)1/2 37 CHAPTER 2. MHD SIMULATION OF THE INNER-HELIOSPHERIC MAGNETIC FIELD can be introduced so that substituting Equations (2.8) and (2.12) into Equation (2.11) and solving for vϕ gives vϕ (r) = ΩW D r MA2 L/(ΩW D r2 ) − 1 . MA2 − 1 (2.13) Here, MA is a steadily rising function of r, reaching unity at the Alfvén critical radius rA , at which the denominator in Equation (2.13) vanishes. In order to keep 2 vϕ finite at rA , the numerator must vanish as well, so that L = ΩW D rA . From Equations (2.1), (2.2), and (2.12) it can be deduced that MA2 1 = 2 vr r 2 vA (rA )rA (2.14) is a constant (evaluated at rA ), finally allowing to find vϕ (r) = ΩW D r 1 − vr /vA (rA ) 1 − MA2 (2.15) and when substituted into Equation (2.8) ΩW D r 1 − (r/rA )2 Bϕ (r) = −Br vA (rA ) 1 − MA2 (2.16) from which the BM solution can easily be obtained by inserting Equation (2.7). The asymptotic behavior is qualitatively derived as follows: the radial velocity becomes almost a constant for large r, and with the equation of continuity follows MA ∝ r so that Bϕ , as well as vϕ , vary as 1/r according to Parker’s prediction. Following BM’s approach (i.e. using initial values at the solar surface), the above derivation can be adopted for the co-rotating frame of reference as well. The initial value for the azimuthal velocity vϕ (R ) = 0 has to be used here, and the equation of motion is extended by the fictitious force term (in the ecliptic plane), such that ρ Br vr ∂r (rvϕ ) = ∂r (rBϕ ) + 2ρΩvr r2 r µ0 r (2.17) instead of Equation (2.9). This, after following the above steps, eventually yields 2 1 R B0,ϕ rA 2 Br (r) 2 vϕ (r) = 2 v0,r 1 − MA − MA Ωr 1 − 2 (2.18) MA − 1 r Br (r) Br (rA ) r and (R /r)v0,r B0,ϕ + vϕ (r)Br (r) . (2.19) vr (r) Here, for the latter vϕ has to be substituted, which was not done here to avoid lengthy expressions. Bϕ (r) = The final equations for the inertial (2.15 & 2.16) as well as for the co-rotating frame of reference (2.18 & 2.19) are used to validate the numerical realizations of these reference frames in the preceding publication. 38 2.4. ADDENDA 2.4.2 Divergence-free initialization for arbitrary input data It is crucial to initialize the magnetic field in a divergence-free manner in Cronos with its constrained transport scheme, because it maintains the respective initial values of ∇ · B. For simple HMFs such as a Parker spiral with a flat current-sheet, this can be easily achieved, as such fields can also be described in terms of vector potentials, whose curl automatically fulfills the solenoidality condition (see Equation (1.10)). For arbitrary inner boundary conditions for the radial magnetic field as used here, additional care must be taken as there are no analytic descriptions available, and the situation is further complicated by the staggered-grid layout (see Section 1.2.2) employed in Cronos, where the magnetic field components are not cell-centered but centered on respective cell-faces. Starting from an arbitrary radial magnetic field distribution at the previous radial cell-face of the inner radial grid boundary Br (r0 − ∆r/2, ϑ, ϕ), an estimate for the azimuthal magnetic field component of the first cell layer is calculated via Bϕ (r0 , ϑ, ϕ ± ∆ϕ/2) = Bravg avg v . vravg ϕ (2.20) Due to the offset locations of these quantities (see Figure 1.11), respective averaging is required at the location of Bϕ . While vr,ϕ are known throughout the computational domain from respective simple radial initializations (see Section 4.3 in W13), which do not have to fulfill such severe restrictions as the magnetic field, the next cell’s radial magnetic field is actually not known at this point. This is estimated initially by a r−2 dependence, which can be augmented by repeated execution of the following algorithm. A further necessary, but valid approximation that has to be made is that the ϑ-components vanish, which is in agreement with open and radial field lines resulting at the outer boundary of potential field models of the solar corona. The idea for a divergence-free radial initialization is to go through all cell layers in order of increasing radial position and adapting Br such that solenoidality holds in the subsequent cells. This can be achieved by solving the discretized solenoidality condition for the radial magnetic field Br (ri+1/2,j,k ) in the next cell (notation according to Equation (1.9)), while the previous one Br (ri−1/2,j,k ), as well as the azimuthal magnetic fields Bϕ (ri,j,k−1/2 ) and Bϕ (ri,j,k+1/2 ) at the respective cell boundaries are known at this point. It has also to be considered that the divergence itself is computed at the cell center at ri,j,k . Starting with the solenoidality condition, where Bϑ is already set to zero, 0=∇·B= 1 1 2 ∂ (r B ) + ∂ϕ (Bϕ ) , r r r2 r sin(ϑ) 39 (2.21) CHAPTER 2. MHD SIMULATION OF THE INNER-HELIOSPHERIC MAGNETIC FIELD and discretizing with finite differences yields 2 2 1 ri+1/2,j,k Br (ri+1/2,j,k ) − ri−1/2,j,k Br (ri−1/2,j,k ) 0= 2 ri,j,k ∆r Bϕ (ri,j,k+1/2 ) − Bϕ (ri,j,k−1/2 ) 1 + , kri,j,k k sin(ϑj ) ∆ϕ (2.22) so that solving for Br (ri+1/2,j,k ) finally gives kri−1/2,j,k k 2 Br (ri+1/2,j,k ) = Br (ri−1/2,j,k ) kri+1/2,j,k k kri,j,k k∆r(Bϕ (ri,j,k+1/2 ) − Bϕ (ri,j,k−1/2 )) . + r2i−1/2,j,k sin(ϑj )∆ϕ (2.23) The steps described above have to be performed for each cell layer. This causes difficulties when running on multiple CPUs: in order to compute the next cell layer, one needs to know the quantities from the previous one, but this is not necessarily the case here as usually each processor has only access to its respective sub-grid. The solution is to devise a global array and have each CPU calculate the initialization on the whole grid and then to assign only the respective sub-grid values to the single CPUs. It has been made sure that the solenoidality condition is satisfied to the degree of machine precision. The following code fragment depicts the actual implementation in the User class: void OBS::div-free-init (Data &gdata, int ibeg[3], int iend[3]) { /* The globally declared array containing the initial fields looks like REAL global[3][3][rad_cells+7][190][370]; with dimensional sizes due to: 2 for sets in time between which to average +1 for initial divergence-free-init, 3 for vrad,brad,bphi , rest due to resolution */ //apply div-free-init to time-interpolated global array if using time int p;//p corresponds to first dimension in global array if (time_usage == 1){p=2;} else {p=0;} //some declarations REAL B_r_neu = 0; REAL delta_Bphi =0; REAL v_phi = 0; REAL v_r = 0; REAL B_r = 0; REAL B_phi = 0; // apply solenoidality condition while initializing Bphi: // go through whole grid, starting at inner radial grid boundary for (int i = -1; i < gdata.global_mx[0]+3; i++){ // set Bphi first: 40 2.4. ADDENDA for (int j = ibeg[1]; j <= gdata.global_mx[1]+3; j++){ for (int k = ibeg[2]; k <= gdata.global_mx[2]+3; k++){ REAL xx = gdata.get_x_global(i,0); //radial coordinate REAL theta = gdata.get_y_global(j,0); int i_min = i - 1; //previous radial index //assume v_phi=0 in inertial frame and transform to co-rotating one v_phi = -w*xx*sin(theta); //averaged quantities at location of Bphi are required to compute it if (k==gdata.global_mx[2]+3){ //periodic treatment at outer phi boundary v_r = (global[p][0][i+3][j+3][k+3] + global[p][0][i+3][j+3][ibeg[2]+7+3])/2; B_r = (global[p][1][i+3][j+3][k+3] + global[p][1][i_min+3][j+3][k+3] + global[p][1][i+3][j+3][ibeg[2]+7+3] + global[p][1][i_min+3][j+3][ibeg[2]+7+3])/4; } else{ v_r = (global[p][0][i+3][j+3][k+3] + global[p][0][i+3][j+3][k+1+3])/2; B_r = (global[p][1][i+3][j+3][k+3] + global[p][1][i_min+3][j+3][k+3] + global[p][1][i+3][j+3][k+1+3] + global[p][1][i_min+3][j+3][k+1+3])/4; } //B_phi follows from the vanishing electric field (E = v x B) in the co-rotating frame //and neglected theta-components B_phi = B_r/v_r * v_phi; global[p][2][i+3][j+3][k+3] = B_phi; //store bphi in global array } } // set Br accordingly: for (int j = ibeg[1]; j <= gdata.global_mx[1]+3; j++){ for (int k = ibeg[2]; k <= gdata.global_mx[2]+3; k++){ REAL r_i = gdata.get_x_global(i,1); //radial distance at current Br-position REAL r_ip = gdata.get_x_global(i-1,1); //radial distance of previous Br-position REAL B_rp = global[p][1][i-1+3][j+3][k+3]; //value of previous Br REAL dr = gdata.getCen_dx_global(0,i); //radial cell size REAL dphi = gdata.getCen_dx_global(2,k);//azimuthal cell size REAL theta = gdata.get_y_global(j,0); //theta coordinate if (k==ibeg[2]) //periodic treatment at lower azimuthal boundary delta_Bphi = -global[p][2][i+3][j+3][gdata.global_mx[2]+3-7+3] + global[p][2][i+3][j+3][k+3]; else delta_Bphi = -global[p][2][i+3][j+3][k-1+3] + global[p][2][i+3][j+3][k+3]; //difference in Bphi at current cell REAL r_phi = gdata.get_x_global(i,0); //radial distance at cell center B_r_neu = power(r_ip/r_i,2)*B_rp - power(r_phi/r_i,2)*dr*delta_Bphi/(r_phi*sin(theta)*dphi); //solenoidality condition solved for next Br global[p][1][i+3][j+3][k+3] = B_r_neu; // store in global array } } } } 41 CHAPTER 2. MHD SIMULATION OF THE INNER-HELIOSPHERIC MAGNETIC FIELD 42 Chapter 3 Modeling background solar wind 3.1 Overview The work described in this chapter based on the publication Cosmic Ray transport in heliospheric magnetic structures. I. Modeling background solar wind using the Cronos magnetohydrodynamic code (Wiengarten et al., 2014, W14) was done in the framework of a joint project with the Universities of Kiel and Potchefstroom, South Africa. The project focused on the influence of CIRs on the propagation of energetic particles, where the respective CIR disturbed background solar wind modeled in my MHD simulations is used as input to a respective SDE model for the transport of energetic particles in the heliosphere. The idea behind this work is that the transport of charged energetic particles such as GCRs or Jovian electrons is largely governed by the properties of the plasma that these particles traverse. Due to the large amount of simultaneous in-situ measurements of the heliospheric plasma environment as well as of energetic particles, the heliosphere can be considered as a natural laboratory to investigate these transport processes. Meanwhile, the transport in the Galaxy is also important, but modeling suffers from the very small number of measurements there and at the sites of acceleration to high energies such as supernovae remnants. Of particular importance for the heliospheric modulation is the influence of transient structures such as CIRs and CMEs (with a focus on CIRs in this case) on the energetic particle propagation, because these structures are associated with magnetic field and turbulence enhancements, as well as shocks that can serve as diffusion barriers (see Figure 3.1), but also as accelerators of low-energy particles. A fully three-dimensional and possibly also time-dependent description of the heliospheric plasma is necessary for an adequate energetic particle modeling. Although simple analytical prescriptions such as the Parker (1958) model (undisturbed solar wind) or the Giacalone et al. (2002) model (CIR-disturbed wind) are useful starting points for such investigations, a much more coherent approach is to use MHD models that can self-consistently account for stream interactions and the formation of shocks, and which can be driven by observations-based input data. Therefore, the work described in the previous 43 CHAPTER 3. MODELING BACKGROUND SOLAR WIND Figure 3.1: Pseudo-particle trajectories from a feasibility study of the SDE code using a toy-model heliosphere with an embedded CIR structure resulting from an MHD simulation. Particles follow the magnetic field spiral pattern, but hardly diffuse across the CIR structure.[Courtesy: A. Kopp] chapter is principally suited, but as mentioned before, the input data I used thus far were found to be incorrect to some extent, so that a new approach had to be found as described below. A first mandatory task is the demonstration of the MHD code’s ability to correctly model CIR structures. This is the topic of Section 3 of W14, where analytic inner boundary conditions of alternating fast and slow streams were prescribed, and the resulting CIR configuration – including pressure enhancements, stream interfaces between fast and slow wind, as well as waves steepening into shocks traveling away from the interface – was successfully compared to previous work by Pizzo (1982). The further evolution of interacting – initially distinct – CIRs to eventually form merged interaction regions with increasing heliocentric distance was also analyzed. This successful validation allowed for an investigation of these structures in a re44 3.1. OVERVIEW alistic solar wind, i.e. in one that is based on observational data. To achieve this, a more suitable model was required for the coronal magnetic field that drives the WSA model to derive inner boundary conditions for the MHD simulations as described in Section 4 of W14. Since the empirical formulas in the WSA approach are tuned to be applied to the coupled PFSS (Altschuler and Newkirk, 1969) and SCS (Schatten, 1971) models in conjunction with GONG magnetograms, this road was followed by implementing these models in Matlab. The PFSS solution describes the coronal magnetic field out to the source surface at 2.5R , where the field is forced to become open and radial. A drawback is that the resulting radial magnetic field distribution at the source surface does not exhibit a thin current sheet as opposed to observations. This is addressed by subsequently involving the SCS model in the region from the source surface out to the inner boundary location of the MHD simulations at 0.1 AU. Both PFSS and SCS usually solve the Laplace equation for a scalar potential ∆Ψ = 0 for the coronal magnetic field by means of spherical harmonics. Such an approach was initially taken by myself as well, but for grid resolutions approaching the magnetograms’ one is faced with numerical artifacts in the reconstructed magnetograms, which is the so-called ringing effect related to Gibb’s phenomenon occurring at sharp or discontinuous reconstructions of periodic signals with Fourier series. This can be avoided by solving the Laplace equation with a finite differences approach, and therefore I used a respective public code (FDIPS1 ) for the PFSS model henceforth. A respective coupling of such a finite-difference approach to the coronal potential field with the WSA model as used in this work had not been undertaken previously. The set of empirical formulas based on the WSA approach was slightly adapted in this work, but is otherwise the same as in W13. However, particularly for the resulting solar wind speeds, which rely on the topology of the coronal potential field via the flux tube expansion factor (EF) and the distance to the nearest coronal hole boundary (FPD), it is necessary to perform field line tracings (FLT) in the solutions for the coronal magnetic field. Previously, these data were provided by the collaborators at the MPS, so that I had to perform the FLT by myself now. The basic idea in FLT for a given vector field is simply to take a small step tangentially to the field at a predefined starting point and repeat it from the newly reached point until a predefined boundary is reached. This can be done numerically by simple forward integration schemes, however, the tricky part is the size of the small steps. I implemented a method of taking adaptive step sizes, namely the embedded Runge-Kutta 4/5 algorithm (e.g. as described in Press et al. (2007)), which takes adequately large or small steps depending on the fields topology at the current point. The algorithm in application to the field line tracing in the coronal potential field was validated by comparing with respective data from the GONG web page, and consequently – after computing the EF and FPD as in W13 – the complete set of inner boundary conditions for all MHD quantities could be set. 1 http://csem.engin.umich.edu/tools/FDIPS/ 45 CHAPTER 3. MODELING BACKGROUND SOLAR WIND While this procedure was principally applicable to any period of time for which GONG magnetograms are available (back to 2006), a time period of descending solar activity (August 2007) was chosen, because of recurrent and stable CIR structures observed at that time. Furthermore, a unique spacecraft constellation was present then with Ulysses crossing the equatorial plane during its third fast-latitude scan and the recently launched Stereo twin spacecrafts leading and trailing Earth (see Section 1 of W14). MHD simulations during that time (Carrington Rotations 2059-2061) were performed on new computer clusters at our institute comprising of 64 and 48 cores, which put typical run-times at several hours to a week (depending on resolution). The results (Section 5 of W14) were compared to respective spacecraft data provided by our collaborators of the University of Kiel. This required me to perform an interpolation of my results on the numerical grid along the spacecraft trajectories. In order to obtain results as close as possible to the reference data, the tuning parameters in the empirical formulas for the inner boundary conditions were adapted, and the final results (Figures 9–11 in W14) show usually good agreement. Particularly the timing and strength of high-speed streams and resulting CIR structures could be reproduced. Nevertheless, some discrepancies remain and reasons for this are discussed in the conclusions (Section 6 of W15). For the comparison with spacecraft data, the simulation box was restricted by an outer boundary at 2 AU to save computing time while still including the trajectories of the considered spacecraft with Ulysses at approximately 1.3 AU at that time. For further utilization in the SDE modeling, the simulation box has to be extended to at least include Jupiter (at about 5 AU), but ideally to much larger distances. Here, additional simulations have been performed out to 10 AU and the evolution of the CIR structures with heliocentric distance was analyzed in a similar fashion as for the validation case described above. 46 3.2 Wiengarten et al. (2014) Cosmic Ray transport in heliospheric magnetic structures. I. Modeling background solar wind using the Cronos magnetohydrodynamic code The Astrophysical Journal 47 The Astrophysical Journal, 788:80 (16pp), 2014 June 10 C 2014. doi:10.1088/0004-637X/788/1/80 The American Astronomical Society. All rights reserved. Printed in the U.S.A. COSMIC RAY TRANSPORT IN HELIOSPHERIC MAGNETIC STRUCTURES. I. MODELING BACKGROUND SOLAR WIND USING THE CRONOS MAGNETOHYDRODYNAMIC CODE T. Wiengarten1 , J. Kleimann1 , H. Fichtner1 , P. Kühl2 , A. Kopp2 , B. Heber2 , and R. Kissmann3 2 1 Institut für Theoretische Physik IV, Ruhr-Universität Bochum, Germany Institut für Experimentelle und Angewandte Physik, Christian-Albrecht-Universität zu Kiel, Germany 3 Institut für Astro- und Teilchenphysik, Universität Innsbruck, Austria Received 2014 February 20; accepted 2014 April 29; published 2014 May 23 ABSTRACT The transport of energetic particles such as cosmic rays is governed by the properties of the plasma being traversed. While these properties are rather poorly known for galactic and interstellar plasmas due to the lack of in situ measurements, the heliospheric plasma environment has been probed by spacecraft for decades and provides a unique opportunity for testing transport theories. Of particular interest for the three-dimensional (3D) heliospheric transport of energetic particles are structures such as corotating interaction regions, which, due to strongly enhanced magnetic field strengths, turbulence, and associated shocks, can act as diffusion barriers on the one hand, but also as accelerators of low energy CRs on the other hand as well. In a two-fold series of papers, we investigate these effects by modeling inner-heliospheric solar wind conditions with a numerical magnetohydrodynamic (MHD) setup (this paper), which will serve as an input to a transport code employing a stochastic differential equation approach (second paper). In this first paper, we present results from 3D MHD simulations with our code CRONOS: for validation purposes we use analytic boundary conditions and compare with similar work by Pizzo. For a more realistic modeling of solar wind conditions, boundary conditions derived from synoptic magnetograms via the Wang–Sheeley–Arge (WSA) model are utilized, where the potential field modeling is performed with a finite-difference approach in contrast to the traditional spherical harmonics expansion often utilized in the WSA model. Our results are validated by comparing with multi-spacecraft data for ecliptical (STEREO-A/B) and out-of-ecliptic (Ulysses) regions. Key words: magnetohydrodynamics (MHD) – methods: numerical – shock waves – solar wind – Sun: heliosphere – Sun: magnetic fields Online-only material: color figures structures are the corotating interaction regions (CIRs) that are formed during the interaction of fast solar wind streams from coronal holes with preceding slow solar wind and usually persist for several solar roations (e.g., Balogh et al. 1999, and references therein). These structures not only lead to particle acceleration, but also to the modulation of GCR and Jovian electron spectra (Richardson 2004). Indeed Ulysses measurements, as described in, e.g., Marsden (2001), at high heliolatitudes during the first orbit of the spacecraft around the Sun indicated that CIRs represent the major constituent for the three-dimensional heliospheric structure close to solar minimum (e.g., Heber et al. 1999). A major surprise came from the measurements of accelerated energetic particles and GCRs that showed clear periodic signals even above the poles of the Sun where no in situ signals of CIRs were registered by the plasma or magnetic field instrument (Kunow et al. 1995). Electron measurements, however, indicate no variation at these region. It can be speculated that the differences are caused by the fact that both GCRs and locally accelerated particles have an extended source, while MeV electrons in the inner heliosphere originate from a point-like source as mentioned above (Chenette 1980). Thus, there is a different influence of CIRs on GCR and Jovian electron flux variations as investigated recently by Kühl et al. (2013). If the electron source is a well-localized point-like source in the heliosphere—namely the Jovian magnetosphere—it is mandatory to describe the structure of the plasma stream in the whole inner heliosphere up to several AU and not only at the location of different spacecraft measuring these particles. 1. INTRODUCTION 1.1. Cosmic Ray Transport The understanding and appropriate modeling of the transport of charged energetic particles in turbulent magnetic fields remains one of the long-standing challenges in astrophysics and space physics. The many simultaneous in situ observations of both the heliospheric magnetic field and different energetic particle populations make the heliosphere a natural laboratory for corresponding studies. Of particular interest in this context are galactic cosmic rays (GCRs) and Jovian cosmic ray electrons. While the former traverse the whole three-dimensionally structured heliosphere and allow for studies of its large-scale variations as well as of the evolution of heliospheric turbulence (e.g., Zank et al. 1996; Heber et al. 2006; Potgieter 2013a; Bruno & Carbone 2013), the latter represent a point-like source and are, thus, well suited for analyses of anisotropic spatial diffusion (e.g., Ferreira et al. 2001a, 2001b; Sternal et al. 2011; Strauss et al. 2013). In order to exploit this natural laboratory in full, it is necessary to reproduce the measurements with simulations that do contain as much as possible of the three-dimensional structure of the plasma background within which the cosmic rays are propagating. Significant progress has been made regarding the modeling of the quiet solar wind (e.g., Potgieter 2013b), but much remains to be done to implement the many features that are structuring the solar wind and the heliospheric magnetic field, into transport models of cosmic rays. Particularly interesting 1 3.3. FURTHER DEVELOPMENT AND FUTURE DIRECTIONS 3.3 Further development and future directions The companion paper of W14 describing the utilization of my MHD simulations’ results in the SDE code is currently in preparation. For the respective modeling of the charged energetic particle propagation, it is desirable to have MHD results within the whole heliosphere and the surrounding disturbed ISM. This would require on the one hand an MHD model of the solar corona, and on the other hand the coupling of my setup with global heliospheric models, which are also performed in our group by J. Kleimann and K. Scherer (see Chapter 6). A coronal MHD model is a challenging task as mentioned before, and for the SDE modeling focusing on the influence of CIRs on GCRs and Jovian electrons, it is not vital to have such a model, as CIRs develop at much greater heliocentric distances only. However, for a study of SEPs such a model would be desirable. A coupling of coronal MHD simulations for rather simple configurations (but, e.g. including a CME) with SEP models is already going on in the framework of a more elaborate future project with W. Dröge at the University of Würzburg. Meanwhile, a coupling with global heliospheric models is principally possible, but on the one hand, the necessary computer resources become huge due to the long propagation times and different required scales, and on the other hand, the SDE modeling beyond the TS in the heliosheaths is more difficult due to the present sub-Alfvénic flows and unknown turbulence evolution there. Therefore, for the current SDE modeling, MHD simulations within the TS have been the next desired step. As the propagation time of typical solar wind parcels from the Sun to the termination shock is about a year, it is unreasonable to assume stationary inner boundary conditions due to the changes in the coronal field, and a time-dependent approach is necessary. This has been achieved by respective interpolation in time between incrementing data sets for the inner boundary values of adjacent Carrington rotations. To reduce computational costs, the grid has been coarsened with increasing radial distance, and the polar axis’ have been included with one large cell each to avoid small time steps (see Appendix A). Furthermore, it was found that the implementation of solar rotation via the fictitious forces as source terms becomes problematic for large heliocentric distances beyond about 10 AU. Instead, a respective conservative scheme was implemented as described in the next chapter (specifically Section 4.4). Due to the still enormous computational costs even on a new 256 core cluster at our institute, simulations were performed out to 20 AU only, and since there is no principle change in the resulting configuration it is not further presented here. Another way of improving the subsequent SDE modeling is to include a turbulence model to be solved self-consistently alongside the MHD equations, because in the SDE model transport coefficients are employed that actually depend not only on the large-scale MHD quantities, but particularly the diffusion coefficients depend 49 CHAPTER 3. MODELING BACKGROUND SOLAR WIND on the magnetic fluctuations. A respective implementation of such a turbulence transport model is the topic of the next two chapters. 50 Chapter 4 Implementing turbulence transport 4.1 Overview From the early measurements of the solar wind it was already apparent that there are random fluctuations super-imposed onto the large-scale magnetic field, solar wind velocity and density, so that it has to be considered as an inherently turbulent plasma (Coleman, 1968; Belcher and Davis, 1971). While astrophysical plasmas should be considered as turbulent media in general due to the very high Reynolds numbers, the numerous available spacecraft measurements made in the heliosphere make the latter a natural laboratory for respective investigations. It is assumed that sources for turbulence (see, e.g., the review by Zank, 1999) are located not only in the corona from where it is advected into interplanetary space, but that it is also generated in-situ by the turbulence itself at density gradients throughout the heliosphere. In the inner heliosphere (within about 10 AU) the main drivers for turbulence are shear at differently fast solar wind stream interfaces, as well as transients such as CIRs and CMEs with their associated shocks. These, however, weaken with increasing heliocentric distance as these structures become merged interaction regions. In the outer heliosphere isotropization of newly born pickup ions is the main agent for generating fluctuations, during which the initially ring-distributed pickup ions isotropize and generate Alfvén waves, which cascade to ever smaller scales and eventually dissipate and heat the solar wind as measured by the Voyager 2 spacecraft (see Figure 4.1). Thus, the temperature does not behave as expected for an adiabatically expanding and cooling gas, but there is actually even an increase. Amongst the further important effects of turbulence (see, e.g., the review by Miesch et al. (2015)) are its catalytic impact on reconnection (e.g. Yokoi et al., 2013), which is one explanation for the heating of the solar corona, as well as its relevance for the propagation of energetic particles. Although particles mainly follow the mean magnetic field, the magnetic field fluctuations, which are usually perpendicular to 51 CHAPTER 4. IMPLEMENTING TURBULENCE TRANSPORT Figure 4.1: Temperature evolution in the heliosphere as measured by the Voyager 2 spacecraft (black solid line) compared to an adiabatic model (blue dashed). The Voyager measurements were taken from the omniweb interface [http://omniweb.gsfc.nasa.gov/], and the adiabatic model is T = T0 (r/r0 )−4/3 . the mean one, lead to a perpendicular diffusion of energetic particles towards neighboring field lines. In the first attempts to involve these effects in models for the transport of energetic particles in the heliosphere, diffusion parameters have been approximated by expressing the fluctuations in terms of the large-scale fields (see the review of Potgieter, 2013). Some recent so-called (ab-initio models (e.g. Engelbrecht and Burger, 2013)), however, formulate the diffusion coefficients explicitly based on the small-scale fluctuations, on which new insights have been gained both from improved analysis of solar wind measurements (e.g. Horbury and Osman, 2008; Horbury et al., 2012) on the one hand, but particularly by new modeling approaches on the other hand (e.g., Zank et al., 1996; Breech et al., 2008; Oughton et al., 2011; Zank et al., 2012). These latter models describe the evolution of turbulent energy density, cross-helicity, and correlation lengths for low-frequency turbulence, as observed in the solar wind, with increasing distance from the Sun. These integral quantities describe the large-scale behavior of turbulence and wave modes, which is a useful approach to circumvent the difficulties that would arise for a complete description of the fluctuations, which is on much smaller scales than typical structures in the heliosphere. Particularly for numerical studies such an approach is desirable because of the reduced required computer resources. The following publication Implementing turbulence transport in the Cronos framework and application to the propagation of CMEs (Wiengarten et al., 2015, W15) describes my efforts to extend the Cronos setup with a respective turbulence trans52 4.1. OVERVIEW port model. Besides the general aim to model and understand the evolution of turbulence in the heliosphere, another aim is to utilize the results to augment the calculation of transport parameters in the SDE modeling. All the above-mentioned turbulence transport models have the drawback that they use prescribed solar wind conditions (usually constant speed and Parker field geometry), but it is desirable – particularly in the light of my previous work – to have a self-consistently computed background solar wind, i.e. a coupling of turbulence transport and MHD equations as described in Section 2 of W15. The basic idea is to decompose the magnetic and velocity fields in the MHD equations into respective mean and fluctuating components, which after applying appropriate averaging schemes leads to the so-called Reynolds-averaged MHD equations (see Usmanov et al. (2011) for a thorough description). In comparison with the usual MHD equations they also comprise of terms involving the fluctuations, which requires additional evolution equations for these quantities. Namely, these are the total turbulent (magnetic plus kinetic) energy density and the energy density difference, from both of which the particular magnetic and kinetic components can be computed. Furthermore, the so-called cross-helicity describes the ratio between inward and outward propagating modes. Turbulence cascades in the solar wind and is eventually dissipated into heat, where the rate of dissipation is controlled by the so-called correlation lengths, whose evolution is also described by additional equations in these models. While most modeling attempts concentrated on the resulting turbulence distribution in the heliosphere for idealized solar wind conditions, thereby neglecting a self-consistent coupling with the Reynolds-averaged MHD equations, such a coupled approach was first taken by Usmanov et al. (2011). To become familiar with this subject and test Cronos’ ability to be correspondingly extended (Section 2 of W15), I started with a feasibility study and compared with these authors’ results (Section 3 of W15). In due course and after discussions with A. Usmanov, some errors were found in their model, which I subsequently corrected in my paper. Another inherent drawback of the Usmanov model was its restriction of applicability to regions in the heliosphere, where the solar wind speed is much higher than the Alfvén speed. As in my previous work employing the WSA model this condition is not fulfilled towards the inner boundary, I sought to extend the model further to be applicable to these regions as well. A respective more general turbulence transport model was recently presented by Zank et al. (2012) with further refinements in Dosch et al. (2013) and Adhikari et al. (2015). This model has been implemented in a still simplified manner here, but see Chapter 5 for a full implementation of that model. As compared to the model of Usmanov et al. (2011) it allows for the inner boundary location to be located closer to the Sun, and it properly incorporates turbulence generation via shear (Section 4 of W15). However, the Zank model has thus far been solved for idealized solar wind conditions only. Therefore, a significant 53 CHAPTER 4. IMPLEMENTING TURBULENCE TRANSPORT improvement provided in W15 is the self-consistent coupling of this model with the Reynolds-averaged MHD equations. These extensions allowed for a subsequent novel application of the model to a CME propagation scenario (Section 5 of W15). The results revealed that, on the one hand, turbulence does not act back strongly on the large-scale quantities describing the CME, so that general studies of CME propagation need not be extended in this direction. On the other hand, CMEs have a strong effect on the turbulence quantities. These effects are as manifold as the CME structure is complex, with partially increased and decreased turbulence levels. A subsequent investigation of the effects of such complicated turbulence distributions on the energetic particle propagation will be quite interesting. Another extension applied to the code setup is a conservative treatment of the co-rotating reference frame as opposed to the formerly applied realization with the required fictitious forces as source terms. This is described in more detail in an addendum in Section 4.4. 54 4.2 Wiengarten et al. (2015) Implementing turbulence transport in the Cronos framework and application to the propagation of CMEs The Astrophysical Journal 55 The Astrophysical Journal, 805:155 (15pp), 2015 June 1 doi:10.1088/0004-637X/805/2/155 © 2015. The American Astronomical Society. All rights reserved. IMPLEMENTINGTURBULENCE TRANSPORT IN THE CRONOS FRAMEWORK AND APPLICATION TO THE PROPAGATION OF CMEs T. Wiengarten1, H. Fichtner1, J. Kleimann1, and R. Kissmann2 1 Institut für Theoretische Physik IV, Ruhr-Universität Bochum, Germany; [email protected] 2 Institut für Astro- und Teilchenphysik, Universität Innsbruck, Austria Received 2015 February 2; accepted 2015 April 1; published 2015 May 28 ABSTRACT We present the implementation of turbulence transport equations in addition to the Reynolds-averaged magnetohydrodynamicequations within the Cronosframework. The model is validated by comparisons with earlier ﬁndings before it is extended to be applicable to regions in the solar wind that are not highly super-Alfvénic. We ﬁnd that the respective additional terms result in absolute normalized cross-helicity to decline more slowly, while a proper implementation of the mixing terms can even lead to increased cross-helicities in the inner heliosphere. The model extension allows us to place the inner boundary of the simulations closer to the Sun, where we choose its location at 0.1 AU for future application to the Wang–Sheeley–Arge model. Here, we concentrate on effects on the turbulence evolution for transient events by injecting a coronal mass ejection (CME). We ﬁnd that the steep gradients and shocks associated with these structures result in enhanced turbulence levels and reduced cross-helicity. Our results can now be used straightforwardly for studying the transport of charged energetic particles, where the elements of the diffusion tensor can now beneﬁt from the self-consistently computed solar wind turbulence. Furthermore, we ﬁnd that there is no strong back-reaction of the turbulence on the large-scale ﬂow so that CME studies concentrating on the latter need not be extended to include turbulence transport effects. Key words: magnetohydrodynamics (MHD) – methods: numerical – shock waves – solar wind – Sun: magnetic ﬁelds – turbulence Supporting material: animations After the pioneering papers by Tu et al. (1984) and Zhou & Matthaeus (1990a, 1990b, 1990c) and the ﬁrst systematic studies of the radial evolution of turbulence quantities (Zank et al. 1996), major progress has only been achieved in recent years. First, Breech et al. (2008) improved the previous modeling by considering off-ecliptic latitudes. Second, in Usmanov et al. (2011) this turbulence model was extended to full time dependence and three spatial dimensions. In the same paper the turbulence evolution equations were simultaneously solved along with a large-scale magnetohydrodynamic (MHD) model of the supersonic solar wind. While Usmanov et al. (2012, 2014) extended this study further with the incorporation of pickup ions and electrons as separate ﬂuids, as well as an eddy-viscosity approximation to self-consistently account for turbulence driven by shear, Oughton et al. (2011) undertook another extension of the modeling by considering not only one but two, mutually interacting, incompressible components, namely quasi-two-dimensional turbulent and wave-like ﬂuctuations. The effect on turbulence due to changing solar wind conditions in the outer heliosphere during the solar cycle have recently been addressed by Adhikari et al. (2014). All of the mentioned models have one limitation in common, namely the condition that the Alfvén speed VA is signiﬁcantly lower than that of the solar wind U, which precludes their application to the heliocentric distance range below about 0.3 AU. The turbulence in this region and its consequences for the dynamics of the solar wind have been studied on the basis of the somewhat simpliﬁed transport equations for the wave power spectrum and the wave pressure (e.g., Tu & Marsch 1995; Hu et al. 1999; Vainio et al. 2003; Shergelashvili & Fichtner 2011). The Alfvén speed limitation in the more detailed models was very recently removed with a new, comprehensive approach to the general problem of the 1. INTRODUCTION In a recent publication (Wiengarten et al. 2014), we presented results for inner-heliospheric solar wind conditions from simulations based on observational boundary conditions derived from the Wang–Sheeley–Arge (WSA; Arge & Pizzo 2000) model. The obtained three-dimensional (3D) conﬁgurations provide the basis for subsequent studies of energetic particle transport. For a long time, the standard approach to couple the solar wind dynamics to the transport of charged energetic particles has been to parameterize all transport processes in terms of “background” quantities such as the large-scale solar wind velocity and heliospheric magnetic ﬁeld strength. This was particularly true for the treatment of spatial diffusion (see, e.g., Potgieter 2013). Only in recent years were so-called ab initio models developed, in which the diffusion coefﬁcients are formulated explicitly as functions of smallscale, low-frequency turbulence quantities such as the magnetic ﬁeld ﬂuctuation amplitude or the correlation length (e.g., Parhi et al. 2003; Pei et al. 2010; Engelbrecht & Burger 2013). This signiﬁcant improvement has been possible, on the one hand due to the development of turbulence transport models that, most often, describe the evolution of turbulent energy density, cross-helicity, and correlation lengths for low-frequency turbulence, as observed in the solar wind, with increasing distance from the Sun;see, e.g., the review in Zank (2014). On the other hand, our present-day knowledge about turbulence in the solar wind (see, e.g., the reviews by Matthaeus & Velli 2011; Bruno & Carbone 2013) has increased tremendously since the ﬁrst measurements (Coleman 1968; Belcher & Davis 1971) thanks to highly sophisticated analyses (e.g., Horbury & Osman 2008; Horbury et al. 2012). 1 4.3. FURTHER DEVELOPMENT 4.3 Further development As described in more detail in the following chapter, there has been a development after discussions with S. Oughton suggesting that the choice of structural similarity parameters for 2D axisymmetric turbulence stated by Zank et al. (2012) might be incorrect, and should be a = b = 1/2 instead of a = 1/2, b = 0. Since the model presented here was derived in simplifying the respective equations with the latter choice, the severity of the error thus made is investigated and revised results are presented. However, it turned out that the topological changes are small and to some extent even negligible. Figure 4.2: Comparison of large-scale MHD quantities during a CME passage between the cases b = 0 (red lines) and b = 1/2 (black) in the same format as Figure 10 in W15, i.e. at two different co-latitudes of 75◦ (dashed) and 85◦ (solid). The revised set of equations is obtained in analogy with Appendix A and B of W15, but now for b = 1/2 and using Equation (5.10) as motivated in the next chapter. This gives an additional term (fourth term on right hand side) in the evolution 57 CHAPTER 4. IMPLEMENTING TURBULENCE TRANSPORT equation for the cross-helicity Z 2 σC ∂t (Z σC ) + ∇ · (UZ σC + VA Z ) = ∇ · U + 2VA · ∇Z 2 + Z 2 σD ∇ · VA 2 Z 2 σD αZ 3 f − (σC ) + 2b √ B̂ · (B̂ · ∇)B − ρ λ 2 2 2 + hz+ · S+ i − hz− · S− i , (4.1) while the remaining turbulence transport as well as large-scale MHD equations remain unchanged on keeping a = 1/2. Figure 4.3: Meridional contour plots of large-scale MHD and turbulence quantities for b = 1/2 in the same format as Figure 8 in W15. As the additional term is proportional to the Alfvén-velocity, it drops out for the model validation performed in Section 3 of W15, which neglected such terms for the case kUk kVA k valid beyond some 0.3 AU of heliocentric distance. The term would have appeared in the following Section 4 of W15 though, where the model was extended. However, as the extended model was first applied to regimes of highly 58 4.3. FURTHER DEVELOPMENT super-Alfvénic winds too, where these terms have little effect, the additional term hardly changed the results at all so that they are not shown here. It becomes more interesting for the simulations of the very inner heliosphere including a CME (Section 5 in W15), as here the Alfvén speed is not negligible compared to the solar wind speed. Figures 4.2 and 4.3 show the newly obtained results, which are in the same format as Figures 10 and 8 in W15, respectively. The line plots compare the large-scale MHD quantities only, and there is obviously almost no change at all. This is in agreement with the findings in W15 that the turbulence quantities have little to no effect on the large-scale flow, so that small changes in the turbulence quantities are even less recognizable in the large-scale quantities. The contour plots also remain remarkably similar, but on close inspection the absolute values of the turbulence quantities have changed noticeably, so that Figure 4.4 shows ratios for the turbulence quantities calculated via (X1/2 /X0 − 1) · 100, where X stands for the respective quantity with the subscripts denoting the value of b, and the resulting value is in per cent. It can be seen that the total turbulent energy density is about 35% lower away from the CME, while close-by changes are much smaller (about 1%). This transfers directly to δB/B as well, where the difference is only about 20% though. Also, the (absolute) cross helicity is about 10% lower, but can spike within the CME, and the correlation length is about 20% larger. Figure 4.4: Meridional contour plots of the ratios of turbulence quantities between b = 1/2 and b = 0 (see text). It can be concluded that the additional term has quite an effect in the very inner heliosphere, but not beyond some 0.3 AU. However, especially since the results directly at the CME are only slightly affected, the general findings presented in W15 remain valid. 59 CHAPTER 4. IMPLEMENTING TURBULENCE TRANSPORT 4.4 Addendum - Conservative treatment of the co-rotating frame of reference It is often desirable to perform the simulations in a co-rotating frame of reference to simplify boundary conditions. Implementing the resulting fictitious (Coriolis-, centrifugal-, and Euler) forces as source terms can lead to instabilities and nonconservation of angular momentum (Kley, 1998; Mignone et al., 2012), in particular for large heliocentric distances, because the azimuthal velocity component increases linearly with distance. The force terms can, however, be partially included in the time derivative and flux function of the hyperbolic momentum equation, which amounts to a quasi-conservative treatment and results in the following set of equations. These are equivalent to the original ideal MHD equations with fictitious forces as source terms, such that the new set ∂t ρ + ∇ · (ρv) = 0 ∂t (ρu) + ∇ · ρvu + (p + kBk2 /2)I − BB = −ρΩ × u ∂t e + ∇ · ev + (p + kBk2 /2)u − (u · B)B = 0 (4.2) (4.3) (4.4) replaces Eqs. (1.1–1.3), and Ohm’s law (1.6) is now E+v×B=0 . (4.5) Here, v denotes the velocity in the co-rotating frame, which is related to the inertial frame velocity u via v =u−Ω×r . (4.6) This can be derived as follows: using the latter relation and the vector identities ∇ · (FG) = G(∇ · F) + (F · ∇)G, (F · ∇)(Ω × r) = Ω × F, where F and G are arbitrary vectors, and ∇ · (Ω × r) = 0, it can be shown that ∂t (ρv) + ∇ · [ρvv] + 2ρΩ × v + ρΩ × (Ω × r) = ∂t (ρu) + ∇ · [ρvu] + ρΩ × u , (4.7) so that the non-conservative momentum equation is reproduced. On using Equation (4.3) one gets for the time derivative of the inertial frame total energy e = p/(γ − 1) + ρkuk2 /2 + kBk2 /2 the above energy equation ∂t e + ∇ · ev + (p + kBk2 /2)u − (u · B)B = 0 (4.8) instead of the time derivative of the co-rotating frame total energy e = p/(γ − 1) + ρkvk2 /2 + kBk2 /2, which yields ∂t e + ∇ · (e + p + kBk2 /2)v − (v · B)B = − v · (2ρΩ × v + ρΩ × (Ω × r)) . 60 (4.9) Chapter 5 Extended turbulence transport model The next step towards a more sophisticated turbulence transport model is taken by implementing the full set of equations as given by Zank et al. (2012), which now also involves a variable energy difference σD , as well as a corresponding correlation length λD and separate correlation lengths for the forward and backward propagating modes λ± . Section 5.1 describes the necessary alterations to the original equations in order to be implemented in the Cronos framework for a more adequate formulation exploiting the conservative features of the code. The implementation is validated in Section 5.2 by comparing with the work of Adhikari et al. (2015), who applied a number of simplifications, however. After successful validation, it is then possible to extend the model further by considering additional terms that have been left out thus far (Section 5.3). 5.1 Model formulation In analogy with the previous chapter, these equations are reformulated in the corotating frame of reference. The assumed symmetry of the turbulence is arbitrary for now by maintaining the general form of the structural similarity parameters a, b. Furthermore, instead of the evolution equations for the correlation lengths λD and λ± , it was found that the corresponding equations for the integral scales defined as L± ≡ λ± Z 2 (1 ± σC ) , LD ≡ λD Z 2 σD (5.1) are numerically better suited and are used for the computations. This was found to be particularly the case for LD , because in some cases λD → ∞ while σD → 0 (see Figure 5.3). The turbulence transport equations solved in Cronos can then 61 CHAPTER 5. EXTENDED TURBULENCE TRANSPORT MODEL be written as ∂t Z 2 + ∇ · (VZ 2 + VA Z 2 σC ) = (Z 2 /2)∇ · U + 2VA · ∇(Z 2 σC ) − (2a − 1/2)Z 2 σD ∇ · U + 2aZ 2 σD n̂ · (n̂ · ∇)U √ √ q 1 + σC 1 − σC 3 2 − Z 1 − σC + λ+ λ− + hz+ · S+ i + hz− · S− i ∂t (Z 2 σC ) + ∇ · (VZ 2 σC + VA Z 2 ) = (Z 2 σC /2)∇ · U + 2VA · ∇Z 2 + Z 2 σ D ∇ · VA √ + (2bZ 2 σD / ρ)n̂ · (n̂ · ∇)B √ √ q 1 + σC 1 − σC 3 2 − Z 1 − σC − λ+ λ− + hz+ · S+ i − hz− · S− i ∂t (Z 2 σD ) + ∇ · (VZ 2 σD ) = (Z 2 /2)(σD − 4a + 1)∇ · U − Z 2 σC ∇ · VA + 2aZ 2 n̂ · (n̂ · ∇)U (σC VA · ∇Z 2 − VA · ∇(Z 2 σC )) p − 1 − σC2 √ − (2bZ 2 σC / ρ)n̂ · (n̂ · ∇)B √ √ 1 + σC 1 − σC 3 + − Z σD λ− λ+ + hz− · S+ i + hz+ · S− i ∂t (ρL± ) + ∇ · (VρL± ) = ρ[ ± VA · ∇L± − ∇ · (U/2 ± VA )L± − ((a − 1/4)∇ · U ∓ ∇ · VA /2)LD + (an̂ · (n̂ · ∇)U √ ± (b/ ρ)n̂ · (n̂ · ∇)B)LD ] ∂t (ρLD ) + ∇ · (VρLD ) = −ρ[ ∇ · U(LD /2) + ∇ · VA (L+ − L− ) r r L+ L− − + V · ∇L − VA · ∇L+ A L− L+ + [(2a − 1/2)∇ · U − 2an̂ · (n̂ · ∇)U](L+ + L− ) √ + (2b/ ρ)[n̂ · (n̂ · ∇)B](L+ − L− )] . (5.2) (5.3) (5.4) (5.5) (5.6) As some of the involved quantities have been introduced in the previous chapter already, there are a few repetitions in what follows in order to have a complete description of the model here. √ The Elsässer variables are z± := u ± b/ ρ with u and b denoting the fluctuations about the mean inertial velocity (U) and magnetic (B) fields, while ρ is the mass density and all quantities and equations are given in their normalized form. The 62 5.1. MODEL FORMULATION following moments of the Elsässer variables are used above: hz+ · z+ i + hz− · z− i = hu2 i + hb2 /ρi 2 hz+ · z+ i − hz− · z− i √ Z 2 σC := = 2hu · b/ ρi 2 2 + − Z σD := hz · z i = hu2 i − hb2 /ρi . Z 2 := (5.7) (5.8) (5.9) Thus, Z 2 is twice the total (kinetic plus magnetic) energy per unit mass of the fluctuations, σC is the normalized cross-helicity, and σD is twice the normalized energy difference in the fluctuations per unit mass, also known as residual energy. √ Furthermore, VA = B/ ρ is the Alfvén velocity, and terms involving sources of turbulence S± are discussed in subsequent sections. Note that the dissipation term (second last line) in the equation for Z 2 σD is not the one used by Zank et al. (2012) (based on a model of Müller and Grappin (2005)), but the ones introduced by Dosch et al. (2013), which has also been used in Adhikari et al. (2015), because the original dissipation term was found to lead to unphysical results both in the present setup as well as in others (N. Pogorelov, priv. comm). The terms involving the structural similarity parameter b in the original equations of Zank et al. (2012) have been rewritten above (S. Oughton, priv. comm.) on using the identity 1 2∇ · VA + [VA · ∇ρ − (n̂ · VA )(n̂ · ∇ρ)] ρ 2 =2n̂ · (n̂ · ∇)VA − √ n̂ · (n̂ · ∇)B , ρ √ which follows from using the product rule on VA = B/ ρ via combined with (5.10) 1 1 n̂ · (n̂ · ∇)VA = n̂ · (n̂ · ∇)B − (n̂ · VA )(n̂ · ∇)ρ ρ 2ρ (5.11) 1 2∇ · VA = − VA · ∇ρ . ρ (5.12) Equation (5.10) can be applied to Equations (42) and (43) in Zank et al. (2012) for the energy density of the forward and backward propagating modes z ± , which allows for a direct comparison with the corresponding Equations (21) and (22) of Matthaeus et al. (1994). On performing this comparison, it becomes obvious that the case a = b = 1/2 corresponds to 2D axisymmetric turbulence, while Zank et al. (2012) states that a = 1/2, b = 0 in that case. Regardless of this, the axisymmetry is usually imposed along the direction of the mean magnetic field so that n̂ = B̂. For isotropic turbulence, a = 1/3, b = 0 holds, while n̂ = 0 in this case, which, however, is not considered in the following. 63 CHAPTER 5. EXTENDED TURBULENCE TRANSPORT MODEL The turbulence transport Equations (5.2) – (5.6) are implemented in the framework of the Cronos code as additional equations to be solved alongside the Reynoldsaveraged MHD equations in the co-rotating frame of reference with respective coupling terms to account for the effects of turbulence on the large-scale MHD quantities (see Wiengarten et al., 2015, (W15)): ∂t ρ + ∇ · (ρV) = 0 (5.13) ∂t (ρU) + ∇ · [ρVU + p̄ 1 − ηBB] = −ρ (g + Ω × r) (5.14) ∂t B + ∇ · [VB − BV] = 0 (5.15) 2 2 ∂t e + ∇ · eV + (p + |B| /2) U − (U · B)B − VA ρZ σC /2 + qH = Z 2 σC − ρV · g − U · ∇pw − VA · ∇ρ − ρVA · ∇(Z 2 σC ) 2 √ √ q ρ 3 1 + σ 1 − σ C C + Z 1 − σC2 + 2 λ+ λ− (5.16) + U · (B · ∇)[(η − 1)B] + ρZ 2 σD (a − 1/2)∇ · U with p̄ = (p + |B|2 /2 + pw ), pw = (σD + 1)ρZ 2 /4 and η = 1 + σD ρZ 2 /(2B 2 ) for transverse and axisymmetric turbulence, while pw = (σD + 3)ρZ 2 /12 and η = 1 for isotropic turbulence. Equations (5.13) – (5.16) are equivalent to those of W15, while the energy equation now also retains the structural similarity parameter a. Furthermore, V = U − Ω × r is the fluid velocity in the co-rotating frame, p is the scalar thermal pressure, g = (GM /r2 )r̂ is the solar gravitational acceleration, and Ω = Ωez is the Sun’s angular rotation velocity with Ω = 14.71◦ /d (Snodgrass and Ulrich, 1990). e = ρU 2 /2 + B 2 /2 + p/(γ − 1) is the energy density without the turbulent energy component ρZ 2 /2, which is instead included by means of the source terms in Equation (5.16). An adiabatic equation of state is used with γ = 5/3 while, due to the inclusion of Hollweg’s heat flux (Hollweg, 1974, 1976) qH = (3/4)pV, the effective value of the adiabatic index (see also Usmanov et al., 2011) is γef f = 13/9, which is close to observationally inferred values. 5.2 Model validation In W15 the code was validated by comparing with the results of Usmanov et al. (2011), who solved a reduced set of turbulence transport equations alongside the Reynolds-averaged MHD quantities. As the coupling of the turbulence to the MHD equations from my previous work is retained, only the implementation of the new extended set of turbulence transport equations (5.2) – (5.6) had to be validated. These equations have been solved numerically by Adhikari et al. (2015), where the calculations were restricted to the ecliptic plane ϑ = π/2, neglected Alfvén velocity for the outer heliosphere calculations, constant solar wind velocity U = 400km/s r̂, and the shear-mixing term B̂ · (B̂ · ∇)U was left out. Furthermore, the structural similarity parameters were adopted as a = 1/2, b = 0 in that work. Therefore, to 64 5.2. MODEL VALIDATION validate the implementation these simplifying assumptions are also taken for now. In accordance with Adhikari et al. (2015) the following expressions for the turbulence source terms due to stream shear and isotropization of pickup ions are used: ∆Ushear 2 Z (1 + σC ) + r ∆Ushear 2 hz− · S− i = Cshear Z (1 − σC ) + r ∆Ushear 2 Z σD hz− · S+ i + hz− · S+ i = 2Cshear r hz+ · S+ i = Cshear Ėpui 2 Ėpui 2 (5.17) (5.18) (5.19) with a tuning factor Cshear = 7.25, the difference between fast and slow speed streams at 0.29 AU is ∆Ushear = 350 km/s and −Lcav f D nH U V A exp (5.20) Ėpui = nsw τion r where fD = 0.25 is the fraction of pickup ion energy transferred into excited waves, nH = 0.1 cm−3 is the interstellar neutral hydrogen density, τion = 106 s is the neutral ionization time at 1 AU, Lcav = 8 AU is the characteristic scale of the ionization cavity of the Sun, and nsw = 5 cm−3 is the solar wind density at 1 AU adopted in this model. Although neglected in the governing equations, VA = 50 km/s here. The fixed inner boundary values are Z02 = 7134 (km/s)2 , σC,0 = 0.9, σD,0 = −0.004, − λD,0 = 0.0204 AU, λ+ 0 = 0.000779 AU and λ0 = 0.00143 AU. Note: Adhikari et al. (2015) use a different definition for the cross-helicity (see Equation (10) there) resulting in a positive cross-helicity for a positive (outwardly directed) radial magnetic field. Here, an inwardly directed field with the usual definition for cross-helicity as given above is adopted, thus also getting positive cross-helicity values. However, the associated correlation lengths have to be reversed, i.e. λ± = λ∓ Adhikari . The computational domain covers the radial range from 0.29 AU to 100 AU with 2400 cells of increasing radial cell size ∆r ∈ [0.4, 43]R . This extremely high resolution in particular at the inner boundary was necessary to retain a stable solution due to the imposed strong shear terms. Rotational symmetry allows to use just one cell covering the azimuth ϕ, and since for this test case the Alfvén velocity (and, thus, the magnetic field) is neglected, there is also no preferred direction along the polar angle, so that a fairly low resolution can be applied in ϑ. Steep gradients near the inner boundary due to large dissipation terms require (global) time-steps that have to be considerably smaller than the usual constraint of the CFL condition calculated from the smallest cell-size and the characteristic velocities vc there. While usually a value of cf lthresh = 0.4 can be chosen in dt = cflthresh · min(∆xi /vc ) it was necessary to reduce it to cflthresh = 0.03. Such a time-step constraint can also be computed self-consistently with methods such as sub-cycling, where several sub-steps for the source terms are taken in-between a larger MHD time-step. However, since a timestep reduction was found only to be necessary in the first few inner radial cells, it 65 CHAPTER 5. EXTENDED TURBULENCE TRANSPORT MODEL was chosen to increase cflthresh after a steady state was reached in these cells, which takes a relatively short time only. Figure 5.1: Comparison of my results (black curves) with Adhikari et al. (2015) (red) for the turbulent energy Z 2 (top left), the normalized cross-helicity σC (top right), turbulent energy density for forward and backward propagating modes z ± (left center) and respective correlation lengths λ± (right center), the residual energy σD (bottom left) and its correlation lengths λD (bottom right). The results (black lines) and the reference solution of Adhikari et al. (2015) (red lines) are shown in Figure 5.1 for the turbulent energy density Z 2 , the normalized cross-helicity σC , turbulent energy density for forward and backward propagating modes z ± and respective correlation lengths λ± , residual energy σD and its correlation length λD . The turbulent energy dissipates with increasing radial distance, and is transferred into heat (temperature not shown here). The rate of dissipation is controlled by the correlation lengths λ± . Because of the low boundary values of the latter, the initial dissipation is very strong. In the outer heliosphere beyond the ionization cavity turbulence is generated by isotropization of pickup ions and the turbulent energy density becomes almost constant. The cross-helicity gives the ratio of forward (z + ) to backward (z − ) propagating modes, where the orientation is anti-parallel to the mean magnetic field. Thus, for an inwardly directed mean magnetic field, only backward propagating modes can escape the sub-Alfvénic solar wind regime below the Alfvén critical point, and the cross-helicity σC . 1 towards 66 5.3. MODEL EXTENSIONS the inner boundary. With the generation of turbulence in the supersonic solar wind, also forward propagating modes are generated and equipartition is reached with increasing radial distance, i.e. σC → 0 and z + = z − . The residual energy gives the difference between the kinetic and magnetic fluctuations’ strength. In case of equipartition the turbulence is called Alfvénic and σD ≈ 0. An interesting result of Adhikari et al. (2015) is that the magnetic fluctuations become increasingly dominant out to about 10 AU, whereas beyond this point the fluctuations become more Alfvénic again. The correlation lengths control the rate of dissipation of turbulent energy. The correlation length of the forward propagating modes λ− is initially longer (i.e. the dissipation is weaker) than for the backward propagating modes, because the backward modes contain the higher energy. With increasing radial distance and turbulence generation the correlation lengths both rise gradually and become equal beyond about 10 AU due to the equipartition of the respective modes. The correlation length of the residual energy is steadily increasing, gradually out to 10 AU and more rapidly thereafter, as the residual energy tends toward zero, and increasingly less energy is ’dissipated’ into the energy difference. The excellent agreement between my results and the reference values demonstrates the correct implementation, so that it can now be proceeded to include additional terms left out thus far. 5.3 Model extensions The calculations can easily be extended on the basis of the general equations of the previous section that have been fully implemented. In the following, the general effects of including these terms are demonstrated. For a sophisticated tuning of the parameters (boundary conditions, strength of the shear and pickup ion driving) a thorough comparison with spacecraft data would be required. This could not be achieved during the work for this thesis, but it should be a mandatory next step. An arguable approach was taken by Adhikari et al. (2015) concerning the shear driving (Equations (5.17)–(5.18)), as the chosen form for the shear- generated turbulence is respectively proportional to the forward and backward propagating modes, so that in this case the initially stronger backward propagating modes are further enhanced. This leads to the initial increase in the cross-helicity (top right panel in Figure 5.1). If instead a shear driving term that is proportional to the total turbulent energy density is chosen, the respective source terms become ∆Ushear 2 Ėpui Z + (5.21) r 2 ∆Ushear 2 Ėpui hz− · S− i = Cshear Z + . (5.22) r 2 Otherwise maintaining the previous setup, the results with the new shear terms are shown in Figure 5.2, which are quite similar to the reference case. However, the hz+ · S+ i = Cshear 67 CHAPTER 5. EXTENDED TURBULENCE TRANSPORT MODEL Figure 5.2: Same as Figure 5.1, but with shear driving proportional to the total turbulent energy density. cross-helicity is now monotonically decreasing and the steep gradients at the inner boundary are more gradual now, which seems to be more realistic. Also, the value for the energy difference is now closer to −1/3 at 1 AU as seen in observations (see W15 and references therein). Therefore, the new source terms are probably better suited, however, a thorough tuning of Cshear is still necessary as well as a latitude-dependent value for the difference between fast and slow streams ∆Ushear , because this term mimics CIRs that are mainly present at equatorial latitudes only. Furthermore, the pickup ion source term can also be augmented (Isenberg, 2005). Such advances should be incorporated in future modeling. Another simulation has been performed, in which as a first step the magnetic field, and thus terms involving VA , are taken into account. The magnetic field is prescribed as in W15, i.e. in the form of a Parker spiral, but in a mono-polar fashion (inwardly directed at all latitudes) to avoid the influence of a current sheet. It comes as no surprise that the results are barely changing (and are, thus, not shown) in this outer-heliospheric scenario, where the approximation U VA is valid. However, this result is a respective validation of this frequently made assumption. In the following, terms involving VA are always taken into account. A more interesting extension (again otherwise maintaining the setup from the previous section, i.e. old shear terms) involves the shear-mixing term B̂ · (B̂ · ∇)U, for 68 5.3. MODEL EXTENSIONS Figure 5.3: Same as Figure 5.1 with the shear-mixing term B̂ · (B̂ · ∇)U. which the magnetic field is also not neglected. (The magnitude of B is irrelevant for this term, however). It was already shown in W15 that this term is latitudedependent, so that in Figure 5.3 the radial evolution at selected co-latitudes differs now. This becomes obvious from simplifying this term for the present underlying assumption of a constant solar wind speed and Parker spiral magnetic field, yielding B̂ · (B̂ · ∇)U ≈ (Bϕ /B)2 (U/r) so that via Bϕ a dependence on latitude is introduced. Accordingly, the effects of this term are most clearly seen at equatorial latitudes (dashed-dotted). Furthermore, while within about 2 AU the azimuthal magnetic field component is still small, it becomes dominant beyond – due to its 1/r dependence as opposed to the 1/r2 one of the radial component – which explains that differences occur mainly beyond 2 AU. As compared with the reference case, major discrepancies are now found for the cross-helicity, which does not approach zero in the outer heliosphere, the stronger so towards the ecliptic. Thus, an equipartition between forward and backward propagating modes is not reached, as can also be seen from the center panels for the respective turbulent energy densities and correlation lengths. Furthermore, the energy difference quickly goes to zero as soon as pickup ion driving sets in, and the respective correlation length diverges (which made the computation with LD instead of λD necessary, as mentioned above). A more elaborate analysis of the effects of this additional term is again beyond the scope of this work, however, its influence is evident and not negligible. 69 CHAPTER 5. EXTENDED TURBULENCE TRANSPORT MODEL Figure 5.4: Same as Figure 5.1 with b = 1/2. As outlined in Section 5.1, there is doubt about the chosen value b = 0 for 2D axisymmetric turbulence, and instead a value of b = 1/2 is deemed appropriate. This introduces a similar term as the previously described one into play, namely ρ−1/2 B̂ · (B̂ · ∇)B, and yet another simulation was performed incorporating only this change as compared to the reference case. As with the previous term, its effects become evident only beyond about 2 AU, where the magnetic field becomes increasingly perpendicular to the flow direction. The results as shown in Figure 5.4 are quite intriguing as they show increases in the turbulent energy density, but no equipartition between forward and backward propagating modes (indeed, σC → 1) and also no equipartition between kinetic and magnetic fluctuations as σD → const. 6= 0. Finally, a setup was considered in which all the additional terms presented in this section were included, and the results are shown in Figure 5.5. As one was tempted to believe that such a seemingly complete model would result in reasonable results, the found behavior is rather disappointing. Particularly a strong increase in the turbulent energy density and no equipartition between forward and backward propagating modes is in conflict with the common conception based on observational data. Furthermore, this simulation became unstable beyond some 70 AU (see center panels). Besides the aforementioned next steps that should be taken, there is also criticism to 70 5.3. MODEL EXTENSIONS Figure 5.5: Same as Figure 5.1 with b = 1/2, shear-mixing term and new source terms. these kind of models taking into account 2D turbulence only, because especially the pickup ion isotropization in the outer heliosphere rather drives Alfvén waves (slab turbulence). It is often assumed for modeling of diffusion coefficients based on the magnetic fluctuations that the ratio between the two is δB2D /δBslab ≈ 4 (i.e. 80% 2D, 20% slab; Bieber et al., 1996), so that the slab component is smaller but not negligible. Models of turbulent transport usually consider only one or the other, but there are also approaches for combined theories, such as the two component model byOughton et al. (2011). In the framework of a planned upcoming project with S. Oughton, this model will also be implemented in Cronos and investigated in the near future. Nevertheless, considerable progress is expected to be made from this general implementation of the model of Zank et al. (2012), as in the present setup it can be solved for the first time in 3D and also in a time-dependent manner, and it can also be coupled self-consistently to more realistic background solar wind conditions, particularly to the observations based simulations driven by the WSA model (Wiengarten et al., 2014). For the latter, it will be necessary to devise empirical formulas for the turbulence values at the inner boundary in a similar fashion as in the current WSA setup, i.e. depending on coronal conditions quantified by, e.g., the flux tube expansion factor. 71 CHAPTER 5. EXTENDED TURBULENCE TRANSPORT MODEL 72 Chapter 6 Modeling of astrospheres 6.1 Overview Global heliospheric models – focusing on the interaction between the solar wind and the ISM – and respective MHD simulations with the Cronos code are another major topic in our working group. For the heliosphere such a numerical setup was developed by J. Kleimann, and in a recent paper (Röken et al., 2015) these simulations were compared to a new analytic description for the magnetic field in the vicinity of the heliopause. In recent years, observational advancements (e.g. Hubble Space Telescope, Widefield Infrared Survey Explorer (WISE) mission) allowed for all-sky surveys measuring dust emission (mostly infrared; e.g., Groenewegen et al., 2011) or the Lyman-α line (Linsky and Wood, 2014) associated with bow-shocks around other stars. It can be concluded that these stars are also emitting a stellar wind interacting with the ISM, consequently forming a respective astrosphere depicting the same bullet-like shape as the heliosphere, and also containing structures such as a termination shock, an astropause, astrosheaths and the already observed bow-shock (an overview is given in Section 2 of Scherer et al., 2015, S15). Therefore, many findings about the solar wind and our heliosphere can be transferred to other stars and their astrospheres, and the numerical setup devised for the heliosphere is already principally suited for respective investigations. Care must be taken though, as the interstellar environment around those stars is most likely different from our local environment, and the processes driving a stellar wind may be very different from the thermally/wave-driven solar wind (see Lamers and Cassinelli (1999) for an overview). For example, the wind emitted by hot (O/B-)stars is line driven, and the energy transfer within these stars is completely radiative, so that it is still an open question whether a respective stellar magnetic field is present (although some examples have recently been found (Hubrig et al., 2015)). The resulting stellar wind velocities can be as high as several thousand km/s, and these stars’ astrospheres can be several parsecs wide. Interestingly, the 73 CHAPTER 6. MODELING OF ASTROSPHERES much cooler M-dwarf stars have also been shown to emit such fast stellar winds, and as these stars have an almost completely convective interior, they are magnetically very active and a considerable amount of the mass loss is probably attributable to flares and CMEs. For nearby stars, it is even possible to infer their photospheric magnetic field distribution (Donati et al., 2008). Respective data has been used in, e.g., Vidotto et al. (2011, 2015); Vidotto (2014) on the basis of solar and heliospheric numerical models with a focus on M-dwarfs with orbiting exoplanets to investigate possible interactions between these stellar winds and planetary magnetospheres. An effective shielding against the stellar wind’s influence should be regarded as an important factor for habitability. Similarly, the size of a respective astrosphere and its capability to shield planets from cosmic rays may also be important, so that global astrospheric models for these cases are a possible interesting future direction in our group. A different goal was pursued in the publication Cosmic rays in astrospheres (Scherer et al., 2015) presented in this chapter, where the focus is on hot stars, because the associated astrospheres are huge cavities that modulate the cosmic ray flux through them and might therefore be responsible for tiny-scale anisotropies in all-sky maps of the GCR flux. This work was done in co-operation with colleagues of the Astronomical Institute at the Ruhr University Bochum, who provided observational data for the investigated O-star λ Cephei, from which some of its stellar wind parameters could be derived, while some educated guesses have to be made for observationally inaccessible quantities like the wind’s temperature (Section 4 of S15). Together with these inner boundary conditions, outer boundary conditions in the chosen spherical setup need to be chosen according to the ISM values (also educated guesses), and special attention has to be paid to make sure that outflow in the tail direction is ensured. The chosen set of equations to describe this scenario in this case is hydrodynamic and assumes a single (proton) fluid. More sophisticated models including magnetic fields and a multi-fluid approach – particularly accounting for interstellar neutrals – are currently being implemented. One important effect that neutrals or partially ionized heavy elements have on the plasma flow in the outer astrosheath is that they can effectively cool the plasma via Coulomb collisions, a process for which the length scale is smaller than the extent of the astrosheath and, thus, has to be taken into account (as opposed to the heliosphere, where this effect is negligible). Respective cooling terms can be introduced in the hydrodynamic equations without the necessity for a explicit multi-fluid approach. Meanwhile, heating terms were also considered, such as photo-ionization due to the large photon flux emitted by the star. Both heating and cooling are discussed in Section 3 of S15. The resulting MHD configuration was then used as input to an SDE model to analyze cosmic ray fluxes modulated by this astrosphere (Section 5 of S15). This was done in co-operation with our colleagues from South Africa. The results demon74 6.1. OVERVIEW strated that cosmic rays are strongly influenced at lower energies, but also at very high energies in the TeV range a small modulation of a few per mille was found as compared to cosmic rays not encountering this astrospherical obstacle. Such tiny-scale anisotropies in all-sky maps of measured high-energy cosmic ray fluxes at Earth are observed, e.g. with the IceCube experiment (see Figure 6.1), so that one possible explanation might be the modulation via nearby large astrospheres as considered here. Figure 6.1: IceCube cosmic ray map for the 400 TeV energy band showing tiny-scale anisotropies (taken from Abbasi et al., 2012). 75 6.2 Scherer et al. (2015) Cosmic rays in astrospheres Astronomy & Astrophysics 76 Astronomy & Astrophysics A&A 576, A97 (2015) DOI: 10.1051/0004-6361/201425091 c ESO 2015 Cosmic rays in astrospheres K. Scherer1,2 , A. van der Schyff3 , D. J. Bomans4,2 , S. E. S. Ferreira3 , H. Fichtner1,2 , J. Kleimann1 , R. D. Strauss3 , K. Weis4 , T. Wiengarten1 , and T. Wodzinski4 1 2 3 4 Institut für Theoretische Physik IV: Weltraum- und Astrophysik, Ruhr-Universität Bochum, 44780 Bochum, Germany e-mail: [kls;hf;jk;tow]@tp4.rub.de Research Department, Plasmas with Complex Interactions, Ruhr-Universität Bochum, 44780 Bochum, Germany Center for Space Research, North-West University, 2520 Potchefstroom, South Africa e-mail: [12834858;Stefan.Ferreira;DuToit.Strauss]@nwu.ac.za Astronomisches Institut, Ruhr-Universität Bochum, 44780 Bochum, Germany e-mail: [bomans;kweis;thomas.wodzinski]@astro.rub.de Received 1 October 2014 / Accepted 10 February 2015 ABSTRACT Context. Cosmic rays passing through large astrospheres can be efficiently cooled inside these “cavities” in the interstellar medium. Moreover, the energy spectra of these energetic particles are already modulated in front of the astrospherical bow shocks. Aims. We study the cosmic ray flux in and around λ Cephei as an example for an astrosphere. The large-scale plasma flow is modeled hydrodynamically with radiative cooling. Methods. We study the cosmic ray flux in a stellar wind cavity using a transport model based on stochastic differential equations. The required parameters, most importantly, the elements of the diffusion tensor, are based on the heliospheric parameters. The magnetic field required for the diffusion coefficients is calculated kinematically. We discuss the transport in an astrospheric scenario with varying parameters for the transport coefficients. Results. We show that large stellar wind cavities can act as sinks for the Galactic cosmic ray flux and thus can give rise to small-scale anisotropies in the direction to the observer. Conclusions. Small-scale cosmic ray anisotropies can naturally be explained by the modulation of cosmic ray spectra in huge stellar wind cavities. Key words. stars: winds, outflows – cosmic rays – hydrodynamics 1. Introduction Recently, simulations of astrospheres around hot stars have gained new interest, see for example Decin et al. (2012), Cox et al. (2012), Arthur (2012), van Marle et al. (2014). These authors modeled astrospheres using a (magneto-)hydrodynamic approach, either in 1D or 2D. In this work, such astrospheric models are used for the first time to estimate the cosmic ray flux (CRF) through it. Because of the large spatial extent of O star astrospheres (wind bubbles), these objects can efficiently cool the spectrum of Galactic cosmic rays (GCR). Especially λ Cephei is an interesting example, being the brightest runaway Of star in the sky (type O6If(n)p). We estimate the CRF at different energies for λ Cephei as an example of an O star astrosphere using stochastic differential equations (SDE; Strauss et al. 2013) to solve the GCR transport equation. Runaway O and B stars are common and part of a sizable population in the Galaxy, a significant number of which show bow-shock nebulae (e.g., Huthoff & Kaper 2002). In Sect. 3 we discuss the radiative cooling functions, while in Sect. 4 we show the astrosphere model results. In Sect. 5 we estimate the CRF. 2. Large-scale structure of astrospheres Winds around runaway stars, or in general, stars with a nonzero relative speed with respect to the surrounding interstellar medium (ISM), develop bullet-shaped astrospheres. The hydrodynamic large-scale structure is sketched in Fig. 1, the notation of which is described below. The hypersonic stellar wind (Mach numbers Ma 1) undergoes a shock transition to subsonic velocities at the termination shock (TS) in the inflow direction. Then a tangential discontinuity, the astropause (AP), is formed between the ISM and the stellar wind, where the velocity normal to it vanishes: there is no mass transport through the AP. Other quantities such as the tangential velocity, temperature, and density are discontinuous, while the thermal pressure is the same on both sides. If the relative speed, or the interstellar wind speed as seen upwind in the rest frame of the star, is supersonic in the ISM, a bow shock (BS) exists. If the relative speed is subsonic, there will be no BS, see Table 1 for the stellar parameters and for the stellar-centric model distances of the TS, AP, and BS. The region between the BS and AP is called outer astrosheath, the region between the AP and TS the inner astrosheath. The AP around the inflow direction at the stagnation line is sometimes called the nose, while the region beyond the downwind TS is called the astrotail. The latter can extend deep into the ISM. The region inside the TS is called the inner astrosphere. In the downwind direction, the termination shock forms a triple point (T), from which the Mach disk (MD) extends down to the stagnation line; this is the line through the stagnation point and the star. A tangential discontinuity (TD) emerges down into the tail. A reflected shock (not shown here) also extends from the triple point toward the TD. Article published by EDP Sciences A97, page 1 of 7 CHAPTER 6. MODELING OF ASTROSPHERES 78 Chapter 7 Summary and outlook This thesis focuses on the modeling of the inner-heliospheric solar wind by numerically solving the equations of magnetohydrodynamics (MHD). Consequently, the first chapter starts with a descriptive overview of the respective protagonists, namely the Sun with its magnetic field, emanating and shaping the turbulent solar wind that interacts with the interstellar medium to form the heliosphere, and the application of the solar findings to other stars and their astrospheres (Section 1.1). Furthermore, aspects of sub-dividing the heliospheric environment for modeling purposes are given, before the main numerical tool used in this work – the state-of-the-art MHD code Cronos – is described in some detail in Section 1.2. The remaining chapters are respectively based on the publications I authored during my PhD, with respective introductions, additional details and afterthoughts leading to subsequent working areas. Chapter 2 – including my initial first-author publication MHD simulation of the inner-heliospheric magnetic field (Wiengarten et al., 2013) – presents the application of the Cronos code to the modeling of the solar wind based on observational input data. The required extensions to the code for adequately describing the solar wind include the introduction of solar rotation and solving the full energy equation, which was validated by reproducing the analytic models of Parker (1958) and Weber and Davis (1967). Although these latter models have been very successful in describing the large-scale solar wind and the heliospheric magnetic field, they are too simple for reproducing more complex conditions, especially those involving transient events such as corotating interaction regions (CIRs) or coronal mass ejections (CMEs). Therefore, in a collaboration with colleagues from the Max Planck Institute for Solar System Research (MPS) in Katlenburg-Lindau (now in Göttingen), maps of the radial magnetic field at the heliobase for solar minimum and maximum conditions served as inner boundary conditions for my subsequent MHD simulations. These maps are the result of a coupled solar surface flux transport (SFT) and a coronal potential field model as described in Jiang et al. (2010), which allow for an adequate description 79 CHAPTER 7. SUMMARY AND OUTLOOK of the usually unobservable far side of the Sun, as well as a reconstruction of the solar magnetic field for times before the space age by using sunspot group data from earlier times. Previously, the Cronos code was operated with analytic prescriptions for initial and inner boundary conditions only, so that new routines to read-in and interpolate the input data on arbitrary grids had to be devised. Moreover, a method for a divergence-free initialization of the magnetic field was developed. Besides the directly provided magnetic field data, the remaining MHD quantities had to be derived as well, which follow from the coronal magnetic field topology – fast wind from coronal holes and slow wind from their boundaries – as described by respective models relying on the inverse relation between flux-tube expansion factor and resulting solar wind speed (Wang and Sheeley, 1990). Here, the model by Detman et al. (2011) was adopted. Simulations were performed on clusters at our institute for typical solar minimum and maximum conditions. Initially, the obtained results were deemed as realistic as they were roughly compared to spacecraft data and agreed with the general topologies expected for the respective phases in the solar cycle. Later on, a more detailed comparison with spacecraft data showed that the magnetic field data used in this work were unsuitable, and the coronal magnetic field modeling was subsequently performed in a different manner by myself, as described next. The work presented in Chapter 3 was done in the framework of a joint project with the Universities of Kiel and Potchefstroom, South Africa, on the influence of CIRs on the propagation of energetic particles, with the goal to use a CIR disturbed background solar wind configuration resulting from my MHD simulations as input to a stochastic differential equation (SDE) model for the transport of energetic particles in the heliosphere. Whereas the publication Cosmic Ray transport in heliospheric magnetic structures. I. Modeling background solar wind using the Cronos magnetohydrodynamic code (Wiengarten et al., 2014) describes the respective MHD modeling, the subsequent SDE modeling paper on energetic particle propagation is currently in preparation. A first mandatory task was to show that the Cronos code is able to correctly compute CIR structures, which was achieved by reproducing previous numerical results by Pizzo (1982). Our colleagues from Kiel proposed the time period around August 2007 during the ascending phase of solar cycle 23 for a detailed comparison with spacecraft measurements. Besides favorable solar wind conditions with a stable CIR configuration, this time is also preferential due to a unique spacecraft constellation with the just previously launched Stereo-A/B spacecraft trailing and leading Earth on its orbit at 1 AU, and Ulysses crossing the equatorial plane during its third fastlatitude scan. The required potential field modeling of the solar corona providing the boundary conditions for the MHD simulations had to be performed in a different fashion as compared to my previous work, because the input data used there were found to be incorrect. Instead, the potential field source surface (Altschuler and Newkirk, 1969) coupled with the Schatten current sheet model (Schatten, 1971) was implemented, 80 because the employed so-called Wang-Sheeley-Arge (WSA; Arge and Pizzo, 2000; Arge et al., 2003) model for deriving the resulting solar wind speeds is tuned to the usage of GONG magnetograms in these potential field approaches. For the analysis of the resulting coronal field topology I developed a sophisticated field line tracing algorithm, which was also necessary to derive the resulting solar wind speeds. The inner boundary conditions for the remaining MHD quantities were set via empirical formulas that contain a number of tunable parameters, for which a best fit to observations was performed. My results showed good agreement with spacecraft data, allowing to expect insights for energetic particle propagation from the ongoing SDE modeling. For the latter it is desirable to have the MHD results extending to large distances, i.e., to the termination shock in order to remain in the supersonic part of the heliosphere. Such simulations require not only huge computer resources, but it is also mandatory to apply time-dependent boundary conditions to account for changes in the corona during the long propagation time of about a year that the solar wind needs to reach the termination shock. Although first steps in these directions have been taken by applying such boundary conditions, it will be necessary to apply for computing time at supercomputing facilities to keep the required simulation time within reasonable limits. This will be particularly necessary for a coupling with outer-heliospheric models to describe the interaction with the interstellar medium (ISM). While the SDE modeling in our group focuses on Galactic Cosmic Rays and Jovian electrons, similar studies for solar energetic particles are of interest as well. Such an investigation, for which the current results from my MHD simulations are principally suited at least in the supersonic solar wind, is planned in collaboration with W. Dröge at the University of Würzburg. However, a coronal model describing not only the magnetic field (as provided by my potential field solutions), but also the evolving solar wind would be more appropriate, which could be achieved with an MHD model of the solar corona. This could be a future direction of work that would provide a more self-consistent treatment of the corona in comparison with the empirically derived solar wind based on the coronal potential field solutions used in this thesis. Another possibility of exploiting my results is an investigation of space weather effects, i.e. impacts on the Earth’s magnetosphere and man-made objects. The employed transport coefficients in the SDE modeling depend not only on the large-scale MHD quantities, but also on the magnetic fluctuations. These cannot be resolved in large-scale MHD simulations, but there are models describing the evolution of integral turbulence quantities such as the total turbulent energy density and the ratio between inward and outward propagating modes. In order to eventually be able to provide these quantities for an improved SDE modeling, a turbulence transport model was implemented alongside the large-scale MHD equations in the Cronos framework, which is the topic of Chapter 4 based on my publication Implementing turbulence transport in the CRONOS framework and application to 81 CHAPTER 7. SUMMARY AND OUTLOOK the propagation of CMEs (Wiengarten et al., 2015). I started again by validating the implementation by comparing with previous work by Usmanov et al. (2011). During this comparison, some discrepancies were found that – in discussion with A. Usmanov – could be traced back to errors made by these authors. A corrected model was then used, which was further extended by removing the constraint of being only applicable in regions of highly super-Alfvénic solar wind (U VA ), while I followed the goal to have a model that can eventually be coupled to my previous WSA-driven MHD simulations that start in regions where U ≥ VA . This was achieved by appropriately simplifying a more general model of turbulence transport by Zank et al. (2012), which in the form used here describes the evolution of turbulent energy density, cross-helicity, and correlation lengths for low-frequency turbulence, as observed in the solar wind with increasing distance from the Sun. These integral quantities characterize the large-scale behavior of turbulence and wave modes, which is a useful approach to circumvent the difficulties that would arise for a complete description of the fluctuations, which is on much smaller scales than typical structures in the heliosphere. Particularly for numerical studies such an approach is desirable because of the reduced required computer resources. My resulting model was then applied to the propagation of CMEs, for which I found that, on the one hand, turbulence does not act back on the large-scale quantities, but, on the other hand, that the CME is a strong driver of turbulence. Therefore, such a model is not required for studies focusing on large-scale CME quantities alone, but it does provide another interesting scenario for a related energetic particle propagation study. There is only a very limited number of observations of turbulence associated with CMEs, however, my results are in agreement with estimates provided by Subramanian et al. (2009). Chapter 5 describes the subsequent implementation of the complete model of Zank et al. (2012), i.e. with evolution equations for the energy difference between kinetic and magnetic fluctuations, as well as respective correlation lengths. Thus far, this had been solved neither in 3D nor self-consistently coupled with the MHD equations. A simplified implementation in the ecliptic plane only and neglecting some of the involved terms was done by Adhikari et al. (2015), whose results served for validation purposes. Afterwards, the complete model was solved, which gave some yet not fully understood results. Consequently, this work is unpublished as of yet, but ongoing. Mandatory next steps are a thorough comparison with spacecraft data and an analysis of the effect of the additional terms, possibly in collaboration with the initiators of these models. There is also criticism to these kind of models taking into account 2D turbulence only, because especially the pickup ion isotropization in the outer heliosphere rather drives Alfvén waves (slab turbulence). Models of turbulent transport usually only consider one or the other, but there are also approaches for combined theories, such as the two-component model by Oughton et al. (2011), which will also be introduced into the Cronos code and investigated in the near future. Furthermore, the model solving the coupled Reynolds-averaged MHD and turbu82 lence transport equations can be fed with inner-boundary conditions derived via the WSA model to advance from the still idealized solar wind conditions applied here to actual conditions during specific periods of time. To do this it will be necessary to devise empirical formulas for the turbulence values at the inner boundary in a similar fashion as in the current WSA setup, i.e. depending on coronal conditions quantified by, e.g., the flux tube expansion factor. Another line of work in our team is the modeling of the large-scale heliosphere and astrospheres, i.e. the interaction between the solar/stellar wind and the ISM. An analytic model for the interstellar magnetic field in the vicinity of the heliopause was recently presented by Röken et al. (2015), in which it was also compared to numerical results obtained with the Cronos code by J. Kleimann. This numerical setup was altered by K. Scherer and myself for the application to the astrosphere of the O-star λ Cephei, which, due to the large terminal velocities of the stellar wind, is a cavity of several parsecs in diameter and is thus able to modulate high-energy cosmic rays that pass through it. This was the topic of Chapter 6 based on the paper Cosmic rays in astrospheres (Scherer et al., 2015). The required computer resources for the simulations performed here are already enormous, and a migration to supercomputing facilities will have to be undertaken, particularly for the planned extensions of the model to incorporate stellar and interstellar magnetic fields, and even a multi-fluid approach. In conclusion, the work described in this thesis provides the essential extension to the MHD code Cronos for an adequate modeling of the super-Alfvénic solar wind by basing the simulations on observationally derived input data and including a description for the evolution of turbulence in the heliosphere, which is particularly useful for related energetic particle propagation studies. 83 CHAPTER 7. SUMMARY AND OUTLOOK 84 Appendix A Non-equidistantly spaced grids Refining or coarsening the grid underlying a simulation domain is a desirable feature because it can better resolve interesting features of small scales on the one hand, while on the other hand the computational costs can be significantly reduced by using lower resolutions for areas of low variability. A powerful technique to achieve this goal is the adaptive mesh refinement (AMR) algorithm (Berger and Colella, 1989) that refines or coarsens a grid during run-time according to some prescribable criterion in the simulation data (e.g. gradients in density). Implementing AMR is a challenging task, especially for non-Cartesian geometries and highly-parallelized systems. AMR is currently not implemented in the Cronos framework, but may be addressed in a future release. Alternatively, one can consider a grid that, although it is not adapting itself during runtime, can be predefined to be finer or coarser in respective regions depending on the expected scale of changes in the simulation. Currently ongoing is the implementation of so-called logically rectangular grids (Calhoun et al., 2008) into the Cronos framework (B. Krebl, priv. comm.). A simpler approach that has found its way into the Cronos framework already is to prescribe the spacing between neighboring grid points for a given direction not to be constant but varying according to some function of the user’s choice. In the two following sections a general overview of these non-equidistantly spaced grids (NESGs) will be given first, followed by considerations important for the case subject to this thesis, namely spherical coordinates and the application to the solar wind. A.1 General overview Given a certain direction p ∈ {x, y, z}, {ρ, ϕ, z}, {r, ϑ, ϕ} of the chosen coordinate system, the location of a grid point in that direction depending on that grid points integer index i (0 ≤ i ≤ Np − 1), where Np specifies the total number of grid points in that direction, can be calculated via p(i) = pbeg + f (ξ)(pend − pbeg ) . 85 (A.1) APPENDIX A. NON-EQUIDISTANTLY SPACED GRIDS The respective beginning and end of the computational domain in that direction are denoted pbeg and pend , while for convenience ξ = i/(Np − 1) is introduced. It is obviously mandatory for the grid-generating function f to fulfill f (0) = 0 and f (1) = 1 in order to get p(0) = pbeg and p(1) = pend . Furthermore, f (ξ) should be a monotonically increasing function so that 0 ≤ f ≤ 1. Otherwise grid points could on the one hand be located outside the desired grid extent, i.e. beyond pend or below pbeg , or on the other hand a subsequent grid point i2 > i1 could be located before the previous one (p(i2 ) < p(i1 )). The class of equidistantly spaced-grids can be realized by choosing flin (ξ) = ξ , (A.2) which clearly satisfies the above demands and results in a constant grid point distance pend − pbeg . (A.3) ∆p = p(i + 1) − p(i) = Np − 1 This is the standard case in the Cronos setup. Another pre-implemented choice employs an exponential approach via fexp (ξ) = (pend /pbeg )ξ − 1 (pend /pbeg ) − 1 (A.4) resulting in p(i) = pbeg (pend /pbeg )ξ , the latter clearly showing the exponential character. Simple algebra gives ! 1 pend Np −1 ∆p(i) = p(i) · −1 , (A.5) pbeg so that the cell size is linearly increasing in this case. Both the linear and the exponential realization are visualized in Figure A.1, in which a one-dimensional grid from pbeg = 1 to pend = 100 is resolved with Np,lin = 100 grid points, while much less points (see Equation (A.7)) are required for the exponential realization. The left panel shows the grid point location p(i) against the grid point index i, where the linear/exponential character of the grid-generating function is evident. The right panel visualizes the variations of grid point distance ∆p with grid point location p, which was shown to be constant/linear for the linear/exponential grid-generating function. The benefit of such adapted grids is the option to put resolution where it is needed, which results in reduced computational costs due to (i) fewer required grid points and (ii) a possibly enhanced global time-step which is directly proportional to the smallest cell-size in the simulation domain according to the CFL condition ∆t ∝ min{∆pj } (geometrical factors in curvilinear coordinate system have been neglected in this statement). 86 A.2. APPLICATION Figure A.1: Illustration of linear and exponential grid-generating functions. A.2 Application In the light of this work there are two aspects that make NESGs a useful tool: on the one hand, heliospheric structures become smeared out (global merged interaction regions) with increasing distance and a coarser grid is applicable, which also reduces the required computer resources. On the other hand the employed spherical grids can contain the polar axis as singularities resulting in increasingly small grid point distances in the azimuthal direction towards the poles1 , which is an effect that can also be significantly reduced by choosing an appropriate NESG2 . These two aspects can be addressed by NESGs for the radial coordinate r and the polar angle ϑ. Naturally the question arises whether there exists some suitable function for the azimuthal angle ϕ as well. Unfortunately, this is not the case, which can be easily understood by considering the evolution of a CIR (see, e.g., Figure 12 in Wiengarten et al. (2014)), which is obviously not constrained in azimuth but with increasing radial distance covers an ever larger interval, so that a refinement at some specific ϕ would not permanently coincide with the CIR. It should also be kept in mind that the grid-generating function can only depend on the respective directions index, and not on some other directions index, i.e. a refinement in ϕ depending on r or ϑ is not possible, since this would result in non-orthogonal grids. Therefore, for the purposes described in this work, only the r and ϑ directions can be used in the 1 Alternative approaches to overcoming this problem usually involve compositions of several (rotated) spherical grids, such as, e.g., the Yin-Yang grids (Laukert, 2008), the six-component grid (Feng et al., 2010), or the composite grid (Usmanov et al., 2012). A drawback of these implementations is the required increased overhead due to communication between the grids and respectively more difficult boundary conditions, as well as a loss of conservative properties at the interfaces. 2 In the first two publications of this thesis this problem was circumvented by leaving out the polar regions as the required features of NESGs and handling of coordinate singularities were not available at the time. 87 APPENDIX A. NON-EQUIDISTANTLY SPACED GRIDS context of NESGs, as will be shown in the subsequent sections. A.2.1 NESG for the radial direction A modified version of the exponential grid described in Section A.1 is suitable, for which it is also possible to alter the largest cell sizes towards the outer boundary by means of the parameter a > 0, where larger values result in exceedingly coarser grids, i.e. the respectively adapted linearly increasing grid point distance (Equation (A.5)) becomes steeper. The grid-generating function in this case reads fM ODexp (ξ) = (pend /pbeg )aξ − 1 . (pend /pbeg )a − 1 (A.6) It is often desirable to fix the smallest cell size ∆r|rbeg = r(i = 1) − r(i = 0) because the smallest radial cell size usually determines the global time step according to the CFL condition. For fixed a and ∆r|rbeg the number of grid points becomes Nr = a ln rend rbeg −1 ∆r|rbeg [(rend /rbeg )a − 1] , +1 ln rend − rbeg (A.7) which has to be rounded to integer values, so that the chosen ∆r|rbeg will be fulfilled only approximately. Of course a can also be respectively adapted for fixed ∆r|rbeg . Another important feature that such a grid should provide is the possibility to solve the grid-generating function uniquely for the index i. It is easy to show that this is possible for the grids presented in the general overview above, as well as for the grid presented here, yielding a −1 r(i) − rbeg N −1 rend rend ln i(r) = − 1 + 1 ln . (A.8) a rend − rbeg rbeg rbeg This feature is important for the subsequent use of the simulation results in different codes, which need values at intermediate locations, for which the desired quantities have to be interpolated from the stored values at neighboring grid points. To give an impression of the achievable savings as compared to a linear grid a simulation box is considered as used in Wiengarten et al. (2015) of a radial extent r ∈ [0.3, 100] AU and a desired smallest cell size at the inner boundary of ∆r = 10R . A linear grid would require Nr,lin ≈ 2144, while with a = 0.55 the exponential grid comprises of Nr,modEXP = 300 with a largest cell size at the outer boundary of roughly 1 AU. A.2.2 NESG for the polar angle Instead of leaving out the polar singularity, the idea is to have two cones, each consisting of a large cell at the two polar caps, while maintaining a linear grid for 88 A.2. APPLICATION the remaining cells. The cell centers of the larger polar cells are then separated sufficiently from the polar axis to avoid the undesired small azimuthal cell diameter that would lead to small time steps. A further advantage as compared to leaving out the polar regions is the application of more realistic boundary conditions: while formerly either extrapolating or reflecting boundary conditions had to be applied at the polar angle’s outer boundaries (see Section 1.2.2 and Wiengarten (2011)), the formal inclusion of the polar axis uses boundary conditions via the adjacent, azimuthally opposite (shifted by π) cells of the simulation domain as ghost-cells, which is self-consistent. Figure A.2: Same as Figure A.1, but for an NESG for the polar angle ϑ ∈ [0◦ , 180◦ ] with polar cell sizes G = 18◦ and an otherwise constant cell size of 2◦ . The grid-generating function can be derived as follows: first the relative size of the large polar cells G (for gap) with respect to the extent of the grid in that direction is introduced. Then, the ansatz for a linear function f (ξ) = mξ + b is taken, where the slope for the linear part of the grid follows as m= 1 − 2G 1 − 2G ∆f = = Nlin ∆ξ (Nlin − 1)/Nlin − 1/Nlin Nlin − 2 (A.9) where Nlin is the number of grid points that would be required to have the whole polar angle covered with a linear grid of the desired constant cell size, e.g. for a cell size ∆ϑ = 1◦ , Nlin = π/1◦ = 180. The first and last grid point are set manually by requiring them to be separated by G from the polar axis, i.e. f (1/Nlin ) = G and f (1 − 1/Nlin ) = 1 − G, while the grid generating function for the part of constant cell size is 1 − 2G Nlin G − 1 f (ξ) = Nlin ξ+ . (A.10) Nlin − 2 Nlin − 2 The actually required number of grid points then is N = hNlin + 2 with h = f (Nlin − 1/Nlin ) − f (1/Nlin ) = 1 − 2G. 89 APPENDIX A. NON-EQUIDISTANTLY SPACED GRIDS First tests showed that the savings due to the reduced number of grid points can be about 15% in computing time, while the time step can be significantly enhanced up to 500% for a typical value of G = 0.1π and Nϕ = 180. A.3 Implementation The implementation of the new grid-generating functions described above is realized in Cronos by additional classes as described in the Cronos Manual (J. Kleimann, priv. comm.). The respective code fragment can be found below, where especially the mirrored ghost-cells for the NESG for the polar angle are of interest. /*** USER SPECIFIED GRID LAYOUTS HERE ***/ /********************************************************************/ /***** EXP GRID ******/ class GridFunction_modEXP: public Interface_GridFunction { public: //class declaration with public functions GridFunction_modEXP(REAL, REAL, REAL); //constructor virtual REAL get_gridFunc(REAL); //get-function private://private variables REAL domain_begin, domain_end, domain_len; //corresponding to the p’s in Eq.(A.1) REAL a; //parameter a from Section A.2.1 }; GridFunction_modEXP::GridFunction_modEXP(REAL domain_begin, REAL domain_end, REAL domain_len) { //constructor with transfer parameters provided by the pre-defined grid-boundaries //and a is read-in from the User accessible "cat-File" this->domain_begin = domain_begin; this->domain_end = domain_end; this->domain_len = domain_len; this->a = value((char*)"stretch_expGrid"); } // function gives 0 at ratio=0 and 1 at ratio=1 REAL GridFunction_modEXP::get_gridFunc(REAL ratio) { //Function return value according to Eq. (A.6) return (pow(domain_end/domain_begin,a*ratio)-1.0)/(pow(domain_end/domain_begin,a)-1.0); } /********************************************************************/ /***** LIN_POLES GRID with one big cell at poles ******/ class GridFunction_LIN_POLES: public Interface_GridFunction { public: //class declaration with public functions GridFunction_LIN_POLES(REAL, REAL, REAL); virtual REAL get_gridFunc(REAL); private: //private variables REAL domain_begin, domain_end, domain_len; //corresponding to the p’s in Eq.(A.1) //parameters as in Section A.2.2: REAL gap; //(=G) REAL N; //(=N_lin) }; GridFunction_LIN_POLES::GridFunction_LIN_POLES(REAL domain_begin, REAL domain_end, REAL domain_len) { //constructor with transfer parameters provided by the pre-defined grid-boundaries //while gap and N are read-in from the User accessible "cat-File" this->domain_begin = domain_begin; this->domain_end = domain_end; 90 A.3. IMPLEMENTATION this->domain_len = domain_len; this->gap = value((char*)"gap_poleGrid"); /************** size of large polar cell with respect to grid extent e.g. gap=0.1 refers to a polar cell size of 0.1*180.0[deg] = 18[deg] ***************/ this->N = value((char*)"Nlin_poleGrid"); } // function gives 0 at ratio=0 and 1 at ratio=1 REAL GridFunction_LIN_POLES::get_gridFunc(REAL ratio) { int adapted = 0;//switch for adapting mapped (opposite/shifted by pi) polar cells if (domain_begin < 0.01) if (ratio < 0){ //top direction ratio = -ratio; //take mirrored (positive) cells value, but return negative f below adapted = 1; } if(domain_end > 3.141) if (ratio > 1){//bottom direction ratio = 1.0 - (ratio-1.0); //mirror at 1 adapted = 2; } // Taking a linear function leaving out the polar cells gives Eq. (A.10): REAL f=N*(1.-2.*gap)/(N-2.)*ratio + (N*gap-1.)/(N-2.); //but has to be set manually to f(0)=0 and f(1)=1 if (ratio == 0) f=0; if (ratio == 1) f=1; //and mapped polar cells are mirrored: if (adapted == 1) f=-f; if (adapted == 2) f=1.0+(1.0-f); return f; } 91 APPENDIX A. NON-EQUIDISTANTLY SPACED GRIDS 92 Bibliography R. Abbasi, Y. Abdou, T. Abu-Zayyad, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers, M. M. Allen, D. Altmann, K. Andeen, and et al. Observation of Anisotropy in the Galactic Cosmic-Ray Arrival Directions at 400 TeV with IceCube. The Astrophysical Journal, 746:33, February 2012. doi: 10.1088/0004637X/746/1/33. L. Adhikari, G. P. Zank, R. Bruno, D. Telloni, P. Hunana, A. Dosch, R. Marino, and Q. Hu. The Transport of Low-frequency Turbulence in Astrophysical Flows. II. Solutions for the Super-Alfvénic Solar Wind. The Astrophysical Journal, 805: 63, May 2015. doi: 10.1088/0004-637X/805/1/63. S.-I. Akasofu, C. Fry, and K. Hakamada. Solar wind disturbances caused by solar flares - Equatorial plane. Planetary and Space Science, 31:1435–1458, December 1983. doi: 10.1016/0032-0633(83)90018-1. M. D. Altschuler and G. Newkirk. Magnetic Fields and the Structure of the Solar Corona. I: Methods of Calculating Coronal Fields. Solar Physics, 9:131–149, September 1969. doi: 10.1007/BF00145734. C. N. Arge and V. J. Pizzo. Improvement in the prediction of solar wind conditions using near-real time solar magnetic field updates. Journal of Geophysical Research, 105:10465–10480, 2000. doi: 10.1029/1999JA900262. C. N. Arge, D. Odstrcil, V. J. Pizzo, and L. R. Mayer. Improved Method for Specifying Solar Wind Speed Near the Sun. In M. Velli, R. Bruno, F. Malara, & B. Bucci, editor, Solar Wind Ten, volume 679 of American Institute of Physics Conference Series, pages 190–193, September 2003. doi: 10.1063/1.1618574. D. N. Baker. What is space weather? Advances in Space Research, 22:7–16, 1998. doi: 10.1016/S0273-1177(97)01095-8. A. Balogh, J. T. Gosling, J. R. Jokipii, R. Kallenbach, and H. Kunow. Corotating Interaction Regions. Proceedings. ISSI Workshop, Bern (Switzerland), 6 - 13 Jun 1998. Space Science Reviews, 89, 1999. P. K. Barker and J. M. Marlborough. Weber and Davis revisited - Mass losing rotating magnetic winds. The Astrophysical Journal, 254:297–300, March 1982. doi: 10.1086/159733. 93 BIBLIOGRAPHY C. S. Beals. On the nature of Wolf-Rayet emission. Monthly Notices of the Royal Astronomical Society, 90:202–212, December 1929. J. W. Belcher and L. Davis, Jr. Large-amplitude Alfvén waves in the interplanetary medium, 2. Journal of Geophysical Research, 76:3534, 1971. doi: 10.1029/JA076i016p03534. M. J. Berger and P. Colella. Local adaptive mesh refinement for shock hydrodynamics. Journal of Computational Physics, 82:64–84, May 1989. doi: 10.1016/00219991(89)90035-1. J. W. Bieber, W. Wanner, and W. H. Matthaeus. Dominant two-dimensional solar wind turbulence with implications for cosmic ray transport. J. Geophys. Res., 101:2511–2522, 1996. L. Biermann. Kometenschweife und solare Korpuskularstrahlung. Zeitschrift für Astrophysik, 29:274, 1951. T. J. Bogdan and B. C. Low. The three-dimensional structure of magnetostatic atmospheres. II - Modeling the large-scale corona. The Astrophysical Journal, 306:271–283, July 1986. doi: 10.1086/164341. F. Boulanger, P. Cox, and A. P. Jones. Course 7: Dust in the Interstellar Medium. In F. Casoli, J. Lequeux, and F. David, editors, Infrared Space Astronomy, Today and Tomorrow, page 251, 2000. J. U. Brackbill and D. C. Barnes. The effect of nonzero product of magnetic gradient and B on the numerical solution of the magnetohydrodynamic equations. Journal of Computational Physics, 35:426–430, May 1980. doi: 10.1016/00219991(80)90079-0. B. Breech, W. H. Matthaeus, J. Minnie, J. W. Bieber, S. Oughton, C. W. Smith, and P. A. Isenberg. Turbulence transport throughout the heliosphere. Journal of Geophysical Research, 113:8105, 2008. doi: 10.1029/2007JA012711. R. A. Burger and M. Hitge. The Effect of a Fisk-Type Heliospheric Magnetic Field on Cosmic-Ray Modulation. The Astrophysical Journal Letters, 617:L73–L76, 2004. doi: 10.1086/427076. R. A. Burger, T. P. J. Krüger, M. Hitge, and N. E. Engelbrecht. A Fisk-Parker Hybrid Heliospheric Magnetic Field With a Solar-Cycle Dependence. The Astrophysical Journal, 674:511–519, 2008. doi: 10.1086/525039. I. Büsching, A. Kopp, M. Pohl, R. Schlickeiser, C. Perrot, and I. Grenier. CosmicRay Propagation Properties for an Origin in Supernova Remnants. The Astrophysical Journal, 619:314–326, 2005. doi: 10.1086/426537. 94 BIBLIOGRAPHY D. A. Calhoun, C. Helzel, and R. J. Leveque. Logically Rectangular Grids and Finite Volume Methods for PDEs in Circular and Spherical Domains. SIAM Review, 50: 723–752, January 2008. doi: 10.1137/060664094. P. J. Cargill and J. A. Klimchuk. Nanoflare Heating of the Corona Revisited. The Astrophysical Journal, 605:911–920, April 2004. doi: 10.1086/382526. M. Carlsson, V. H. Hansteen, B. de Pontieu, S. McIntosh, T. D. Tarbell, D. Shine, S. Tsuneta, Y. Katsukawa, K. Ichimoto, Y. Suematsu, T. Shimizu, and S. Nagata. Can High Frequency Acoustic Waves Heat the Quiet Sun Chromosphere? Publications of the Astronomical Society of Japan, 59:663, November 2007. doi: 10.1093/pasj/59.sp3.S663. S. Chapman and H. Zirin. Notes on the Solar Corona and the Terrestrial Ionosphere. Smithsonian Contributions to Astrophysics, 2:1, 1957. O. Cohen, I. V. Sokolov, I. I. Roussev, C. N. Arge, W. B. Manchester, T. I. Gombosi, R. A. Frazin, H. Park, M. D. Butala, F. Kamalabadi, and M. Velli. A Semiempirical Magnetohydrodynamical Model of the Solar Wind. The Astrophysical Journal Letters, 654:L163–L166, January 2007. doi: 10.1086/511154. P. J. Coleman, Jr. Turbulence, Viscosity, and Dissipation in the Solar-Wind Plasma. The Astrophysical Journal, 153:371, August 1968. doi: 10.1086/149674. R. Courant, K. Friedrichs, and H. Lewy. Über die partiellen Differenzengleichungen der mathematischen Physik. Mathematische Annalen, 100:32–74, 1928. doi: 10.1007/BF01448839. S. R. Cranmer and A. A. van Ballegooijen. On the Generation, Propagation, and Reflection of Alfvén Waves from the Solar Photosphere to the Distant Heliosphere. The Astrophysical Journal Supplement Series, 156:265–293, February 2005. doi: 10.1086/426507. S. R. Cranmer, A. A. van Ballegooijen, and R. J. Edgar. Self-consistent Coronal Heating and Solar Wind Acceleration from Anisotropic Magnetohydrodynamic Turbulence. The Astrophysical Journal Supplement Series, 171:520–551, August 2007. doi: 10.1086/518001. A. C. Cummings, J. R. Cummings, R. A. Mewaldt, E. C. Stone, B. Blake, M. Fraenz, B. Klecker, D. Hovestadt, and W. R. Webber. Observations of anomalous cosmic rays in the heliosphere from the SAMPEX, Ulysses, Voyager, and Pioneer spacecraft. Advances in Space Research, 16:337, 1995. doi: 10.1016/02731177(95)00357-K. G. Dalakishvili, J. Kleimann, H. Fichtner, and S. Poedts. Magnetic clouds in the solar wind: a numerical assessment of analytical models. Astronomy and Astrophysics, 536:A100, December 2011. doi: 10.1051/0004-6361/201015468. 95 BIBLIOGRAPHY R. Davis, D. S. Harmer, and K. C. Hoffman. Search for Neutrinos from the Sun. Physical Review Letters, 20:1205–1209, May 1968. doi: 10.1103/PhysRevLett.20.1205. M. A. de Avillez, D. Breitschwerdt, A. Asgekar, and E. Spitoni. ISM simulations: an overview of models. Highlights of Astronomy, 16:606–608, March 2015. doi: 10.1017/S1743921314012393. T. Detman, Z. Smith, M. Dryer, C. D. Fry, C. N. Arge, and V. Pizzo. A hybrid heliospheric modeling system: Background solar wind. Journal of Geophysical Research, 111:7102, 2006. doi: 10.1029/2005JA011430. T. R. Detman, D. S. Intriligator, M. Dryer, W. Sun, C. S. Deehr, and J. Intriligator. The influence of pickup protons, from interstellar neutral hydrogen, on the propagation of interplanetary shocks from the Halloween 2003 solar events to ACE and Ulysses: A 3-D MHD modeling study. Journal of Geophysical Research (Space Physics), 116(A15):A03105, March 2011. doi: 10.1029/2010JA015803. J.-F. Donati, J. Morin, P. Petit, X. Delfosse, T. Forveille, M. Aurière, R. Cabanac, B. Dintrans, R. Fares, T. Gastine, M. M. Jardine, F. Lignières, F. Paletou, J. C. Ramirez Velez, and S. Théado. Large-scale magnetic topologies of early M dwarfs. Monthly Notices of the Royal Astronomical Society, 390:545–560, October 2008. doi: 10.1111/j.1365-2966.2008.13799.x. A. Dosch, L. Adhikari, and G. P. Zank. The transport of low-frequency turbulence in astrophysical flows: Correlation lengths. In G. P. Zank, J. Borovsky, R. Bruno, J. Cirtain, S. Cranmer, H. Elliott, J. Giacalone, W. Gonzalez, G. Li, E. Marsch, E. Moebius, N. Pogorelov, J. Spann, and O. Verkhoglyadova, editors, American Institute of Physics Conference Series, volume 1539 of American Institute of Physics Conference Series, pages 155–158, June 2013. doi: 10.1063/1.4811011. N. E. Engelbrecht and R. A. Burger. An Ab Initio Model for the Modulation of Galactic Cosmic-ray Electrons. The Astrophysical Journal, 779:158, December 2013. doi: 10.1088/0004-637X/779/2/158. X. Feng, L. Yang, C. Xiang, S. T. Wu, Y. Zhou, and D. Zhong. Three-dimensional Solar WIND Modeling from the Sun to Earth by a SIP-CESE MHD Model with a Six-component Grid. The Astrophysical Journal, 723:300–319, November 2010. doi: 10.1088/0004-637X/723/1/300. L. A. Fisk. Motion of the footpoints of heliospheric magnetic field lines at the Sun: Implications for recurrent energetic particle events at high heliographic latitudes. Journal of Geophysical Research, 101:15547–15553, 1996. L. A. Fisk and G. Gloeckler. On Whether or Not Voyager 1 has Crossed the Heliopause. The Astrophysical Journal, 789:41, July 2014. doi: 10.1088/0004637X/789/1/41. 96 BIBLIOGRAPHY P. R. Gazis, A. Barnes, J. D. Mihalov, and A. J. Lazarus. Solar wind velocity and temperature in the outer heliosphere. Journal of Geophysical Research, 99: 6561–6573, April 1994. doi: 10.1029/93JA03144. J. Giacalone, J. R. Jokipii, and J. Kóta. Particle Acceleration in Solar Wind Compression Regions. The Astrophysical Journal, 573:845–850, July 2002. doi: 10.1086/340660. R. Grauer, C. Marliani, and K. Germaschewski. Adaptive Mesh Refinement for Singular Solutions of the Incompressible Euler Equations. Physical Review Letters, 80:4177–4180, May 1998. doi: 10.1103/PhysRevLett.80.4177. M. A. T. Groenewegen, C. Waelkens, M. J. Barlow, F. Kerschbaum, P. Garcia-Lario, J. Cernicharo, et al. MESS (Mass-loss of Evolved StarS), a Herschel key program. Astronomy and Astrophysics, 526:A162, February 2011. doi: 10.1051/00046361/201015829. J. Grygorczuk, A. Czechowski, and S. Grzedzielski. Why are the Magnetic Field Directions Measured by Voyager 1 on Both Sides of the Heliopause so Similar? The Astrophysical Journal Letters, 789:L43, July 2014. doi: 10.1088/20418205/789/2/L43. J. V. Hollweg. On electron heat conduction in the solar wind. Journal of Geophysical Research, 79:3845, 1974. doi: 10.1029/JA079i025p03845. J. V. Hollweg. Collisionless electron heat conduction in the solar wind. Journal of Geophysical Research, 81:1649–1658, April 1976. doi: 10.1029/JA081i010p01649. T. S. Horbury and K. T. Osman. Multi-Spacecraft Turbulence Analysis Methods. ISSI Scientific Reports Series, 8:55–64, 2008. T. S. Horbury, R. T. Wicks, and C. H. K. Chen. Anisotropy in Space Plasma Turbulence: Solar Wind Observations. Space Science Reviews, 172:325–342, November 2012. doi: 10.1007/s11214-011-9821-9. Y. Q. Hu, S. R. Habbal, and X. Li. On the cascade process of Alfvén waves in the fast solar wind. Journal of Geophysical Research, 104:24819–24834, November 1999. doi: 10.1029/1999JA900340. S. Hubrig, M. Schöller, L. Fossati, T. Morel, N. Castro, L. M. Oskinova, N. Przybilla, S. S. Eikenberry, M.-F. Nieva, and N. Langer. B fields in OB stars (BOB): FORS 2 spectropolarimetric follow-up of the two rare rigidly rotating magnetosphere stars HD 23478 and HD 345439. Astronomy and Astrophysics, 578:L3, June 2015. doi: 10.1051/0004-6361/201526262. A. J. Hundhausen. Interplanetary Shock Waves and the Structure of Solar Wind Disturbances. NASA Special Publication, 308:393, 1972. 97 BIBLIOGRAPHY W.-H. Ip, A. Kopp, and J.-H. Hu. On the Star-Magnetosphere Interaction of Close-in Exoplanets. The Astrophysical Journal Letters, 602:L53–L56, 2004. doi: 10.1086/382274. P. A. Isenberg. Turbulence-driven Solar Wind Heating and Energization of Pickup Protons in the Outer Heliosphere. The Astrophysical Journal, 623:502–510, April 2005. doi: 10.1086/428609. V. V. Izmodenov, D. B. Alexashov, and M. S. Ruderman. Electron Thermal Conduction as a Possible Physical Mechanism to Make the Inner Heliosheath Thinner. The Astrophysical Journal Letters, 795:L7, November 2014. doi: 10.1088/20418205/795/1/L7. J. Jiang, R. Cameron, D. Schmitt, and M. Schüssler. Modeling the Sun’s Open Magnetic Flux and the Heliospheric Current Sheet. The Astrophysical Journal, 709:301–307, 2010. doi: 10.1088/0004-637X/709/1/301. J. Jiang, R. H. Cameron, and M. Schüssler. Effects of the Scatter in Sunspot Group Tilt Angles on the Large-scale Magnetic Field at the Solar Surface. The Astrophysical Journal, 791:5, August 2014. doi: 10.1088/0004-637X/791/1/5. J. R. Jokipii and J. Kota. The polar heliospheric magnetic field. Geophysical Research Letters, 16:1–4, 1989. M. L. Kaiser. The STEREO mission: an overview. Advances in Space Research, 36: 1483–1488, 2005. doi: 10.1016/j.asr.2004.12.066. R. Kissmann, J. Kleimann, H. Fichtner, and R. Grauer. Local turbulence simulations for the multiphase ISM. Monthly Notices of the Royal Astronomical Society, 391: 1577–1588, December 2008. doi: 10.1111/j.1365-2966.2008.13974.x. R. Kissmann, M. Flaig, and W. Kley. Accretion Disk Turbulence With a Detailed Thermodynamics. In N. V. Pogorelov, E. Audit, and G. P. Zank, editors, 5th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2010), volume 444 of Astronomical Society of the Pacific Conference Series, page 36, October 2011. J. Kleimann, A. Kopp, H. Fichtner, and R. Grauer. A novel code for numerical 3-D MHD studies of CME expansion. Annales Geophysicae, 27:989–1004, 2009. W. Kley. On the treatment of the Coriolis force in computational astrophysics. Astronomy and Astrophysics, 338:L37–L41, October 1998. A. Kopp, I. Büsching, R. D. Strauss, and M. S. Potgieter. A stochastic differential equation code for multidimensional Fokker-Planck type problems. Computer Physics Communications, 183:530–542, March 2012. doi: 10.1016/j.cpc.2011.11.014. 98 BIBLIOGRAPHY A. Kopp, I. Büsching, M. S. Potgieter, and R. D. Strauss. A stochastic approach to Galactic proton propagation: Influence of the spiral arm structure. New Astronomy, 30:32–37, July 2014. doi: 10.1016/j.newast.2014.01.006. A. G. Kosovichev. Advances in Global and Local Helioseismology: An Introductory Review. In J.-P. Rozelot and C. Neiner, editors, Lecture Notes in Physics, Berlin Springer Verlag, volume 832 of Lecture Notes in Physics, Berlin Springer Verlag, pages 3–642, 2011. doi: 10.1007/978-3-642-19928-81. M. W. Kunz, J. M. Stone, and X.-N. Bai. Pegasus: A new hybrid-kinetic particle-incell code for astrophysical plasma dynamics. Journal of Computational Physics, 259:154–174, February 2014. doi: 10.1016/j.jcp.2013.11.035. A. Kurganov and D. Levy. A Third-Order Semi-Discrete Central Scheme for Conservation Laws and Convection-Diffusion Equations. ArXiv Mathematics e-prints, February 2000. H. J. G. L. M. Lamers and J. P. Cassinelli. Introduction to Stellar Winds. Cambridge, UK: Cambridge University Press, ISBN 0521593980, 1999. K. R. Lang. Sun, Earth and Sky. Springer Science+Business Media, LLC, New York, NY, 2006. M. Laukert. MHD-Simulationen auf zusammengesetzten Gittern in sphärischer Geometrie. Diplomarbeit, Theoretische Physik IV, Ruhr Universität Bochum, 2008. J. A. le Roux and H. Fichtner. Global merged interaction regions, the heliospheric termination shock, and time-dependent cosmic ray modulation. Journal of Geophysical Research, 104:4709–4730, 1999. J. L. Linsky and B. E. Wood. Lyman-α observations of astrospheres. ASTRA Proceedings, 1:43–49, August 2014. doi: 10.5194/ap-1-43-2014. K. B. MacGregor and V. J. Pizzo. Mass loss from rotating magnetic stars - Weber and Davis re-revisited. The Astrophysical Journal, 267:340–343, April 1983. doi: 10.1086/160873. W. H. Matthaeus, S. Oughton, D. H. Pontius, Jr., and Y. Zhou. Evolution of energycontaining turbulent eddies in the solar wind. Journal of Geophysical Research, 99:19267, October 1994. doi: 10.1029/94JA01233. D. J. McComas, H. O. Funsten, S. A. Fuselier, W. S. Lewis, E. Möbius, and N. A. Schwadron. IBEX observations of heliospheric energetic neutral atoms: Current understanding and future directions. Geophysical Research Letters, 38:L18101, September 2011. doi: 10.1029/2011GL048763. 99 BIBLIOGRAPHY D. J. McComas, F. Allegrini, M. Bzowski, M. A. Dayeh, R. DeMajistre, H. O. Funsten, S. A. Fuselier, M. Gruntman, P. H. Janzen, M. A. Kubiak, H. Kucharek, E. Möbius, D. B. Reisenfeld, N. A. Schwadron, J. M. Sokól, and M. Tokumaru. IBEX: The First Five Years (2009-2013). The Astrophysical Journal Supplement Series, 213:20, August 2014. doi: 10.1088/0067-0049/213/2/20. N. Meyer-Vernet. Basics of the Solar Wind. Cambridge University Press, January 2007. M. S. Miesch, W. H. Matthaeus, A. Brandenburg, A. Petrosyan, A. Pouquet, C. Cambon, F. Jenko, D. Uzdensky, J. Stone, S. Tobias, J. Toomre, and M. Velli. Large-Eddy Simulations of Magnetohydrodynamic Turbulence in Astrophysics and Space Physics. ArXiv e-prints, May 2015. A. Mignone, M. Flock, M. Stute, S. M. Kolb, and G. Muscianisi. A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code. Astronomy and Astrophysics, 545:A152, September 2012. doi: 10.1051/0004-6361/201219557. R. Mitalas and K. R. Sills. On the photon diffusion time scale for the sun. The Astrophysical Journal, 401:759, December 1992. doi: 10.1086/172103. W.-C. Müller and R. Grappin. Spectral Energy Dynamics in Magnetohydrodynamic Turbulence. Physical Review Letters, 95(11):114502, September 2005. doi: 10.1103/PhysRevLett.95.114502. Y. Nakagawa. Evolution of Magnetic Field and Atmospheric Responses - Part Two - Formulation of Proper Boundary Equations. The Astrophysical Journal, 247: 719, July 1981. doi: 10.1086/159083. D. Odstrčil, P. Riley, and X. P. Zhao. Numerical simulation of the 12 May 1997 interplanetary CME event. Journal of Geophysical Research (Space Physics), 109: A02116, February 2004. doi: 10.1029/2003JA010135. M. Opher and J. F. Drake. On the Rotation of the Magnetic Field Across the Heliopause. The Astrophysical Journal Letters, 778:L26, December 2013. doi: 10.1088/2041-8205/778/2/L26. M. Opher, J. F. Drake, B. Zieger, and T. I. Gombosi. Magnetized Jets Driven By the Sun: the Structure of the Heliosphere Revisited. The Astrophysical Journal Letters, 800:L28, February 2015. doi: 10.1088/2041-8205/800/2/L28. S. Oughton, W. H. Matthaeus, C. W. Smith, B. Breech, and P. A. Isenberg. Transport of solar wind fluctuations: A two-component model. Journal of Geophysical Research (Space Physics), 116:A08105, August 2011. doi: 10.1029/2010JA016365. E. N. Parker. Dynamics of the Interplanetary Gas and Magnetic Fields. The Astrophysical Journal, 128:664, 1958. 100 BIBLIOGRAPHY E. N. Parker. The Hydrodynamic Theory of Solar Corpuscular Radiation and Stellar Winds. The Astrophysical Journal, 132:821, November 1960. doi: 10.1086/146985. E. N. Parker. Topological Dissipation and the Small-Scale Fields in Turbulent Gases. The Astrophysical Journal, 174:499, June 1972. doi: 10.1086/151512. V. J. Pizzo. A three-dimensional model of corotating streams in the solar wind. III Magnetohydrodynamic streams. Journal of Geophysical Research, 87:4374–4394, June 1982. doi: 10.1029/JA087iA06p04374. N. V. Pogorelov, S. N. Borovikov, J. Heerikhuisen, T. K. Kim, and G. P. Zank. Time-dependent Processes in the Sheath Between the Heliospheric Termination Shock and the Heliopause. In N. V. Pogorelov, E. Audit, and G. P. Zank, editors, 8th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2013), volume 488 of Astronomical Society of the Pacific Conference Series, page 167, September 2014. M. Potgieter. Solar Modulation of Cosmic Rays. Living Reviews in Solar Physics, 10:3, June 2013. doi: 10.12942/lrsp-2013-3. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical recipes – The art of scientific computing. Cambridge University Press, 2007. S. Preusse, A. Kopp, J. Büchner, and U. Motschmann. Stellar wind regimes of close-in extrasolar planets. Astronomy and Astrophysics, 434:1191–1200, 2005. doi: 10.1051/0004-6361:20041680. S. Preusse, A. Kopp, J. Büchner, and U. Motschmann. MHD simulation scenarios of the stellar wind interaction with Hot Jupiter magnetospheres. Planetary and Space Science, 55:589–597, April 2007. doi: 10.1016/j.pss.2006.04.037. E. R. Priest. Magnetohydrodynamics of the sun. Cambridge University Press, New York, NY, Jan 2014. URL https://cds.cern.ch/record/1640811. K. Reitberger, R. Kissmann, A. Reimer, and O. Reimer. Simulating Threedimensional Nonthermal High-energy Photon Emission in Colliding-wind Binaries. The Astrophysical Journal, 789:87, July 2014a. doi: 10.1088/0004637X/789/1/87. K. Reitberger, R. Kissmann, A. Reimer, O. Reimer, and G. Dubus. High-energy Particle Transport in Three-dimensional Hydrodynamic Models of Colliding-wind Binaries. The Astrophysical Journal, 782:96, February 2014b. doi: 10.1088/0004637X/782/2/96. P. Riley, J. A. Linker, Z. Mikić, R. Lionello, S. A. Ledvina, and J. G. Luhmann. A Comparison between Global Solar Magnetohydrodynamic and Potential Field Source Surface Model Results. The Astrophysical Journal, 653:1510–1516, December 2006. doi: 10.1086/508565. 101 BIBLIOGRAPHY P. Riley, R. Lionello, J. A. Linker, Z. Mikic, J. Luhmann, and J. Wijaya. Global MHD Modeling of the Solar Corona and Inner Heliosphere for the Whole Heliosphere Interval. Solar Physics, 274:361–377, December 2011. doi: 10.1007/s11207010-9698-x. C. Röken, J. Kleimann, and H. Fichtner. An Exact Analytical Solution for the Interstellar Magnetic Field in the Vicinity of the Heliosphere. The Astrophysical Journal, 805:173, June 2015. doi: 10.1088/0004-637X/805/2/173. T. Sakurai. Magnetic stellar winds - A 2-D generalization of the Weber-Davis model. Astronomy and Astrophysics, 152:121–129, November 1985. K. H. Schatten. Current sheet magnetic model for the solar corona. Cosmic Electrodynamics, 2:232–245, 1971. K. H. Schatten, J. M. Wilcox, and N. F. Ness. A model of interplanetary and coronal magnetic fields. Solar Physics, 6:442–455, March 1969. doi: 10.1007/BF00146478. K. Scherer and H. Fichtner. The Return of the Bow Shock. The Astrophysical Journal, 782:25, February 2014. doi: 10.1088/0004-637X/782/1/25. K. Scherer, H. Fichtner, T. Borrmann, J. Beer, L. Desorgher, E. Flückiger, H. J. Fahr, S. E. S. Ferreira, U. Langner, M. S. Potgieter, B. Heber, J. Masarik, N. Shaviv, and J. Veizer. Interstellar-Terrestrial Relations, The Variable Cosmic Environments, the Dynamic Heliosphere, and Their Imprints on Terrestrial Archives. Space Science Reviews, 127:327–465, 2006. doi: 10.1007/s11214-006-9126-6. K. Scherer, H. Fichtner, B. Heber, S. E. S. Ferreira, and M. S. Potgieter. Cosmic ray flux at the Earth in a variable heliosphere. Advances in Space Research, 41: 1171–1176, 2008. doi: 10.1016/j.asr.2007.03.016. K. Scherer, H. Fichtner, F. Effenberger, R. A. Burger, and T. Wiengarten. Comparison of different analytic heliospheric magnetic field configurations and their significance for the particle injection at the termination shock. Astronomy and Astrophysics, 521:A1, 2010. doi: 10.1051/0004-6361/200913638. K. Scherer, H. Fichtner, H.-J. Fahr, M. Bzowski, and S. E. S. Ferreira. Ionization rates in the heliosheath and in astrosheaths. Spatial dependence and dynamical relevance. Astronomy and Astrophysics, 563:A69, March 2014. doi: 10.1051/00046361/201321151. K. Scherer, A. van der Schyff, D. J. Bomans, S. E. S. Ferreira, H. Fichtner, J. Kleimann, R. D. Strauss, K. Weis, T. Wiengarten, and T. Wodzinski. Cosmic rays in astrospheres. Astronomy and Astrophysics, 576:A97, April 2015. doi: 10.1051/0004-6361/201425091. 102 BIBLIOGRAPHY N. A. Schwadron. An explanation for strongly underwound magnetic field in co-rotating rarefaction regions and its relationship to footpoint motion on the the sun. Geophysical Research Letters, 29(14):140000–1, 2002. doi: 10.1029/2002GL015028. N. A. Schwadron, M. Owens, and N. U. Crooker. The Heliospheric Magnetic Field over the Hale Cycle. ASTRA, 4:19–26, 2008. K. Schwarzschild. Gesammelte Werke Vol.1, Vol.2, Vol.3. Springer, 1992. R. Schwenn and E. Marsch. Physics of the Inner Heliosphere I. Large-Scale Phenomena. Springer Verlag, 1990. B. M. Shergelashvili and H. Fichtner. On the Low-frequency Boundary of Sungenerated Magnetohydrodynamic Turbulence in the Slow Solar Wind. The Astrophysical Journal, 752:142, June 2012. doi: 10.1088/0004-637X/752/2/142. C. W. Smith and J. W. Bieber. Solar cycle variation of the interplanetary magnetic field spiral. The Astrophysical Journal, 370:435–441, 1991. H. B. Snodgrass and R. K. Ulrich. Rotation of Doppler features in the solar photosphere. The Astrophysical Journal, 351:309–316, March 1990. doi: 10.1086/168467. C. W. Snyder and M. Neugebauer. Direct Observations of the Solar Wind by the Mariner II spacecraft. International Cosmic Ray Conference, 1:210, 1963. E. A. Spiegel and J.-P. Zahn. The solar tachocline. Astronomy and Astrophysics, 265:106–114, November 1992. O. Sternal, A. Burger, B. Heber, H. Fichtner, and P. Dunzlaff. The Diffusion Tensor of Energetic Particles in Different HMF Configurations. Proc. 30th Int. Cosmic Ray Conf., Paper SH3.1-885, 2007. M. Stix. The Sun: an introduction. Springer-Verlag, 2004. P. Subramanian, H. M. Antia, S. R. Dugad, U. D. Goswami, S. K. Gupta, Y. Hayashi, N. Ito, S. Kawakami, H. Kojima, P. K. Mohanty, P. K. Nayak, T. Nonaka, A. Oshima, K. Sivaprasad, H. Tanaka, S. C. Tonwar, and Grapes3 Collaboration. Forbush decreases and turbulence levels at coronal mass ejection fronts. Astronomy and Astrophysics, 494:1107–1118, February 2009. doi: 10.1051/0004-6361:200809551. C.-Y. Tu and E. Marsch. MHD structures, waves and turbulence in the solar wind: Observations and theories. Space Science Reviews, 73:1–210, July 1995. doi: 10.1007/BF00748891. 103 BIBLIOGRAPHY C.-Y. Tu, Z.-Y. Pu, and F.-S. Wei. The power spectrum of interplanetary Alfvenic fluctuations Derivation of the governing equation and its solution. Journal of Geophysical Research, 89:9695–9702, November 1984. doi: 10.1029/JA089iA11p09695. A. V. Usmanov and M. L. Goldstein. A tilted-dipole MHD model of the solar corona and solar wind. Journal of Geophysical Research (Space Physics), 108: 1354, September 2003. doi: 10.1029/2002JA009777. A. V. Usmanov, W. H. Matthaeus, B. A. Breech, and M. L. Goldstein. Solar Wind Modeling with Turbulence Transport and Heating. The Astrophysical Journal, 727:84, February 2011. doi: 10.1088/0004-637X/727/2/84. A. V. Usmanov, M. L. Goldstein, and W. H. Matthaeus. Three-dimensional Magnetohydrodynamic Modeling of the Solar Wind Including Pickup Protons and Turbulence Transport. The Astrophysical Journal, 754:40, July 2012. doi: 10.1088/0004-637X/754/1/40. A. V. Usmanov, M. L. Goldstein, and W. H. Matthaeus. Three-fluid, Threedimensional Magnetohydrodynamic Solar Wind Model with Eddy Viscosity and Turbulent Resistivity. The Astrophysical Journal, 788:43, June 2014. doi: 10.1088/0004-637X/788/1/43. R. Vainio, T. Laitinen, and H. Fichtner. A simple analytical expression for the power spectrum of cascading Alfvén waves in the solar wind. Astronomy and Astrophysics, 407:713–723, August 2003. doi: 10.1051/0004-6361:20030914. B. van der Holst, I. V. Sokolov, X. Meng, M. Jin, W. B. Manchester, IV, G. Tóth, and T. I. Gombosi. Alfvén Wave Solar Model (AWSoM): Coronal Heating. The Astrophysical Journal, 782:81, February 2014. doi: 10.1088/0004-637X/782/2/81. A. A. Vidotto. Incorporating magnetic field observations in wind models of low-mass stars. ASTRA Proceedings, 1:19–22, June 2014. doi: 10.5194/ap-1-19-2014. A. A. Vidotto, M. Jardine, M. Opher, J. F. Donati, and T. I. Gombosi. Powerful winds from low-mass stars: V374 Peg. Monthly Notices of the Royal Astronomical Society, 412:351–362, 2011. doi: 10.1111/j.1365-2966.2010.17908.x. A. A. Vidotto, R. Fares, M. Jardine, C. Moutou, and J.-F. Donati. On the environment surrounding close-in exoplanets. Monthly Notices of the Royal Astronomical Society, 449:4117–4130, June 2015. doi: 10.1093/mnras/stv618. Y.-M. Wang and N. R. Sheeley, Jr. Solar wind speed and coronal flux-tube expansion. The Astrophysical Journal, 355:726–732, June 1990. doi: 10.1086/168805. E. J. Weber and L. Davis, Jr. The Angular Momentum of the Solar Wind. The Astrophysical Journal, 148:217–227, 1967. doi: 10.1086/149138. 104 BIBLIOGRAPHY T. Wiegelmann, J. K. Thalmann, and S. K. Solanki. The magnetic field in the solar atmosphere. The Astronomy and Astrophysics Review, 22:78, November 2014. doi: 10.1007/s00159-014-0078-7. T. Wiengarten. Vergleich verschiedener Magnetfeldkonfigurationen im Bereich des heliosphärischen Terminationsschocks. Bachelor’s thesis, Theoretische Physik IV, Ruhr Universität Bochum, 2009. T. Wiengarten. MHD simulation of the inner-heliospheric magnetic field. Master’s thesis, Theoretische Physik IV, Ruhr Universität Bochum, 2011. T. Wiengarten, J. Kleimann, H. Fichtner, R. Cameron, J. Jiang, R. Kissmann, and K. Scherer. MHD simulation of the inner-heliospheric magnetic field. Journal of Geophysical Research (Space Physics), 118:29–44, January 2013. doi: 10.1029/2012JA018089. T. Wiengarten, J. Kleimann, H. Fichtner, P. Kühl, A. Kopp, B. Heber, and R. Kissmann. Cosmic Ray Transport in Heliospheric Magnetic Structures. I. Modeling Background Solar Wind Using the CRONOS Magnetohydrodynamic Code. The Astrophysical Journal, 788:80, June 2014. doi: 10.1088/0004-637X/788/1/80. T. Wiengarten, H. Fichtner, J. Kleimann, and R. Kissmann. Implementing Turbulence Transport in the CRONOS Framework and Application to the Propagation of CMEs. The Astrophysical Journal, 805:155, June 2015. doi: 10.1088/0004637X/805/2/155. B. E. Wood, J. L. Linsky, and G. P. Zank. Heliospheric, Astrospheric, and Interstellar Lyα Absorption toward 36 Ophiuchi. The Astrophysical Journal, 537:304–311, 2000. doi: 10.1086/309026. N. Yokoi, K. Higashimori, and M. Hoshino. Transport enhancement and suppression in turbulent magnetic reconnection: A self-consistent turbulence modela). Physics of Plasmas, 20(12):122310, December 2013. doi: 10.1063/1.4851976. G. Yu. The interstellar wake of the solar wind. The Astrophysical Journal, 194: 187–202, November 1974. doi: 10.1086/153235. G. P. Zank. Interaction of the solar wind with the local interstellar medium: a theoretical perspective. Space Science Reviews, 89:413–688, 1999. G. P. Zank, W. H. Matthaeus, and C. W. Smith. Evolution of turbulent magnetic fluctuation power with heliospheric distance. Journal of Geophysical Research, 101:17093–17107, 1996. G. P. Zank, A. Dosch, P. Hunana, V. Florinski, W. H. Matthaeus, and G. M. Webb. The Transport of Low-frequency Turbulence in Astrophysical Flows. I. Governing Equations. The Astrophysical Journal, 745:35, January 2012. doi: 10.1088/0004637X/745/1/35. 105 X. P. Zhao and J. T. Hoeksema. The Magnetic Field at the Inner Boundary of the Heliosphere Around Solar Minimum. Solar Physics, 266:379–390, 2010. doi: 10.1007/s11207-010-9618-0. Y. Zhou and W. H. Matthaeus. Transport and turbulence modeling of solar wind fluctuations. Journal of Geophysical Research, 95:10291–10311, July 1990. doi: 10.1029/JA095iA07p10291. Acknowledgments The work leading to this thesis has greatly benefited from the great working atmosphere at the Institute for Theoretical Physics IV (TP4) of the Ruhr University Bochum, especially in our Heliophysics science group lead by PD Dr. Horst Fichtner, whom I am deeply grateful for endless helpful discussions and support during my PhD. I cannot imagine a better way of supervising! His infectious good spirit transfers to the whole group, whose members over the years (Dr. Jens Kleimann, Dr. Klaus Scherer, Dr. Frederic Effenberger, and Dr. Andreas Kopp) I also owe truly many thanks. It was a pleasure working with you all! For providing a good-working infrastructure at TP4, I want to thank Bernd Neubacher, Verena Kubiak, Gisela Buhr and Prof. Dr. Reinhard Schlickeiser. Speaking of infrastructure, particular thanks are also in order to Prof. Dr. Ralf Kissmann for allowing me to use his code Cronos, which is the main numerical tool used in this thesis, but also for his interest in my work and his help in maintaining and extending the code. I also had the pleasure to co-operate with many interesting people during the last years’ projects, and I enjoyed working with M.Sc. Patrick Kühl, Prof. Dr. Bernd Heber, Dr. Roelf Du Toit Strauss, Dr. Eugene Engelbrecht, Prof. Dr. Marius Potgieter, Prof. Dr. Stefan Ferreira, Dr. Robert Cameron, Dr. Nick Arge, Prof. Dr. Sean Oughton, Dr. Arcadi Usmanov, M.Sc. Laxman Adhikari, Dr. Alexander Dosch and Prof. Dr. Gary Zank. For her volunteering as second reviewer, I would like to express my thanks to Prof. Dr. Julia Tjus. For everything that is truly important in life: I thank my wife Stefanie, my parents Michael and Heike, my brother Lukas, my grandparents, and everyone I may call my friend, for always being there for me! 107 Ich versichere hiermit, dass ich diese Arbeit selbstständig verfasst und keine anderen als die angegebenen Quellen und Hilfsmittel benutzt sowie Zitate kenntlich gemacht habe. Tobias Wiengarten Curriculum vitae Personal data Name Address E-Mail Date of birth Nationality Famility status Tobias Wiengarten Arnold-Berthold-Weg 3, 59494 Soest [email protected] 11 March 1987 in Soest German married Education since 11/2011 Ph.D., Ruhr-Universität Bochum. Thesis: Numerical modeling of plasma structures and turbulence transport in the solar wind 10/2009 – 11/2011 M.Sc., Ruhr-Universität Bochum, Grade: 1.0. 10/2006 – 08/2009 B.Sc., Ruhr-Universität Bochum, Grade: 1.5. 08/1997 – 06/2006 Abitur, Aldegrever Gymnasium Soest, Grade: 1.6. Thesis: MHD Simulation of the Inner-Heliospheric Magnetic Field Thesis: Vergleich verschiedener Magnetfeldkonfigurationen im Bereich des heliosphärischen Terminationsschocks Publications K. Scherer, A. van der Schyff, D. J. Bomans, S. E. S. Ferreira, H. Fichtner, J. Kleimann, R. D. Strauss, K. Weis, T. Wiengarten und T. Wodzinski. Cosmic rays in astrospheres. Astronomy and Astrophysics, 576:A97, Apr 2015. doi:10.1051/0004-6361/201425091. T. Wiengarten, H. Fichtner, J. Kleimann und R. Kissmann. Implementing Turbulence Transport in the CRONOS Framework and Application to the Propagation of CMEs. The Astrophysical Journal, 805:155, Jun 2015. doi:10.1088/0004-637X/805/2/155. K. Scherer, H. Fichtner, J. Kleimann und T. Wiengarten. Large-scale structures in the heliotail and their implication for cosmic ray transport. In 40th COSPAR Scientific Assembly, Vol. 40 von COSPAR Meeting. 2014. T. Wiengarten, J. Kleimann, H. Fichtner und R. Kissmann. MHD Simulations: Corotating Interaction Regions. In N. V. Pogorelov, E. Audit und G. P. Zank, Hrsg., 8th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2013), Vol. 488 von Astronomical Society of the Pacific Conference Series. Sep 2014. T. Wiengarten, J. Kleimann, H. Fichtner, P. Kühl, A. Kopp, B. Heber und R. Kissmann. Cosmic Ray Transport in Heliospheric Magnetic Structures. I. Modeling Background Solar Wind Using the CRONOS Magnetohydrodynamic Code. The Astrophysical Journal, 788:80, Jun 2014. doi:10.1088/0004-637X/788/1/80. T. Wiengarten, J. Kleimann, H. Fichtner, R. Cameron, J. Jiang, R. Kissmann und K. Scherer. MHD simulation of the inner-heliospheric magnetic field. Journal of Geophysical Research (Space Physics), 118:29–44, Jan 2013. doi:10.1029/2012JA018089. K. Scherer, H. Fichtner, F. Effenberger, R. A. Burger und T. Wiengarten. Comparison of different analytic heliospheric magnetic field configurations and their significance for the particle injection at the termination shock. Astronomy and Astrophysics, 521:A1, Okt 2010. doi:10.1051/0004-6361/200913638. Conference and workshop participations 03/2015 DPG Spring-meeting Wuppertal; Talk: Incorporating Turbulence Transport in the CRONOS MHD framework 02/2015 Workshop on Reconnection, Turbulence, and Particles in the Heliosphere in Queenstown, New Zealand; Talk: Coupling Turbulence Transport and WSA-driven MHD simulations 01/2015 Workshop Cosmic Rays in Astrospheres in Bad Honnef 03/2014 Workshop Cosmic Rays: From the Galaxy to the Heliosphere; a numerical modeling approach in Potchefstroom, Southafrica; Talk: Modeling heliospheric background solar wind using the CRONOS MHD code 01/2014 ISSI Team Meeting Heliosheath Processes and Structure of the Heliopause: Modeling Energetic Particles, Cosmic Rays, and Magnetic Fields in Bern, Switzerland 12/2013 AGU Fall Meeting in San Francisco, USA; Poster: MHD simulations: Corotating Interaction Regions 11/2013 Workshop Wind Bubbles, Astrospheres and the Heliosphere: Environments and Cosmic Rays in Bochum; Poster: MHD simulations: Corotating Interaction Regions 09/2013 Space Science Training Week 2013: Data Driven Modeling and Forecasting in Leuven, Belgium; Book award for best formulated research proposal 07/2013 8th International Conference of Numerical Modeling of Space Plasma Flows (ASTRONUM 2013) in Biarritz, France; Talk: MHD simulations: Corotating Interaction Regions 03/2013 527. Wilhelm und Else Heraeus-Seminar: Plasma and Radiation Environment in Astrospheres and Implications for the Habitability of Extrasolar Planets in Bad Honnef; Poster: MHD simulations of the inner heliosphere 02/2013 DPG Spring-meeting Jena; Poster: MHD simulations of CIRs 9th Heidelberg Summer School: Computational Astrophysics 09/2012 03/2012 DPG Spring-meeting Stuttgart; Talk: MHD Simulation of the Inner-Heliospheric Magnetic Field 07/2011 Joint Space Weather Summer School 2011 in Huntsville, USA Soest, June 18, 2015 Tobias Wiengarten