Magnus expansion terms to any order

The numerical core writes out orders 1 to 6 of the Magnus expansion explicitly in Python, which is fast and keeps them checkable by eye. It generates orders 7 to 10 from the recursion. magnus.expansionterms derives the same terms symbolically, in exact rational arithmetic, at any order.

The symbolic derivation serves two purposes. It lets the hard-coded coefficients be checked rather than trusted — the test suite regenerates them and compares — and it lets you inspect an order beyond the implemented ceiling without running a propagation.

The recursion

Writing \(A(l) = -i \mathbb{H}(l)\), the Magnus expansion sums terms \(\Omega = \sum_n \Omega_n\) generated by the Bernoulli-number recursion [BCOR09], in the \(B_1 = -1/2\) convention:

\[\begin{split}\Omega_1(l) &= \int_{l_0}^{l} A(s)\, ds \\ \Omega_n(l) &= \sum_{j=1}^{n-1} \frac{B_j}{j!} \int_{l_0}^{l} S_n^{(j)}(s)\, ds ,\end{split}\]

with the \(S_n^{(j)}\) defined recursively in terms of the lower-order \(\Omega_m\):

\[\begin{split}S_n^{(1)} &= [\Omega_{n-1}, A] \\ S_n^{(j)} &= \sum_{m=1}^{n-j} [\Omega_m, S_{n-m}^{(j-1)}] , \qquad 2 \leq j \leq n-1 .\end{split}\]

Two facts about this recursion matter in practice.

Odd Bernoulli numbers vanish. The Bernoulli number \(B_j\) vanishes for every odd \(j \geq 3\), so whole commutator groups drop out and only \(j = 1, 2, 4, 6, 8, \ldots\) contribute. The surviving coefficients are the ones the numerical core hard-codes:

\[\frac{B_1}{1!} = -\frac{1}{2} , \quad \frac{B_2}{2!} = \frac{1}{12} , \quad \frac{B_4}{4!} = -\frac{1}{720} , \quad \frac{B_6}{6!} = \frac{1}{30240} , \quad \frac{B_8}{8!} = -\frac{1}{1209600} .\]

Every term is a right-nested chain. Unrolling the recursion shows that each term of \(\Omega_n\) has the form

\[[\Omega_{m_1}, [\Omega_{m_2}, \ldots [\Omega_{m_j}, A] \ldots ]] , \qquad m_1 + m_2 + \cdots + m_j = n - 1 ,\]

carrying the coefficient \(B_j/j!\). The terms of the \(j\)-th group are indexed by the compositions of \(n-1\) into \(j\) positive parts, of which there are \(\binom{n-2}{j-1}\). Hence, the number of terms is

\[\#\,\Omega_n = \sum_{j\,:\,B_j \neq 0} \binom{n-2}{j-1} ,\]

which is what the implementation for orders 7 to 10 iterates over directly, rather than storing their 244 expressions.

How many terms there are

import magnus.expansionterms as et

{n: et.count_terms(n) for n in range(1, 13)}

{1: 1,
 2: 1,
 3: 2,
 4: 3,
 5: 5,
 6: 9,
 7: 17,
 8: 33,
 9: 65,
 10: 129,
 11: 257,
 12: 513}

The count roughly doubles per order. The implemented ceiling, magnus.magnus.MAGNUS_EXP_ORDER_MAX = 10, is chosen for this reason: the terms remain easy to generate far beyond it, but the work per slab grows with their number.

The expansion, written out

The numerical core spells out orders 1 through 6 inline:

et.print_magnus_terms(6)

Omega_1  (1 term)
    int A

Omega_2  (1 term)
    int -1/2 [Om_1, A]

Omega_3  (2 terms)
    int -1/2 [Om_2, A]
    int +1/12 [Om_1, [Om_1, A]]

Omega_4  (3 terms)
    int -1/2 [Om_3, A]
    int +1/12 [Om_1, [Om_2, A]]
    int +1/12 [Om_2, [Om_1, A]]

Omega_5  (5 terms)
    int -1/2 [Om_4, A]
    int +1/12 [Om_1, [Om_3, A]]
    int +1/12 [Om_2, [Om_2, A]]
    int +1/12 [Om_3, [Om_1, A]]
    int -1/720 [Om_1, [Om_1, [Om_1, [Om_1, A]]]]

Omega_6  (9 terms)
    int -1/2 [Om_5, A]
    int +1/12 [Om_1, [Om_4, A]]
    int +1/12 [Om_2, [Om_3, A]]
    int +1/12 [Om_3, [Om_2, A]]
    int +1/12 [Om_4, [Om_1, A]]
    int -1/720 [Om_1, [Om_1, [Om_1, [Om_2, A]]]]
    int -1/720 [Om_1, [Om_1, [Om_2, [Om_1, A]]]]
    int -1/720 [Om_1, [Om_2, [Om_1, [Om_1, A]]]]
    int -1/720 [Om_2, [Om_1, [Om_1, [Om_1, A]]]]

Orders 7 through 10 are generated from the recursion rather than typed in. Here is order 7 in full, and the sizes of the rest:

for term in et.omega_terms(7):
    print("   int " + et.format_term(term))

   int -1/2 [Om_6, A]
   int +1/12 [Om_1, [Om_5, A]]
   int +1/12 [Om_2, [Om_4, A]]
   int +1/12 [Om_3, [Om_3, A]]
   int +1/12 [Om_4, [Om_2, A]]
   int +1/12 [Om_5, [Om_1, A]]
   int -1/720 [Om_1, [Om_1, [Om_1, [Om_3, A]]]]
   int -1/720 [Om_1, [Om_1, [Om_2, [Om_2, A]]]]
   int -1/720 [Om_1, [Om_1, [Om_3, [Om_1, A]]]]
   int -1/720 [Om_1, [Om_2, [Om_1, [Om_2, A]]]]
   int -1/720 [Om_1, [Om_2, [Om_2, [Om_1, A]]]]
   int -1/720 [Om_1, [Om_3, [Om_1, [Om_1, A]]]]
   int -1/720 [Om_2, [Om_1, [Om_1, [Om_2, A]]]]
   int -1/720 [Om_2, [Om_1, [Om_2, [Om_1, A]]]]
   int -1/720 [Om_2, [Om_2, [Om_1, [Om_1, A]]]]
   int -1/720 [Om_3, [Om_1, [Om_1, [Om_1, A]]]]
   int +1/30240 [Om_1, [Om_1, [Om_1, [Om_1, [Om_1, [Om_1, A]]]]]]

{n: et.count_terms(n) for n in (8, 9, 10)}

{8: 33, 9: 65, 10: 129}

Going beyond the implemented ceiling

The generator has no ceiling of its own. Order 12, which the numerical core does not implement, is still inspectable:

terms_12 = et.omega_terms(12)
print(f"Omega_12 has {len(terms_12)} terms; the first three are")
for term in terms_12[:3]:
    print("   int " + et.format_term(term))

Omega_12 has 513 terms; the first three are
   int -1/2 [Om_11, A]
   int +1/12 [Om_1, [Om_10, A]]
   int +1/12 [Om_2, [Om_9, A]]

Coefficients come back as fractions.Fraction, so they are exact at any order and can be compared without a tolerance:

from fractions import Fraction

et.bernoulli_factor(8) == Fraction(-1, 1209600)

True

Checking the implementation against the derivation

This module shares no code with the numerical core, so agreement between the two is an independent check. The hard-coded group factors agree with the derived ones:

import magnus.magnus as mg

{'F1': (mg.F1, float(et.bernoulli_factor(2))),
 'F2': (mg.F2, float(et.bernoulli_factor(4))),
 'F3': (mg.F3, float(et.bernoulli_factor(6))),
 'F4': (mg.F4, float(et.bernoulli_factor(8)))}

{'F1': (0.08333333333333333, 0.08333333333333333),
 'F2': (-0.001388888888888889, -0.001388888888888889),
 'F3': (3.306878306878307e-05, 3.306878306878307e-05),
 'F4': (-8.267195767195768e-07, -8.267195767195768e-07)}

tests/test_expansionterms.py goes further: it evaluates the generated terms numerically on a sampled \(A(l)\) and compares them, order by order, against what magnus.magnus.magnus_expansion() produces internally. They agree to a relative precision of 1e-11 at every order from 1 to 10, the accuracy of the quadrature that evaluates them. This single check covers both the hand-written low orders and the generated high ones.

Choosing an order

A higher order makes the error fall faster with the slab width. The table below was measured on a single slab with a smooth three-flavor Hamiltonian, against a tight-tolerance ODE solution. It uses the cumulative-quadrature method 'simpson', for which magnus_exp_order is the index of the last term kept and the delivered order is higher (see the table in Methodology). The error of a single slab carries one more power of \(h\) than the error over the whole trajectory: order 4, which delivers order 6 over a trajectory, converges as \(h^7\) on one slab.

Magnus order

Observed convergence rate

Cost per slab, relative to order 6

4

\(\sim h^{7.1}\)

0.25

6

\(\sim h^{9.1}\)

1

7

\(\sim h^{9.1}\)

2.7

8

\(\sim h^{11.3}\)

5.0

10

measurement floor

16.9

Two caveats keep this from being a simple “higher is better”.

First, the cost grows faster than the order does, and narrowing the slabs at order 4 or 6 often reaches a given accuracy for less total work than raising the order. Magνs warns (magnus.magnus.MagnusHighOrderCostWarning) when an order above 6 is requested on the quadrature path, for exactly this reason. The collocation schemes are unaffected: order 8 there costs four Hamiltonian evaluations and six commutators.

Second, and more fundamentally, the Magnus series is guaranteed to converge only while \(\int \lVert A \rVert\, dl < \pi\). This bound is sufficient, not necessary. Beyond it, a higher order is not guaranteed to help, and the remedy is narrower slabs — see Methodology and magnus.magnus.MagnusConvergenceWarning.

Note

Odd orders and orders above 8 require integration_method='trapezoid' or 'simpson'. The Gauss–Legendre collocation schemes ('gl', the default) are separately derived integrators [BCR00, BCR02], not products of this recursion, and exist only for orders 2, 4, 6 and 8. Any other order raises a ValueError rather than quietly returning a result of a different order.

API

See magnus.expansionterms in the API reference for the full signatures.