Numerical simulation of oscillating magnetospheres with resistive electrodynamics

We present a model of the magnetosphere around an oscillating neutron star. The electromagnetic fields are numerically solved by modeling electric charge and current induced by the stellar torsional mode, with particular emphasis on outgoing radiation passing through the magnetosphere. The current is modeled using Ohm's law, whereby an increase in conductivity results in an increase in the induced current. As a result, the fields are drastically modified, and energy flux is thereby enhanced. This behavior is however localized in the vicinity of the surface since the induced current disappears outwardly in our model, in which the exterior is assumed to gradually approach a vacuum.


Introduction
Giant flares in Soft Gamma-ray Repeaters (SGRs) are some of the most energetic astrophysical events. The total energy released during a flare is ∼ 10 46 ergs, and flares are considered to be powered by a global reconfiguration of magnetic field > 10 14 G in a magnetar (Mereghetti (2008) for a review). Global seismic oscillations of the stellar crust are likely to be associated with flares (Duncan 1998). Quasi-Periodic Oscillations (QPOs) are observed in the decaying tails of the flares of SGR 1900+14 and SGR 1806-20, and their frequencies are in the range of 10-10 3 Hz, which are interpreted as the frequencies of torsional oscillations with some harmonic numbers (e.g., Israel et al. (2005); Strohmayer & Watts (2006)). The typical frequency of the torsional shear mode is ∼ 10Hz, which is estimated from the shear modulus, density at the crust and stellar radius. More realistic treatment is necessary for the actual fitting. In particular, considering the effects of strong magnetic fields, which couple core and crust, is essential. A number of researchers have been exploring the oscillations of strongly magnetized stars (Piro (2005); Levin (2007); Lee (2008); Sotani et al. (2008); Cerdá-Durán et al. (2009) ;Colaiuda & Kokkotas (2011);van Hoven & Levin (2012); Gabler et al. (2012) and references therein).
In spite of substantial improvements in oscillation theory, there is considerable uncertainty about the processes behind X-and gamma-ray radiation phenomena. Main peaks in the power spectrum may correspond to the oscillation frequency, but their amplitudes and decays should be affected by the surrounding magnetosphere. Timokhin et al. (2000) showed that an oscillating magnetized star should have a magnetosphere filled with charge-separated plasma so that the longitudinal electric field (parallel to the magnetic field) in vacuum is screened 1 . This situation resembles that of a pulsar magnetosphere. Rotation in pulsars is replaced by oscillation in magnetized stars. The induced charge density has been calculated as a generalization of the well-known Goldreich-Julian density for pulsars (Goldreich & Julian 1969). In the almost analytic treatment in Timokhin et al. (2000), the current is assumed to be zero everywhere. In magnetar flares, large current associated with a magnetic twist ∇ × B = 0 is expected to flow out from the surface. Therefore, it is necessary to consider the effects of flowing current on the magnetosphere. The calculations necessary to consistently determine the charge and current densities induced by the oscillation are more complex, requiring time-dependent numerical simulations. One approach is to perform relativistic magnetohydrodynamics (MHD) simulations. However, this approach is rather difficult in the case of strong magnetic regimes, such as that in magnetars, since plasma effects are less important. Recently, a new numerical approach was proposed for the construction of pulsar magnetospheres (Gruzinov 2011;Li et al. 2012;Kalapotharakos et al. 2012). In this approach, the Maxwell equations are solved for the current modeled using Ohm's law, and particle acceleration as well as the generation of radiation are relevant to the electric field component parallel to the magnetic field, that is, B · E = 0. Thus, such dissipative effects of non-ideal MHD are taken into account in the model, which differs from the force-free electrodynamics approach, in which plasma can be approximated as massless and ideal MHD conditions are assumed to hold everywhere.
We adopt the same resistive electrodynamics approach in order to simulate global structures in the flares. The electromagnetic fields outside a star are disturbed by periodic torsional motion, which is horizontal at the surface. There is no justification for the use of resistive electrodynamics since the portion of plasma energy with respect to electromagnetic energy ejected from the surface is unknown. Moreover, the physical conditions relevant to QPO phenomena of the flares is uncertain. Non-thermal X-and gamma-rays in pulsars originate from the effects of non-ideal MHD, namely accelerating electric fields. It is therefore a reasonable first step to simulate the disturbances in the magnetosphere by analogy with pulsars, in whose framework non-ideal MHD is allowed. The results are expected to provide useful insight into observed phenomena and to lead to further improvement of theoretical models. This paper is organized as follows. In Section 2, we discuss the proposed numerical method for the evolution of electromagnetic fields in two-dimensional axially symmetric space. We adopt spherical coordinates (r, θ, φ) and assume axial symmetry (∂/∂φ = 0). However, in general, the azimuthal component of a vector is nonzero. For example, a toroidal magnetic field B φ can be produced by a poloidal current (j r , j θ ). Disturbances associated with shear oscillations are assumed at the inner boundary. In Section 3, we present the results of a numerical simulation based on the proposed model, and Section 4 presents a discussion of the results.

Electromagnetic fields
Maxwell's equations are solved by assuming a charge density ρ e and a current density j: To solve these equations, we use a scalar potential Φ and a vector potential A satisfying the Coulomb gauge, ∇ · A = 0. For axially symmetric fields, the following form given by two functions, F (t, r, θ) and G(t, r, θ), is automatically satisfied with the gauge condition: where e φ is a unit vector in the azimuthal direction. One of the advantages in the potential formalism is that the constraint Eq.(3) is automatically satisfied. Helmholtz's vector decomposition due to the divergence and rotation of j is also clear for these functions as discussed below. The electric field E and magnetic field B are given by the derivatives of the potentials: The function S in B φ is given by where the differential operator D is defined as The electric field (Eq. (6)) given by the time derivative of F is used to facilitate solving the time derivative of Eq. (8), i.e., Substituting these forms (Eqs. (6) and (7)) into Maxwell's equations, we obtain two wave equations for G and S, and a Poisson equation for Φ: The evolution of these functions is based on a model of current density. We briefly introduce the model presented in Li et al. (2012). Ohm's law relating the current and the electric field is where σ denotes electric conductivity and Eq. (14) is a representation in the fluid rest frame. Upon performing transformation into quantities in some inertial frame, the general form of the current model is given by where and E 0 , B 0 are given by the two Lorentz invariants , and E 0 ≥ 0. The drift velocity contains the term E 2 0 in the denominator to account for the non-zero electric field in the fluid rest frame. This correction ensures | v D | ≤ c, even if the magnetic field vanishes. The term Γ in Eq. (15) represents the effects of non-ideal MHD. Various model prescriptions are possible since there is ambiguity concerning the plasma velocity, in which frame Ohm's law holds. For a choice of 'minimal velocity', 2 Γ is given by (Li et al. 2012), For comparison, the expression Γ G proposed by Gruzinov (2011) is where the charge density is assumed to vanish in the fluid rest frame. In addition, note that Γ is not necessarily positive, and it is possible to model the energy transfer to electromagnetic fields in the region where Γ < 0. Thus, at present there is no unique form for non-ideal dissipation. We adopt Eq. (17), which is the simplest form, and easily compare our results with those in the pulsar magnetosphere by (Li et al. 2012). In more general cases, the plasma velocity should vary temporally and spatially, and would provide some complicated structure. Obtaining a realistic form of σ strongly depends on the spatial position. The conductivity is likely to decrease to zero with the radius r since the plasma number density decreases. To derive the behavior toward the vacuum region, we use the simple power law form where σ 0 is a constant and r 0 is the inner boundary radius.
Finally, the time evolution of ρ e is given by the charge density conservation law The numerical procedure for calculating the electromagnetic fields is summarized here. Once the charge density ρ e and current density j are given at a particular time, the time evolutions of G, S and ρ e are found by solving Eqs. (11), (12) and (20) with appropriate boundary conditions. The evolution of these functions is separately governed by the nature of vector field j. More specifically, G is governed by the toroidal component j φ , S is governed by the rotation of poloidal components ∇ × j p , and ρ e is governed by their divergence ∇ · j p . The elliptic equation (13) is solved at each time step for the function Φ with ρ e . There is another elliptic equation (10) for the function ∂F/∂t. Thus, the electromagnetic fields are given by Eqs. (6) and (7) in terms of the functions G, S, Φ and ∂F/∂t, and the current j is also updated by Eq. (15).

Torsional oscillation
We consider disturbances in the electromagnetic fields in the presence of torsional oscillation within a star. Before the onset, the exterior fields are in the form of a purely magnetic dipole given by where µ is the dipole moment and the typical field strength B s at the surface r 0 is B s = µ/r 3 0 . The magnetic flux function is We examine the axially symmetric torsional oscillation, whose motion is horizontal in the azimuthal direction at the surface. In general, the velocity associated with the sinusoidal oscillation with angular frequency ω is given by (Unno et al. 1989) where P l (θ) is a Legendre function of order l, and ǫ denotes the dimensionless amplitude of the oscillation. The frozen-in condition, E + v × B/c = 0, results in a poloidal electric field E p at the surface: In the exterior vacuum region, the torsional mode with l in general induces outgoing radiation with l ± 1 (McDermott et al. 1984). Below, we consider oscillation with l = 2 only in order to examine dipole radiation, in which case the damping is most effective. The azimuthal component E φ is not induced since it belongs to a different type of mode. Electromagnetic wave solutions in vacuum are classified into two modes: the TE mode ( E p , B φ ) and the TM mode ( B p , E φ ) (Jackson 1975). The latter is described solely by the function G, which is given by the continuity of B r . That is, G at r 0 is always fixed as a dipolar function (Eq. (22)). The azimuthal component of the electric field (E φ = 0) also satisfies the continuity across the surface. On the other hand, B θ is not specified at the surface, since it is described by the derivative ∂G/∂r, and it is not in general continuous due to the surface currents.
We discuss the boundary conditions of TE part ( E p , B φ ). The tangential component E θ should be continuous. Using E The radial component E r is not continuous to E r ( Eq.(24) ) due to the surface charges. Remaining boundary condition for B φ (= S/(r sin θ)) is specified by a physical argument. The derivative ∂B φ /∂r is related with ∂E (0) θ /∂t and j θ in Eq.(2). By assuming that the tangential component j θ vanishes, the explicit form of S at r 0 can be described by The assumption j θ = 0 is good at least in the polar regions, since the currents flow along the magnetic field, which is radial near the poles. The boundary condition (26) is not a unique one, since there exist a number of models for the out-flowing currents at r 0 . In exterior vacuum, i.e., j = 0 everywhere, the condition (Eq. (26)) provides some amount of outgoing energy power, which is analytically calculated. The model with Eq.(26) is therefore useful for comparison. We consider the response of the magnetosphere to the disturbances, which are caused with the same boundary value as that for the vacuum. In this paper, we discuss how and where the ideal MHD condition is broken and the resistive effect works in the magnetosphere.
Here, we estimate the outgoing energy flux for the boundary condition (Eq. (26)). The timeaveraged Poynting flux at the wave zone contains dipole (l = 1) and octupole (l = 3) radiation. By matching B φ at r 0 with the wave solution of l = 1 (Jackson 1975;McDermott et al. 1984), we obtain the luminosity as where x 0 = ωr 0 /c. At the low frequency limit, this is reduced to L ∝ (ωr 0 /c) 4 (ǫr 0 B s ) 2 c. The dependence on ω is due to dipole radiation. Similarly, the luminosity of octupole radiation is of the order of (ωr 0 /c) 8 (ǫr 0 B s ) 2 c, and is negligibly small for ωr 0 /c ≤ 1.

Numerical method and model parameters
We adopt the finite difference method in order to solve the partial differential equations. The numerical domain of the spherical coordinates is 1 ≤ r/r 0 ≤ 150 and 0 ≤ θ ≤ π. The typical number of cells on the grid is 300 and 180 in the directions of r and θ, respectively. The grid cell spacing in the radial direction is taken as ∆r ∝ r 1/2 to obtain fine resolution near the inner region, whereas the spacing in the angular direction is constant. The boundary conditions at the inner boundary are already discussed in the previous subsection. The other conditions are regularity conditions at the polar axis (θ = 0, π) and outgoing wave conditions at the outer boundary. The initial conditions assume a pure dipole, that is, the function G is given by Eq. (22), and the other functions are zero: S = Φ = ∂F/∂t = ρ e = j = 0.
In the numerical simulations, the oscillation frequency and amplitude are chosen to be larger than realistic values. The QPO frequencies observed in magnetar flares are of the order of ν = 10 1 -10 3 Hz. In particular, 30 Hz QPOs can be interpreted as nodeless shear oscillations with l = 2 (Strohmayer & Watts 2006), for which the frequency is written as ν ∼ 10 −3 c/R s using the stellar radius R s and c. In our numerical model, the frequency is chosen as ν = ω/(2π) = 0.1c/r 0 (wavelength λ = 10r 0 ) in order to reduce the model size and facilitate the evaluation of the luminosity, which is typically given by Eq. (27). The frequency is scaled up by a factor of 10 2 if the inner radius r 0 is taken as the stellar surface radius R s . The scaling of ν does not hold exactly since the conductivity depends on the spatial position. The dynamics of our simple model is governed by a dimensionless parameter (the ratio of electric conductivity to the frequency), as described below. In future research, it would be necessary to study the dynamics by devising a more realistic conductivity model. The oscillation amplitude is also slightly higher, where ǫ = 0.1 is used in the numerical calculations, whereas the actual amplitude of flares is ǫ ≪ 1 in Eq. (23). The scaling of ǫ is confirmed up to ǫ ≤ O(1).
The electric conductivity σ is key to determining the electromagnetic field structure. For example, in the case of σ = 0 everywhere, the electric current always vanishes when the initial conditions assume a purely magnetic dipole, and hence the state is stationary. The magnitude σ 0 in the form of the power law in Eq. (19) is normalized as p ≡ σ 0 (4π/c)(c/ω) = 2σ 0 /ν, where 4π/c ≈ 377ohms is the characteristic impedance of vacuum 3 . The dimensionless parameter p = 2σ 0 /ν denotes the ratio of the current to the displacement current: p = 4πj/(ωE) = 4πσ 0 /ω. If p ≫ 1, the displacement current is negligible, whereas the limit of p = 0 corresponds to vacuum. In magnetar flares, a sufficient amount of electron-positron pair plasma may be produced, which fills the magnetosphere (Beloborodov & Thompson 2007). Therefore, the case of extremely small p is not applicable, although it is still useful for comparison. Figure 1 shows the luminosity, which is the radial component of the Poynting flux through a surface. From Eq. (27), for our numerical values (νr 0 /c = ǫ = 0.1), the typical value of the luminosity is L/((r 0 B s ) 2 c) ∼ 10 −3 . The luminosities for σ 0 /ν = 5(top panel) and 0.5 (bottom panel) are shown. They are estimated at two different radii (r = λ(= 10r 0 ) and 4λ). Outgoing waves are modified in the presence of resistivity. Similar sinusoidal curves with a period of 5r 0 /c can be seen in the figure 4 . The amplitudes of the waves depend on both the estimated radius and the parameter σ 0 . Regarding the flux at r = λ, the oscillating amplitude for higher conductivity is slightly larger than that for lower conductivity. The amplitudes decrease with the radius. Large current is temporarily induced within the inner region and tends to zero in our model with σ ∝ r −1 . The current significantly modifies the electromagnetic fields and allows for a greater Poynting flux. With the increase of the parameter σ 0 , the induced current also increases. The curves in Fig. 1 are almost completely described by a single Fourier component, and are therefore different from the X-ray light curve corresponding to flares. The complexity in the observed power spectrum, such as broad peak width and rapid time variation in the QPO, is likely to reflect the nonlinearity of the dynamical system and/or the emission mechanism. In the numerical simulation, a more complex profile can be produced, for example, by further increasing the current. The parameter range σ 0 /ν > 10 is more interesting. Actually, we tried and found that the fluctuations in the small scale are easily excited, and sometimes grow. This local instability may originate from our numerical scheme or the current model. Performing a comparison also becomes difficult since care should be taken to ensure the validity and efficiency of analysis in the case of highly complex and time-dependent data. This regime is beyond the scope of this paper, and the parameters are limited to certain values σ 0 /ν = p/2 ≤ O(10).

Results
We examine the dependence of σ 0 on the luminosity in Fig. 2. Time-averaged values integrated at the surface corresponding to r = λ, 4λ and 8λ are denoted by L 1 , L 4 and L 8 (L 1 > L 4 > L 8 ). The luminosities irrespective to the position (r > λ) approach L * ≡ 9.8 × 10 −4 (r 0 B s ) 2 c in the limit of small σ 0 , where L * is the value of dipole radiation in vacuum, estimated in Eq. (27). All the values L k (k = 1, 4, 8) slightly decrease with increasing σ 0 , but tend to increase for σ 0 /ν > 2. The luminosity minimum is located near σ 0 /ν ∼ 2. The dependence of σ 0 can be explained by two effects. One is the Ohmic dissipation:the decay timescale increases with it. The other is flowing currents produced in our model:the magnitude generally increases with σ 0 . In the vacuum limit, σ ≡ 0, there are no currents and no dissipation of the electromagnetic energy. In the small region, the energy is quickly dissipated, and the reduction of the flux is more evident at great distances. With increase of σ 0 , the Ohmic dissipation becomes less effective, and the amount of flowing currents increases. They produce electro-magnetic field structure different from that of vacuum in large σ 0 region. The induced fields dominate, and the outgoing flux is independent of the inner boundary condition as long as it is estimated outside the induction region.  Figure 3 shows time evolution of electromagnetic fields. It is almost periodic, and the evolution during half of the oscillation period (0.5ν −1 = 5r 0 /c) is shown. The top and bottom panels show the results for low conductivity (σ 0 /ν = 0.5) and high conductivity (σ 0 /ν = 5), respectively. The contours of the function S(= B φ r sin θ) clearly represents a wave propagating at speed c. The oscillation behavior becomes apparent in the wave zone r > λ = 10r 0 . The angular dependence at which the amplitude is maximum on the equator forms a dipolar radiation pattern. Although octupole radiation may be involved, its effects are rather small. The magnetic flux function G describing poloidal magnetic fields is also plotted. The pattern is dipolar in the inner region r < 10r 0 but gradually deviates from this as the radius increases. This fact becomes much clearer at high conductivity. The contour of G in Fig. 3 shows remarkable topological changes, i.e., reconnection due to resistivity: a small oval becomes detached from the outer part (r ∼ 15r 0 , θ ∼ π/2) and shrinks beyond that region. In addition to the interesting behavior of poloidal fields, the magnitude of B φ is considerably increased. Thus, poloidal and toroidal components of the magnetic fields are effectively coupled with each other through induced current at high conductivity. The induced current disappears at larger radii in our model (σ → 0), and the magnetic field eventually becomes dipolar with some radiative component at infinity.

Discussion
We studied time variation in a magnetosphere disturbed by torsional shear oscillation by using resistive electrodynamics. In this formalism, the current is described by Ohm's law and the conductivity, and the electromagnetic fields are solved numerically. The results indicate that the induced current increases with the increase in conductivity, the fields are substantially modified, and the energy flux is thereby increased. Beyond a certain critical value of the conductivity σ 0 /ν ∼ O(1), induced components dominate, and the inner boundary condition becomes less important. The behavior is however localized since the induced current disappears in the outer part. The exterior is assumed to approach vacuum in our model, in which the conductivity approaches zero. The extent of energy loss in the form of radiation is almost the same as that in vacuum. It is interesting to compare our results with those for a pulsar magnetosphere by Li et al. (2012), where the conductivity is assumed to be uniform everywhere. The results in Li et al. (2012) show that the luminosity increases with conductivity and approaches that of force-free calculation in the limit of high conductivity (σ 0 /ν ≫ 1). Mathematically, the force-free limit is written as σ → ∞ and E 0 → 0, but σE 0 is finite in Eq. (17). In our numerical simulation, we found that E 0 does not necessarily approach zero everywhere (E 0 B 0 = E · B = 0). The electromagnetic wave components E w and B w satisfy E w · B w = 0, however, there is another dipole component B d , for which E w · B d = 0. It is not clear whether the problem arising in the limit of high conductivity depends on the spatial distribution model or on the numerical method.
The electric conductivity differs depending on the direction in the strongly magnetized regime, -Time evolution of toroidal and poloidal magnetic fields in meridional plane (0 ≤ θ ≤ π and r ≤ 25r 0 = 2.5λ). Color contour represents S(= B φ r sin θ) normalized by B s r 0 . Lines represent the level of the magnetic function G/(B 0 r 2 0 ) for 0.03×n(n = 1, 2, · · · , 7). The innermost contour line is written in black. Contours of larger G within it, which are almost unchanged, are omitted. The top and bottom panels show the results for low conductivity (σ 0 /ν = 0.5), and high conductivity (σ 0 /ν = 5), respectively. From left to right, the panels correspond to ct/r 0 = 50.0, 51.7, 53.3 and 55.0.
in which case the isotropic form of Ohm's law (14) may be replaced by a tensor. As a brief discussion of directional differences, the conductivity is assumed to arise from collisions, and some terms in the generalized Ohm's law are ignored. The conductivity parallel to the magnetic field is given by σ = e 2 n e /(m e ν c ), and that in perpendicular direction is given by σ ⊥ = σ /(1 + (ω B /ν c ) 2 ), where ω B = eB/(m e c) is the magnetic gyro-frequency and ν c is the collision frequency (e.g., Somov (2006)). Using the Goldreich-Julian number density n e = Bν/(ec), the ratio is reduced to σ ⊥ /σ = (1 + (σ /ν) 2 ) −1 . The isotropic form may be no longer valid in the limit of σ /ν ≫ 1, where σ ≫ σ ⊥ . In our numerical simulation, the parameter is limited to σ 0 /ν < 10. The anisotropic form of Ohm's law would be adequate for larger values of to σ 0 /ν. Thus, conductivity is an important factor, but its magnitude and spatial distribution are poorly understood at present. A more elaborate treatment is necessary for constructing a more realistic model. The outgoing energy flux is essential for estimating the damping mechanism of the oscillations. The luminosity L calculated in this paper is of the same order as that in vacuum, although it depends on the conductivity. Appropriate scaling of the shear oscillation with l = 2 leads to L ∼ 2 × 10 37 (ν/30Hz) 4 (ǫ/10 −2 ) 2 (B/10 14 G) 2 ergs s −1 . The damping timescale for the energy E ∼ 10 45 ergs is t d = E/L ∼ a few years. This is a rather long period. It is difficult to significantly increase radiation loss, even if there exists a magnetosphere. The loss to infinity is not so efficient in the exterior vacuum region. The oscillation energy slowly escapes into the radiative zone, which is located beyond the region corresponding to a few wavelengths. In the presence of a magnetosphere, both ingoing and outgoing fluxes are likely to be produced due to interaction inside that zone, and current flows also change direction there. The inward energy flux and return current flows may be non-negligible in magnitude and may be dissipated as Joule heating in the region with higher resistance. Thus, the problem is extended to a globally coupled system between a neutron star/crust and a magnetosphere, which is more complex to model. The dynamics of coupled oscillations is an interesting topic for future research.
A model of a twisted magnetosphere would be more realistic since global electric current flows through the stellar surface (Thompson et al. 2002;Beloborodov & Thompson 2007). The exterior of a star before flares is modeled by a purely dipole field, and in the present model current flows are allowed only during the oscillations. However, the effects of twisting may be significant. An interesting simulation by using a time-dependent force-free electrodynamics model was recently presented in Parfrey et al. (2012), where the evolution of the twisted magnetosphere was calculated. It was found that slow shearing at the surface leads to twisting and an abrupt reconnection. The evolutionary time scale (≫ 1 s) of that model is longer than that of ours (< 1 s). Such a gradual twisting process would also be an interesting topic for future research.