**The inner magnetosphere particle transport and acceleration model (IMPTAM), developed by Ganushkina et al. (2001, 2005, 2006),
follows distributions of ions and electrons with arbitrary pitch angles from the plasma sheet to the inner L-shell regions
with energies reaching up to hundreds of keVs in time-dependent magnetic and electric fields. **

We trace a distribution of particles in the guiding center, or drift, approximation,
in which we can picture the motion of a charged particle as displacements of its guiding center,
or the center of the circular Larmor orbit of a moving particle.
The guiding center theory assumes that the electromagnetic fields are known and can be used in geophysical plasmas,
where the external field is strong and will not be changed much by the motion of the particle themselves.

As guiding center drifts we take into account **E × B** drift,
where** E** and **B** are electric and magnetic fields, respectively,
and magnetic drift, which, in its turn, includes gradient and curvature drifts.
The drift velocity is a combination of the velocity V_{E×B} due to **E × B** drift
V_{E×B} = (**E × B**)/B^{2}
and the velocities of gradient V_{∇} and curvature V_{cur} drifts

V_{∇} + V_{cur} = (mv_{⊥}^{2})/(2qB^{2})(**B**× ∇B)+(mv_{∥}^{2})/(qR_{c}^{2}B^{2})(**R _{c}**×B)
(Roederer, 1970),

where m is the particle mass, q is the particle charge, v

We assume that the first and second adiabatic invariants are conserved.
The first adiabatic invariant for nonrelativistic particles is the particle magnetic moment given by

μ = p_{⊥}^{2}/(2mB),
where p_{⊥} = p is the particle's momentum in the guiding center system (Roederer, 1970).
The magnetic moment of a particle is conserved in all cases, even in non-stationary fields,
as long as the guiding center approximation is valid.
There are two main conditions for guiding center approximation:

1. Spatial variations in the area of Larmor radius ρc are very small, so that ρc(∇_{⊥}B/B) << 1.
Magnetic field B does not vary much along the Larmor orbit.

2. Temporal variations are small in comparison with Larmor period τ_{c}(∇_{⊥}B/B) << 1.

Second adiabatic invariant is related to the bounce motion of a particle along magnetic field line.
The integral J =∲ p_{∥}ds taken along a given fixed magnetic field line for a complete bounce cycle,
where p_{∥} is the particle's momentum parallel to the magnetic field, and ds is the length element of a field line,
is an adiabatic invariant and is conserved during the drift of a trapped particle as long as field variations during
the time of a bounce period τb are small (τ_{b}dB/dt)/B << 1.

The particle's momentum p is constant during one bounce, so J =2pI, where
and S'_{m} and S_{m} are
the mirror points, B(s) is the magnetic field along magnetic field line, B_{m} is the magnetic field at the mirror point.

With the above mentioned assumptions, we consider bounce-average drift velocity after averaging over
one bounce of **E × B** magnetic drift velocities (Roederer, 1970)

where **E**_{0} and **B**_{0} are the electric and magnetic fields in the equatorial plane, respectively,
e_{0} is the unit vector in the direction of the magnetic field **B**_{0}.
In order to follow the evolution of the particle distribution function f and particle fluxes in the inner magnetosphere dependent on
the position R, time t, energy E_{kin}, and pitch angle φ,
it is necessary to specify:

1. particle distribution at initial time at the model boundary

2. magnetic and electric fields everywhere dependent on time

3. drift velocities

4. all sources and losses of particles

Generally, the changes in the distribution function f(R, φ, t, E_{kin}, α),
where R and φ are the radial and azimuthal coordinates in the equatorial plane, respectively, t is the time, E_{kin} is the particle energy,
α is the particle pitch angle, are obtained by solving the following equation:

where V_{φ} and V_{R} are the azimuthal and radial components of the bounce-average drift velocity.

At the beginning of modeling with IMPTAM, the inner magnetosphere is considered empty. In this case,
only the effects of newly entering particles from the plasma sheet are investigated.
The model boundary is set in the plasma sheet at distances, depending on the scientific questions we are trying to answer, from 6.6 to 10 R_{E}.
The particle distribution at the boundary is defined as a Maxwellian or kappa distribution function with parameters obtained from the empirical relations or
from the observations during specific events.

Liouville's theorem is used to gain information of the entire distribution function. If we know the distribution function f(R, φ, t, E_{kin}, α)
of particles at a time moment t_{1},
then we can obtain the distribution function of particles at a time moment t_{2} = t_{1} + λt, by computing the drift velocity of the particles.
The distribution function at t_{2} will not be the same as at t_{1} at the corresponding positions,
since we need to take into account the phase-space-dependent losses (τ_{loss}).
The final distribution function at t_{2} will be f(t_{2}) = f(t_{1}) exp(-λt/τ_{loss}).

Particle loss processes, which are important for modeling the ring current ions, include charge exchange with neutral hydrogen in the upper atmosphere, Coulomb collisions, and convective outflow through the magnetopause. The charge-exchange cross-section is obtained from Janev and Smith [1993]. The thermosphere model MSISE 90 [Hedin, 1991] and the plasmasphere model by Carpenter and Anderson [1992] are used. For the electrons, the loss lifetimes will be calculated via pitch angle scattering by plasma waves, as described below.

An advantage of this model is that it can utilize any magnetic or electric field model, including self-consistent magnetic field.
At present, several versions of the Tsyganenko magnetic field model, including the TS05 model [Tsyganenko and Sitnov, 2005], can be used.
For the electric field, the Volland-Stern [Volland, 1975; Stern, 1975], Boyle [Boyle et al., 1997], and Weimer [Weimer, 1996, 2001, 2005] models are incorporated.
In addition to the large-scale fields, transient fields associated with the dipolarization process in the magnetotail during substorm onset were modeled as
an earthward propagating electromagnetic pulse of localized radial and longitudinal extent [Li et al., 1998; Sarris et al., 2002].
The magnetic field disturbance from this dipolarization process was obtained from Faraday's law.
Several pulses were launched at substorm onset times [Ganushkina et al., 2001, 2005, 2006].
**One of the important results obtained from IMPTAM modeling is the ability of the model
to reproduce the observed amount of ring current protons with energies > 80 keV during a storm recovery phase [Ganushkina et al., 2006],
which was not possible to obtain by other models using a dipole model for the magnetic field and large-scale convection electric field.**

The evolution of modeled proton distributions was followed using combinations of several representations for magnetic and electric fields with different boundary conditions. We use a dipole for the internal magnetic field. For the external magnetic field four different representations were used:

1. no external field sources (dipole only);

2. Kp-dependent Tsyganenko T89c (T89 abbreviation is used throughout the paper) (Tsyganenko, 1989; Peredo et al., 1993);

3. T96 (Tsyganenko, 1995) with Dst, Psw, IMF By and Bz as input parameters;

4. Tsyganenko and Sitnov TS04 (Tsyganenko and Sitnov, 2005) with Dst, Psw, IMF By and Bz and six variables W_{i} ,i = 1,6 as input parameters.
The variables W enter in the six magnitude coefficients for the magnetic fields from each source and calculated as time integrals dependent on solar wind and IMF parameters
from the moment in time when IMF Bz turns southward.

The electric field representations include:

1. Kp-dependent Volland-Stern VS (Volland, 1973; Stern, 1975) convection electric field;

2. Boyle et al. (1997) polar cap potential dependent on so- lar wind and IMF parameters applied to Volland-Stern type convection electric field field.

Four types of boundary proton distributions used in the model are the following:

BC1: Maxwellian distribution function at 6.6 R_{E} with n=0.5cm^{–3} and T =5keV;

BC2: Kappa-type distribution function (κ = 5) with the observed parameters T (< T >= 7.8 keV) and n (< n >= 0.97 cm^{–3}) by LANL MPA (Bame et al., 1993)
in the 3–45 keV energy range at 6.6 R_{E} on the nightside part of the boundary (18:00–06:00 MLT) in the equatorial plane;

BC3: Maxwellian distribution function at 10 R_{E} with parameter n given by the empirical relation between the plasma sheet number density
and the solar wind number density nps = 0.025 nsw + 0.395 (Ebihara and Ejiri, 2000) and T = 5 keV;

BC4: Maxwellian distribution function at 10 R_{E} with n and T given by the empirical model derived from Geotail data by Tsyganenko and Mukai (2003).