Appendix A: Monte Carlo simulations of particle propagation through self-generated waves

The present numerical code is modified from our previously published Monte Carlo code (Vainio et al. 2000). Similar transport codes, including also non-linear effects, have been developed previously (e.g., Ellison et al. 1996), but to our knowledge, ours is the first Monte Carlo code to employ self-generated waves.

When the Monte Carlo particles (numbered by j) are injected in the simulation, they are a given a weight w_j that normalizes their injection rate to the physical value given by Q. The particles are treated in the guiding center approximation, which in the present case of constant magnetic field strength means that during the Monte Carlo time step, $\delta t$, we move the particles in the spatial coordinate (measured along the field lines) as $x_{j}\leftarrow x_{j}+v_{j}\mu _{j}\, \delta t$, where v_j and $\mu _{j}$ are the particle speed and pitch-angle cosine of the jth particle as measured in the fixed frame, where the background plasma is assumed to be stationary. In addition, the particles suffer scatterings from two Alfvén wave fields propagating parallel and anti-parallel to the magnetic field. In scatterings, performed after each Monte Carlo time step modeling pitch-angle diffusion, the particles are subsequently scattered (elastically) in the two wave frames, first Lorentz-transforming the particle velocity to the wave frame, then using

$$\mu _{j\pm }\leftarrow \mu _{j\pm }\cos \vartheta _{j\pm }+(1... ...^2_{j\pm })^{1/2}\sin \vartheta _{j\pm }\cos \varphi _{j\pm },$$ (A.1)

to compute the new wave-frame pitch-angle cosine $\mu _{j\pm }$, and finally transforming back to the fixed frame. Here $\vartheta _{j\pm }$ and $\varphi _{j\pm }$ are the scattering angles measured from and around the particle's velocity vector $\vec{v}_{\pm }$ prior to the scattering. As described in Vainio et al. (2000), $\vartheta ^{2}_{j\pm }$ and $\varphi _{j\pm }\in [0,2\pi )$ are random numbers picked up from exponential and uniform distributions, respectively, and $\langle \vartheta _{j\pm }^{2}\rangle =2\nu _{j\pm }\, \delta t$, where $\nu _{j\pm }\propto U_{\pm }(x_j,t)$ is the scattering frequency related to the wave-energy density by Eq. (1). Each Monte Carlo time step is chosen small enough to keep $\langle \vartheta _{j\pm }^{2}\rangle < 0.2$ for all j.

We keep track of the wave energy densities on a spatial grid with N=150elements numbered by i=1,...,N, and with a spacing of $\Delta x=L/N$ and central coordinates x=X_i. The value of $U_{\pm }(x)$ is taken to be constant inside each grid cell. Thus, $U_{\pm }(x_{j})=U_{\pm }(X_{i_{j}})$ and i_j is the index of the grid cell that contains x_j, i.e., $x_{j}\in I_{i_{j}}$, where $I_{i}=[X_{i}-\Delta x/2,X_{i}+\Delta x/2)$. During each Monte Carlo time step, a change of the wave energy density due to wave-particle interactions is computed from

$${U_{\pm }(X_{i})\leftarrow}U_{\pm }(X_{i})+\frac{1}{A\, \Delta x}... ...i}}w_{j}\, \nu _{\pm }(X_{j})\delta t\, V_{\mathrm{A}}\, p_{j\pm }\mu _{j\pm },$$ (A.2)

where A is the cross-sectional area of the flux tube. This means that each Monte Carlo particle is taken to represent a large number of particles with the same parallel momentum, and the energy increment of the waves is computed from the mean value of the plasma-frame energy loss in wave-particle interactions of such particles. This method is much more stable computationally than an alternative method, where one computes the plasma-frame energy change in the scatterings of the Monte Carlo particles directly, and give this to the waves. After each Monte Carlo time step, we check that the wave-energy densities are above the minimum level given as a parameter.

In addition to the wave-particle interactions, the waves are convected on the grid by $U_{\pm }(X_{i})\leftarrow U_{\pm }(X_{i\mp 1})$ if $1\leq i,\, i\mp 1\leq N$ and U_+[-](X_1[N])=U₀, each time the simulation time t has elapsed an amount of $\Delta x/V_\mathrm{A}$. At boundaries, the wave energy convected out of the grid is lost.

As an output, the simulation code saves the momentum, pitch-angle cosine and escape time of all particles leaving the simulation box. In addition, the wave-energy densities and the energetic particle pressure inside the simulation box is saved after every time the waves are moved on the grid, i.e., at t mod  $\Delta x/V_{\mathrm{A}}=0$. To illustrate the typical development of a simulation, we have plotted in Fig. A.1 the wave-energy densities and particle pressures in a few frames for simulation in Case B with source position at $x=0.1\, L$ (see Fig. 3 for the escaping particle flux).

Figure A.1: Snapshots from a simulation with an asymmetric position of the source (Case B), $x_{0}=0.1\, L$. The solid ( dashed) curve gives U+(-) and the dashed-dotted curve the energetic-particle pressure $P_{\rm c}$ inside the loop as a function of position at indicated times after the start of the injection. The dotted lines give the positions $x=x_0\pm V_{\mathrm{A}}t$.

As an outline of the future work, we note that a generalization to a spectrum of Alfvén waves would mean that the wave-energy density grid would contain another dimension (wavenumber), and that the resonance condition would have to be taken into account in deciding which particles contribute the the growth of the waves in the particular grid element. If one uses the full quasi-linear resonance condition with $\mu$-dependence, one has to also take this into account when modeling the scattering frequency, which has to be allowed a $\mu$-dependence. This, naturally, affects also the growth rate of the waves. The effects of non-constant magnetic field to the particle transport are easy to take into account (Vainio et al. 2000); for wave transport, one has to use an equation that employs the diverging field effects as well as the effects of a non-constant group speed.