Notebooks · posts/the-boris-push-gyrates-and-drifts

A gyro-orbit and the E×B drift by the Boris push

Adam Zettel · September 25, 2026 · written with an agent

The Boris push moves a charged particle through electric and magnetic fields in almost every particle-in-cell code. In a magnetic field alone it keeps the speed to rounding and puts every position on one circle; the gyroperiod comes out 6.28449 instead of $2 π$ , longer by $(ω_{c} Δ t)^{2} /12$ . With an electric field added, a particle started at $E \times B / B^{2}$ moves at exactly that velocity forever, and any other start averages to it.

θ = 2 arctan (\frac{ω _{c} Δ t}{2}), v ˉ = \frac{E \times B}{B ^{2}}

A charged particle in a uniform magnetic field circles the field lines at the cyclotron frequency $ω_{c} = q B / m$ . This circling is gyration, and its circle is the gyro-orbit. Add a uniform electric field at right angles to $B$ and the center of the circle, the guiding center, moves at the $E \times B$ drift velocity $E \times B / B^{2}$ , whatever the particle's charge or mass.

A particle-in-cell (PIC) code moves millions of particles through these fields for millions of steps, so the step it uses must not add or remove energy on its own. The Boris push (J. P. Boris, 1970) is the standard choice. It splits the step into half an electric kick, a rotation about $B$ , and the other half kick. This notebook runs it with $q = m = B = 1$ , so $ω_{c} = 1$ and the exact gyroperiod is $2 π$ , and checks each property against a closed form. posts/refs/boris_push.py computes the same numbers with a separate numpy implementation.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 28–47

DT = 0.05
# 20000 steps of 0.05 is about 159 gyroperiods
NSTEP = 20000
PI = 3.141592653589793

#| One Boris step for a velocity v = (vx, vy) in B = ẑ and
#| E = (ex, ey), with q = m = 1. Half kick, rotate by tz = dt/2
#| and sz = 2·tz/(1 + tz²), half kick. Updates v in place.
def boris_step(v: Arr, ex, ey, dt) {
    tz = 0.5 * dt
    sz = 2.0 * tz / (1.0 + tz * tz)
    mx = v[0] + 0.5 * dt * ex
    my = v[1] + 0.5 * dt * ey
    wx = mx + my * tz
    wy = my - mx * tz
    px = mx + wy * sz
    py = my - wx * sz
    v[0] = px + 0.5 * dt * ex
    v[1] = py + 0.5 * dt * ey
}

Gyration keeps the speed and the circle

With no electric field the two kicks vanish and the step is a pure rotation of $v$ , so $∣ v ∣$ cannot change except by rounding. The cell runs 20000 steps from $v = (0, 1)$ and records the largest and smallest speed.

The positions are updated with the new velocity, $x_{n + 1} = x_{n} + Δ t v_{n + 1}$ . Each step turns $v$ by the same angle $θ$ , so each position step is the previous one turned by $θ$ : the positions are the corners of a regular polygon and all lie on one circle. The cell takes the circle through the first three positions and measures every other position's distance from its center. Its radius is $Δ t / (2 sin (θ /2))$ , which is $1 + (ω_{c} Δ t)^{2} /4 v / ω_{c} = 1.000312$ , not the exact $v / ω_{c} = 1$ .

[2]lines 66–115

vel = [0.0, 1.0]
xs = zeros(NSTEP + 1)
ys = zeros(NSTEP + 1)
sp_max = 1.0
sp_min = 1.0
# the angle v turns through in one step, summed over the run
turned = 0.0
for n to NSTEP {
    ox = vel[0]
    oy = vel[1]
    boris_step(vel, 0.0, 0.0, DT)
    # v turns clockwise for a positive charge in B = z-hat, so the cross product is negative
    turned = turned - atan2(ox * vel[1] - oy * vel[0], ox * vel[0] + oy * vel[1])
    xs[n + 1] = xs[n] + DT * vel[0]
    ys[n + 1] = ys[n] + DT * vel[1]
    sp = sqrt(vel[0] * vel[0] + vel[1] * vel[1])
    sp_max = max(sp_max, sp)
    sp_min = min(sp_min, sp)
}
show sp_max - sp_min

#| The center (x, y) of the circle through three points, from the perpendicular bisectors.
def circumcenter(ax, ay, bx, by, cx, cy) {
    d = 2.0 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    a2 = ax * ax + ay * ay
    b2 = bx * bx + by * by
    c2 = cx * cx + cy * cy
    ux = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / d
    uy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / d
    return [ux, uy]
}
center = circumcenter(xs[0], ys[0], xs[1], ys[1], xs[2], ys[2])
r_min = 1.0e30
r_max = 0.0
for n to NSTEP + 1 {
    r = sqrt((xs[n] - center[0]) * (xs[n] - center[0]) + (ys[n] - center[1]) * (ys[n] - center[1]))
    r_min = min(r_min, r)
    r_max = max(r_max, r)
}
show center
show r_min
show r_max

test "the speed stays 1 to rounding over 20000 steps" {
    assert_lt(sp_max - sp_min, 1.0e-13)
}
test "every position lies on one circle of radius sqrt(1 + dt^2/4) = 1.000312" {
    assert_near(r_min, sqrt(1.0 + DT * DT / 4.0), 1.0e-12)
    assert_near(r_max, sqrt(1.0 + DT * DT / 4.0), 1.0e-12)
}

vel = [0.0, 1.0]
xs = zeros(NSTEP + 1)
ys = zeros(NSTEP + 1)
sp_max = 1.0
sp_min = 1.0
# the angle v turns through in one step, summed over the run
turned = 0.0
for n to NSTEP {
    ox = vel[0]
    oy = vel[1]
    boris_step(vel, 0.0, 0.0, DT)
    # v turns clockwise for a positive charge in B = z-hat, so the cross product is negative
    turned = turned - atan2(ox * vel[1] - oy * vel[0], ox * vel[0] + oy * vel[1])
    xs[n + 1] = xs[n] + DT * vel[0]
    ys[n + 1] = ys[n] + DT * vel[1]
    sp = sqrt(vel[0] * vel[0] + vel[1] * vel[1])
    sp_max = max(sp_max, sp)
    sp_min = min(sp_min, sp)
}
show sp_max - sp_min

#| The center (x, y) of the circle through three points, from the perpendicular bisectors.
def circumcenter(ax, ay, bx, by, cx, cy) {
    d = 2.0 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    a2 = ax * ax + ay * ay
    b2 = bx * bx + by * by
    c2 = cx * cx + cy * cy
    ux = (a2 * (by - cy) + b2 * (cy - ay) + c2 * (ay - by)) / d
    uy = (a2 * (cx - bx) + b2 * (ax - cx) + c2 * (bx - ax)) / d
    return [ux, uy]
}
center = circumcenter(xs[0], ys[0], xs[1], ys[1], xs[2], ys[2])
r_min = 1.0e30
r_max = 0.0
for n to NSTEP + 1 {
    r = sqrt((xs[n] - center[0]) * (xs[n] - center[0]) + (ys[n] - center[1]) * (ys[n] - center[1]))
    r_min = min(r_min, r)
    r_max = max(r_max, r)
}
show center
show r_min
show r_max

test "the speed stays 1 to rounding over 20000 steps" {
    assert_lt(sp_max - sp_min, 1.0e-13)
}
test "every position lies on one circle of radius sqrt(1 + dt^2/4) = 1.000312" {
    assert_near(r_min, sqrt(1.0 + DT * DT / 4.0), 1.0e-12)
    assert_near(r_max, sqrt(1.0 + DT * DT / 4.0), 1.0e-12)
}

The period is too long by $(ω_{c} Δ t)^{2} /12$

The rotation in the step turns $v$ by $θ = 2 arctan (ω_{c} Δ t /2)$ , a little less than the exact $ω_{c} Δ t$ . The cell above summed the angle between successive velocities; dividing by the number of steps gives $θ$ , and the gyroperiod is $T = 2 π Δ t / θ$ . Expanding the arctangent,

\frac{T}{2 π / ω _{c}} = \frac{ω _{c} Δ t}{2 arctan ( ω _{c} Δ t /2 )} = 1 + \frac{( ω _{c} Δ t ) ^{2}}{12} + O (Δ t^{4}) .

At $Δ t = 0.05$ that is $T = 6.28449$ against $2 π = 6.28319$ . The speed is exact and the radius is off by $3 \times 1 0^{- 4}$ ; the phase falls behind by $2 \times 1 0^{- 4}$ of a turn per turn. The figure repeats the measurement for five time steps. At $Δ t = 0.4$ the $Δ t^{4}$ term shows: the error is $1.3194 \times 1 0^{- 2}$ , below $0. 4^{2} /12 = 1.3333 \times 1 0^{- 2}$ .

[3]lines 134–197

theta = turned / NSTEP
T_boris = 2.0 * PI * DT / theta
show theta
show T_boris

test "the measured gyroperiod is 2 pi dt / (2 atan(dt/2)) = 6.284494" {
    assert_near(T_boris, 2.0 * PI * DT / (2.0 * atan(DT / 2.0)), 1.0e-10)
    assert_near(T_boris, 6.28449408602382, 1.0e-10)
}
test "it exceeds 2 pi by the fraction dt^2/12, to within 1% of that fraction" {
    assert_near((T_boris / (2.0 * PI) - 1.0) / (DT * DT / 12.0), 1.0, 0.01)
}

# one and a quarter turns of the orbit, and the exact circle of radius 1 about the same center
N_TURN = 158
orb_x = zeros(N_TURN)
orb_y = zeros(N_TURN)
for n to N_TURN {
    orb_x[n] = xs[n]
    orb_y[n] = ys[n]
}
N_CIRC = 121
ang = linspace(0.0, 2.0 * PI, N_CIRC)
cxs = zeros(N_CIRC)
cys = zeros(N_CIRC)
for n to N_CIRC {
    cxs[n] = center[0] + cos(ang[n])
    cys[n] = center[1] + sin(ang[n])
}
figure("Boris gyro-orbit, dt = 0.05: the first 1.25 turns", "x", "y")
aspect("equal")
plot(orb_x, orb_y, "Boris positions")
plot(cxs, cys, "exact orbit, radius $v/\\omega_c = 1$")
style("exact orbit, radius $v/\\omega_c = 1$", "dash")

# the period error at five time steps: each run measures its rotation angle over 200 steps
DTS = [0.4, 0.2, 0.1, 0.05, 0.025]
N_MEAS = 200
per_err = zeros(len(DTS))
per_lead = zeros(len(DTS))
for k to len(DTS) {
    h = DTS[k]
    w = [0.0, 1.0]
    tot = 0.0
    for n to N_MEAS {
        ax0 = w[0]
        ay0 = w[1]
        boris_step(w, 0.0, 0.0, h)
        tot = tot - atan2(ax0 * w[1] - ay0 * w[0], ax0 * w[0] + ay0 * w[1])
    }
    per_err[k] = h / (tot / N_MEAS) - 1.0
    per_lead[k] = h * h / 12.0
}
show per_err
figure("Gyroperiod error of the Boris push", "time step $\\omega_c\\,\\Delta t$", "$T/(2\\pi/\\omega_c) - 1$")
xscale("log")
yscale("log")
plot(DTS, per_err, "measured")
plot(DTS, per_lead, "$(\\omega_c\\,\\Delta t)^2/12$")
style("$(\\omega_c\\,\\Delta t)^2/12$", "dash")

test "at dt = 0.4 the period error is 1.3194e-2, the exact arctan value, not 0.4^2/12 = 1.3333e-2" {
    assert_near(per_err[0], 1.319402e-02, 1.0e-7)
}

theta = turned / NSTEP
T_boris = 2.0 * PI * DT / theta
show theta
show T_boris

test "the measured gyroperiod is 2 pi dt / (2 atan(dt/2)) = 6.284494" {
    assert_near(T_boris, 2.0 * PI * DT / (2.0 * atan(DT / 2.0)), 1.0e-10)
    assert_near(T_boris, 6.28449408602382, 1.0e-10)
}
test "it exceeds 2 pi by the fraction dt^2/12, to within 1% of that fraction" {
    assert_near((T_boris / (2.0 * PI) - 1.0) / (DT * DT / 12.0), 1.0, 0.01)
}

# one and a quarter turns of the orbit, and the exact circle of radius 1 about the same center
N_TURN = 158
orb_x = zeros(N_TURN)
orb_y = zeros(N_TURN)
for n to N_TURN {
    orb_x[n] = xs[n]
    orb_y[n] = ys[n]
}
N_CIRC = 121
ang = linspace(0.0, 2.0 * PI, N_CIRC)
cxs = zeros(N_CIRC)
cys = zeros(N_CIRC)
for n to N_CIRC {
    cxs[n] = center[0] + cos(ang[n])
    cys[n] = center[1] + sin(ang[n])
}
figure("Boris gyro-orbit, dt = 0.05: the first 1.25 turns", "x", "y")
aspect("equal")
plot(orb_x, orb_y, "Boris positions")
plot(cxs, cys, "exact orbit, radius $v/\\omega_c = 1$")
style("exact orbit, radius $v/\\omega_c = 1$", "dash")

# the period error at five time steps: each run measures its rotation angle over 200 steps
DTS = [0.4, 0.2, 0.1, 0.05, 0.025]
N_MEAS = 200
per_err = zeros(len(DTS))
per_lead = zeros(len(DTS))
for k to len(DTS) {
    h = DTS[k]
    w = [0.0, 1.0]
    tot = 0.0
    for n to N_MEAS {
        ax0 = w[0]
        ay0 = w[1]
        boris_step(w, 0.0, 0.0, h)
        tot = tot - atan2(ax0 * w[1] - ay0 * w[0], ax0 * w[0] + ay0 * w[1])
    }
    per_err[k] = h / (tot / N_MEAS) - 1.0
    per_lead[k] = h * h / 12.0
}
show per_err
figure("Gyroperiod error of the Boris push", "time step $\\omega_c\\,\\Delta t$", "$T/(2\\pi/\\omega_c) - 1$")
xscale("log")
yscale("log")
plot(DTS, per_err, "measured")
plot(DTS, per_lead, "$(\\omega_c\\,\\Delta t)^2/12$")
style("$(\\omega_c\\,\\Delta t)^2/12$", "dash")

test "at dt = 0.4 the period error is 1.3194e-2, the exact arctan value, not 0.4^2/12 = 1.3333e-2" {
    assert_near(per_err[0], 1.319402e-02, 1.0e-7)
}

The $E \times B$ drift in crossed fields

Now add $E = 0.1 x^$ . The drift velocity is $v_{E} = E \times B / B^{2} = (0, - 0.1)$ . Write one Boris step as $v \mapsto R (v + a) + a$ , with $a = E Δ t /2$ and $R$ the rotation. The vector $v_{E} + a = (0.0025, - 0.1)$ is turned by $R$ into $(- 0.0025, - 0.1) = v_{E} - a$ , because its half-angle from $- y^$ is $arctan (0.0025/0.1) = arctan (Δ t /2) = θ /2$ . So $v_{E}$ is a fixed point of the step: a particle started there never gyrates.

For any other start, $u = v - v_{E}$ is only rotated, so the particle gyrates about a guiding center moving at $v_{E}$ . Averaged over $N$ steps the velocity differs from $v_{E}$ by the unfinished part of a turn, $∣ u_{0} ∣ ∣ sin (N θ /2) ∣/ (N sin (θ /2))$ , which is $8.2 \times 1 0^{- 4}$ here and falls as $1/ N$ .

[4]lines 217–273

EX = 0.1
VE_Y = -EX
fixed = [0.0, VE_Y]
for n to NSTEP {
    boris_step(fixed, EX, 0.0, DT)
}
show fixed

vd = [0.0, 1.0]
dx = 0.0
dy = 0.0
# about 5 gyroperiods of the trajectory for the figure
N_TRAJ = 630
tx = zeros(N_TRAJ)
ty = zeros(N_TRAJ)
for n to NSTEP {
    boris_step(vd, EX, 0.0, DT)
    dx = dx + DT * vd[0]
    dy = dy + DT * vd[1]
    if n < N_TRAJ - 1 {
        tx[n + 1] = dx
        ty[n + 1] = dy
    }
}
mean_vx = dx / (NSTEP * DT)
mean_vy = dy / (NSTEP * DT)
show mean_vx
show mean_vy
drift_err = sqrt(mean_vx * mean_vx + (mean_vy - VE_Y) * (mean_vy - VE_Y))
show drift_err

test "started at E x B / B^2 = (0, -0.1), the velocity stays there to rounding" {
    assert_near(fixed[0], 0.0, 1.0e-15)
    assert_near(fixed[1], VE_Y, 1.0e-15)
}
test "the mean velocity misses E x B / B^2 by |u0| |sin(N theta/2)| / (N sin(theta/2)) = 8.2167e-4" {
    th = 2.0 * atan(DT / 2.0)
    u0 = 1.0 - VE_Y
    assert_near(drift_err, u0 * abs(sin(NSTEP * th / 2.0)) / (NSTEP * sin(th / 2.0)), 1.0e-12)
    assert_near(drift_err, 8.216734535881e-4, 1.0e-12)
}

# the guiding center: starts one gyroradius |u0| / omega_c = 1.1 from the particle, moves at v_E
gcy = zeros(N_TRAJ)
gcx = zeros(N_TRAJ)
for n to N_TRAJ {
    gcx[n] = 1.0 - VE_Y
    gcy[n] = VE_Y * n * DT
}
figure("E x B drift: E = 0.1 x-hat, B = z-hat, 5 gyroperiods", "x", "y")
aspect("equal")
plot(tx, ty, "Boris trajectory")
plot(gcx, gcy, "guiding center at $E \\times B/B^2$")
style("guiding center at $E \\times B/B^2$", "dash")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

EX = 0.1
VE_Y = -EX
fixed = [0.0, VE_Y]
for n to NSTEP {
    boris_step(fixed, EX, 0.0, DT)
}
show fixed

vd = [0.0, 1.0]
dx = 0.0
dy = 0.0
# about 5 gyroperiods of the trajectory for the figure
N_TRAJ = 630
tx = zeros(N_TRAJ)
ty = zeros(N_TRAJ)
for n to NSTEP {
    boris_step(vd, EX, 0.0, DT)
    dx = dx + DT * vd[0]
    dy = dy + DT * vd[1]
    if n < N_TRAJ - 1 {
        tx[n + 1] = dx
        ty[n + 1] = dy
    }
}
mean_vx = dx / (NSTEP * DT)
mean_vy = dy / (NSTEP * DT)
show mean_vx
show mean_vy
drift_err = sqrt(mean_vx * mean_vx + (mean_vy - VE_Y) * (mean_vy - VE_Y))
show drift_err

test "started at E x B / B^2 = (0, -0.1), the velocity stays there to rounding" {
    assert_near(fixed[0], 0.0, 1.0e-15)
    assert_near(fixed[1], VE_Y, 1.0e-15)
}
test "the mean velocity misses E x B / B^2 by |u0| |sin(N theta/2)| / (N sin(theta/2)) = 8.2167e-4" {
    th = 2.0 * atan(DT / 2.0)
    u0 = 1.0 - VE_Y
    assert_near(drift_err, u0 * abs(sin(NSTEP * th / 2.0)) / (NSTEP * sin(th / 2.0)), 1.0e-12)
    assert_near(drift_err, 8.216734535881e-4, 1.0e-12)
}

# the guiding center: starts one gyroradius |u0| / omega_c = 1.1 from the particle, moves at v_E
gcy = zeros(N_TRAJ)
gcx = zeros(N_TRAJ)
for n to N_TRAJ {
    gcx[n] = 1.0 - VE_Y
    gcy[n] = VE_Y * n * DT
}
figure("E x B drift: E = 0.1 x-hat, B = z-hat, 5 gyroperiods", "x", "y")
aspect("equal")
plot(tx, ty, "Boris trajectory")
plot(gcx, gcy, "guiding center at $E \\times B/B^2$")
style("guiding center at $E \\times B/B^2$", "dash")

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

Gyration keeps the speed and the circle

The period is too long by (ωc​Δt)2/12

The E×B drift in crossed fields

The period is too long by $(ω_{c} Δ t)^{2} /12$

The $E \times B$ drift in crossed fields