Adiabatic losses

This notebook demonstrates a basic GT simulation workflow for numerically integrating the equation of motion of a charged particle in the Heliosphere with the inclusion of adiabatic losses of energy.

Goals of this example:

• Define the Parker magnetic field and initial conditions for protons with kinetic energy in range from 300 MeV to 30 GeV, from starting from \((1, 0, 0)\) AU.

• Compute the trajectory with Buneman–Boris (BB) and an energy loss due to adiabatic losses.

Expected output:

• A plot of energy loss dependence over time for protons of considered initial energies.

from datetime import datetime

import matplotlib.pyplot as plt
import numpy as np
from gtsimulation.Algos import BunemanBorisSimulator
from gtsimulation.Global import Regions, Units
from gtsimulation.MagneticFields.Heliosphere import Parker
from gtsimulation.Particle import Generators, Flux

Problem setup: field, particle, and integration parameters

This section defines the simulation setup:

• date (used by time-dependent field models);

• the Heliosphere region and the Parker field model;

• in region parameters if we activate additional energy losses;

• particle initial conditions: 10 batches of 10 protons in energy range from 300 MeV to 30 GeV with an initial coordinate \((1, 0, 0)\) AU;

• integration settings: number of steps n_steps, and time step dt;

• which quantities to store in the track (coordinates, velocities, energy, and time along the trajectory);

• simulation using BB scheme.

date = datetime(2008, 1, 1)
region = Regions.Heliosphere
# including the adiabatic energy losses into the calculation
region.value.set_params(CalcAdditionalEnergy=True)
bfield = Parker()

medium = None
use_decay = False
nuclear_interaction = None

energy_array = np.geomspace(3e2, 3e4, 10) * Units.MeV
flux = Flux(
    Spectrum=Generators.Spectrums.UserInput(energy=energy_array),
    Distribution=Generators.Distributions.SphereSurf(
        Center=np.array([1, 0, 0]) * Units.AU,
        Radius=0
    ),
    Names="proton",
    Nevents=energy_array.size
)

save = [100, {'Energy': True, 'Clock': True}]

n_steps = 100_000
dt = 1

# n_steps = 1_000_000
# dt = 0.1

forward_tracking = 1
n_files = 1
output = None

break_conditions = {"Rmax": 80 * Units.AU}
verbose = 0

simulator_BB = BunemanBorisSimulator(
    Bfield=bfield,
    Region=region,
    Medium=medium,
    Particles=flux,
    InteractNUC=nuclear_interaction,
    UseDecay=use_decay,
    Date=date,
    Step=dt,
    Num=n_steps,
    Nfiles=n_files,
    BreakCondition=break_conditions,
    Save=save,
    Output=output,
    ForwardTrck=forward_tracking,
    Verbose=verbose
)

To reduce runtime, the computation is parallelized over batches.

from concurrent.futures import ProcessPoolExecutor

def worker(seed):
    np.random.seed(seed)
    track_list = simulator_BB()[0]
    return track_list

with ProcessPoolExecutor() as executor:
    futures = [executor.submit(worker, i) for i in range(10)]
    results = [f.result() for f in futures]

Energy loss vs time

Here we plot the relative energy loss \(\left(-\dfrac{\Delta T}{T_0}\right)\) dependence on time for protons with considered initial energies.

We have an ensemble of initial conditions, as we considered 10 batches of particles. This allows to estimate the median relative energy loss (the solid line) and the \(15\% \div 85\%\), i.e. \(1\sigma\) spread (the shaded area).

adiabatic_losses = {energy: [] for energy in energy_array}
time_energy = {energy: [] for energy in energy_array}

for event_list in results:
    for event in event_list:
        T0 = event["Particle"]["T0"]
        E = event["Track"]["Energy"]
        t = event["Track"]["Clock"] / Units.day
        adiabatic_losses[T0].append(E)
        time_energy[T0].append(t)

time_min = 2e-3 # days

fig = plt.figure()
ax = fig.subplots()

for T in energy_array:
    i_max = min(arr.size for arr in time_energy[T])
    time = time_energy[T][0][:i_max]
    energy = np.array([arr[:i_max] for arr in adiabatic_losses[T]])
    # limit arrays from below by time
    energy = energy[:, time > time_min]
    time = time[time > time_min]
    # calculate median vaalues and quantiles
    med = np.median(energy, axis=0)
    low = np.quantile(energy, q=0.15, axis=0)
    high = np.quantile(energy, q=0.85, axis=0)
    # plot values
    ax.plot(time, (T - med) / T, label=f"T = {int(T)} MeV")
    ax.fill_between(time, y1=(T - low) / T, y2=(T - high) / T,  alpha=0.5)
ax.set_xscale('log')
ax.set_xlabel('Time [days]')
ax.set_ylabel(r"$-\frac{\Delta T}{T_0}$")
ax.legend()

plt.show()