# Novel frame changes for quantum physics

Pierre-Louis Giscard, Omid Faizy, Christian Bonhomme

► To cite this version:

Pierre-Louis Giscard, Omid Faizy, Christian Bonhomme. Novel frame changes for quantum physics. 2026. hal-05470150

HAL Id: hal-05470150

https://hal.sorbonne-universite.fr/hal-05470150v1

Preprint submitted on 21 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Novel frame changes for quantum physics

Pierre-Louis Giscard†,*, Omid Faizy†,*, and Christian Bonhomme*

†Laboratoire de Mathématiques Pures et Appliquées Joseph Liouville,
Université du Littoral Côte d'Opale, 50 rue Ferdinand Buisson, CS 80699, 62228 Calais, France. and

*Laboratoire de Chimie de la Matière Condensée de Paris, UMR CNRS 7574,
Sorbonne Université, 4, place Jussieu, 75252 Paris, France.

(Dated: October 7, 2025)

We present novel, exotic types of frame changes for the calculation of quantum evolution operators. We detail in particular the biframe, in which a physical system's evolution is seen in an equal mixture of two different standard frames at once. We prove that, in the biframe, convergence of all series expansions of the solution is quadratically faster than in 'conventional' frames. That is, if in laboratory frame or after a standard frame change the error at order n of some perturbative series expansion of the evolution operator is on the order of $\epsilon^n$, $0 < \epsilon < 1$, for a computational cost $C(n)$ then it is on the order of $\epsilon^{2n+1}$ in the biframe for the same computational cost. We demonstrate that biframe is one of an infinite family of novel frames, some of which lead to higher accelerations but require more computations to set up initially, leading to a trade-off between acceleration and computational burden.

# I. INTRODUCTION

In this work we consider physical systems whose properties or dynamics is described by a system of ordinary linear differential equations with possibly time-dependent coefficients,

$$ \frac{d}{dt}U(t) = A(t)U(t), \quad U(0) = \mathrm{Id}. \qquad (1) $$

Chief examples include closed quantum systems obeying Schrödinger equation with $A(t) = -iH(t)$, $H(t)$ being the quantum Hamiltonian; as well as the partition function $Z(\beta)$ of statistical physics, which being the exponential of the Hamiltonian operator $\exp(-\beta H)$, solves the ordinary differential equation $\dot{Z}(\beta) = -\mathcal{H}Z(\beta)$.

Calculating analytically the solution $U$ (called evolution operator) of Eq. (1) when $A(t)$ does not commute with itself at different times is notoriously difficult. Formally, the solution is designated as a time-ordered exponential,

$$ \begin{aligned} U(t) &= \mathcal{T}e^{\int_0^t A(\tau)d\tau} \\ &= \mathrm{Id} + \int_0^t A(\tau)d\tau + \int_0^t \int_0^{\tau_1} A(\tau_1)A(\tau_2)d\tau_2d\tau_1 + \dots \end{aligned} $$

with $\mathcal{T}$ the time-ordering operator. The lack of self-commutativity $A(t)A(s) \neq A(s)A(t)$ hinders simplifications in the Dyson series above and leads to a break down in the invariance of the evolution operator under time translations $U(t,s) \neq U(t-s,0)$. That is, evolving the solution from a time $s$ to a time $t > s$ is different from evolving it from $0$ to $t-s$. This critical observation indicates that the correct general mathematical framework to determine $U$ is that of bivariate functions and matrices, i.e. depending on two time-variables, even though only $U(t) := U(t,0)$ is usually desired in physics. Volterra and Pérès were the first to consider the algebraic structures necessitated by this two-times approach to differential equations in a pioneering study completed in the early 1920s [15]. Lacking a proper

understanding of distributions at the time, their mathematics encountered many profound obstacles and their work was subsequently overlooked. In a series of recent developments the operation Volterra introduced [14] to deal with iterated integrals (the Volterra composition) has been given its proper context and generalization within the theory of distributions, yielding the $\star$ product,

$$ (f \star g)(t, s) := \int_{-\infty}^{+\infty} f(t, \tau)g(\tau, s)d\tau. \qquad (2) $$

Here $f,g \in \mathcal{D}$ belong to a certain set of distributions which we do not need to explicit, see [12] for a rigorous presentation and Appendix A. The $\star$ product extends naturally to matrices in $A, B \in \mathcal{D}^{N \times N}$ with

$$ (A \star B)(t, s) = \int_{-\infty}^{+\infty} A(t, \tau)B(\tau, s)d\tau. \qquad (3) $$

This product has a unit, $I_{\star} := \mathrm{Id} \delta(t-s)$. For our purposes here the most important fruit of this formalism is that even if $A(t)$ does not commute with itself at all times the differential system of Eq. (1) has Green's function

$$ G(t, s) = (I_{\star} - A(t)\Theta(t - s))^{*-1}, \qquad (4) $$

and $U(t,s) = \int_s^t G(\tau, s)d\tau$ is the integral of $G$. In this expression, $\Theta(t-s)$ is the Heaviside theta function, $\Theta(x) = 1$ if $x \ge 0$ and 0 otherwise. Interestingly, it arises out of mathematical necessity but encodes a physical reality: causality, $s$ being the time at which initial conditions are imposed and the drive by $A$ starts, while $t$ is the time of observation. The crucial observation is that $G$ is a true resolvent of $A$: this means that it is amenable to all techniques from ordinary linear algebra, so long as $\star$ products are substituted in place of ordinary matrix products.

An implication of this observation concerns frame changes in physics. By frame change, we here mean in a narrow sense those implemented mathematically on differential systems like Eq. (1). Concretely it is well known that the evolution operator $U$ solution of Eq. (1) can be recast as

$$ U(t) = U_B(t) \mathcal{T} e^{\int_0^t U_B^\dagger(\tau) A(\tau) U_B(\tau) d\tau}, $$

<footnote>* giscard@univ-littoral.fr</footnote>
with, for discrete systems, $B(t)$ a matrix of the same size as $A$ and $U_B$ the evolution operator solution of $\dot{U}_B = B \cdot U$. This corresponds to seeing the dynamics driven by $A$ from a new frame whose movement is due to $B$. We shall see that this technique is in fact the simplest of infinitely many novel types of frame changes which offer progressively higher computational rewards as they become increasingly intricate.

## II. FRAME CHANGES AS LINEAR ★-ALGEBRA

### A. Ordinary frame change

Instead of directly aiming for frame changes, let us consider the following problem: suppose that the matrix $A(t)$ can be expressed as the sum of two parts

$$A(t) = A_0(t) + A_1(t)$$

and that one wishes to relate the evolution operator $U$ solution of $\dot{U} = A(t)U$ to the evolution operators $U_0$ and $U_1$ solving the equations $\dot{U}_0 = A_0(t)U_0$ and $\dot{U}_1 = A_1(t)U_1$, respectively. Given that the Green's function $G$ is a ★-resolvent as shown by Eq. (4), the question is the same as expressing, for two matrices $M_0$ and $M_1$, the quantity $R := (I - M_0 - M_1)^{-1}$ in terms of the 'isolated' resolvents $R_0 := (I - M_0)^{-1}$ and $R_1 := (I - M_1)^{-1}$. Standard linear algebra provides many such representations. For example,

$$\frac{1}{I - M_0 - M_1} = \frac{1}{I - M_1} \frac{1}{I - M_0 \frac{1}{I - M_1}}, \quad (5a)$$

which we wrote in fraction form for improved readability but which should be understood as the statement that

$$R = R_1 \cdot (I - M_0 \cdot R_1)^{-1}. \quad (5b)$$

This equation remains valid should $M_0$ and $M_1$ not commute.

Let us see what this implies for Green's functions. Turning back to $A(t) = A_0(t) + A_1(t)$, let $G_1 = (I_* - A_0 \Theta)^{-1}$ designate the Green's function associated to $A_1$. Then Eq. (5b) is equivalent to the statement that

$$G = G_1 \star (I_* - A_0 \Theta \star G_1)^{-1},$$

and the full evolution operator $U(t)$ obeys,

$$U = U_1 \star (I_* - A_0 \Theta \star G_1)^{-1}. \quad (6)$$

Evaluating the ★-products and inverses explicitly in the above yields (Appendix B),

$$\mathcal{F}(t, s) = \int_s^t U_1^{-1}(\tau) A_0(\tau) U_1(\tau) d\tau, \quad (7a)$$

$$U(t, s) = U_1(t) \mathcal{T} e^{\mathcal{F}(t,s)} U_1^{-1}(s), \quad (7b)$$

which is the standard frame change formula. This proof of the formula is neither better nor simpler than the classical proofs relying on differential calculus, but the point is we obtained it from the purely linear algebraic statement

Eq. (5b). This suggests that changing frame is one of a very large reservoir of differential transformations which are ★-linear algebraic results in disguise. Indeed, plenty more such formulas must exist since resolvents of ordinary matrices satisfy a host of relations in the spirit of Eq. (5b), all of which ought to correspond to valid results concerning Green's functions and evolution operators.

We put this reasoning to the test by looking for an hitherto unknown types of frame change, in particular one that does not break the symmetric roles played by $A_0$ and $A_1$ in $A$: the biframe. As an additional bonus, we will show that contrary to the standard frame change, advanced frame changes such as the biframe intrinsically accelerate the convergence of perturbation series.

### B. Biframe change

Coming back to Eq. (5b), we see that the formula treats $M_0$ and $M_1$ differently. One might instead wish for a more symmetrical treatment of both parts that does not give undue importance to one over the other. This is achieved by the following formula (which is certainly not the only possible solution to the quest for symmetry!), which is just as easily verified using first year calculus rules:

$$\frac{1}{I - M_0 - M_1} = \frac{1}{I - M_0} \frac{1}{I - \frac{M_1}{I - M_1} \frac{M_0}{I - M_0}} \frac{1}{I - M_1}, \quad (8a)$$

more rigorously,

$$R = R_0 \cdot (I - M_1 \cdot R_1 \cdot M_0 \cdot R_0)^{-1} \cdot R_1. \quad (8b)$$

Again this remains true when $M_0$ and $M_1$ do not commute.

Let us now see what this implies for the evolution operator $U$. Taking $A(t) = A_0(t) + A_1(t)$, replacing matrix products by ★-products and ordinary resolvents $R$ by Green's functions one gets,

$$G = G_0 \star (I_* - A_1 \Theta \star G_1 \star A_0 \Theta \star G_0)^{-1} \star G_1, \quad (9)$$

and from there

$$U = U_0 \star (I_* - A_1 \Theta \star G_1 \star A_0 \Theta \star G_0)^{-1} \star G_1. \quad (10)$$

In spite of appearances with e.g. $A_1$ on the left of $A_0$ inside the ★-resolvent, both parts play exactly equivalent roles in these results. For example, noting that $\dot{U}_i = A_i \Theta \star G_i$, Eqs. (9, 10) are equivalent to the following result for the derivative $\dot{U}$ (Appendix C),

$$\dot{U} = \dot{U}_0 + \dot{U}_1 + \dot{U}_0 \star \dot{U}_1 + \dot{U}_1 \star \dot{U}_0 + \dot{U}_0 \star \dot{U}_1 \star \dot{U}_0 + \dot{U}_1 \star \dot{U}_0 \star \dot{U}_1 + \dot{U}_0 \star \dot{U}_1 \star \dot{U}_0 \star \dot{U}_1 + \dots,$$

which is manifestly invariant under the exchange of indices $0 \leftrightarrow 1$.

Evaluating the ★-products in the ★-resolvent and keeping in mind that such a resolvent is a time-ordered exponential, Eq. (10) is equivalent to

$$\mathcal{B}(t, s) := A_1(t)U_1(t) \int_s^t U_1^{-1}(\tau) A_0(\tau)U_0(\tau) d\tau U_0^{-1}(s), \quad (11a)$$

$$U = U_0 \star \mathcal{T} e^{\mathcal{B}(t,s)} \star G_1. \quad (11b)$$
This is the biframe change. A detailed proof of this result is presented in Appendix D. An alternative though completely equivalent form of the biframe is obtained on rearranging the terms in the formula above, yielding (proof in Appendix E),

$$ \mathcal{B}_2(t, s) := A_1(t)U_0(t) \int_s^t U_0^{-1}(\tau)A_0(\tau)U_1(\tau)d\tau U_1^{-1}(s), \quad (12a) $$

$$ U = U_0 \star G_1 \star \mathcal{T} e^{\mathcal{B}_2(t,s)}. \quad (12b) $$

In Appendix F we present simplifications of the above results should A be a constant matrix. Before we see how to implement the biframe in a concrete example, we show that it is worth the effort.

### C. Intrinsic acceleration of perturbation series

While standard frame changes can yield structurally simpler Hamiltonians, the frame change itself is not expected to inherently speed up perturbation series. To see this, coming back to Section II A, we showed that the standard frame change is equivalent to the ★-linear algebraic statement Eq. (6). Expanding this expression as a Dyson series shows that the mth order approximation $U_{\text{std-frame}}^{[m]}$ of U after a standard frame change defined as

$$ U_{\text{std-frame}}^{[m]} := U_1 \star \sum_{k=0}^{m} (A_0 \Theta \star G_1)^{*k}, $$

correctly reproduces up to order m of the original Dyson series expansion of U in the laboratory frame. In other terms, the change of frame did not alter the approximation order m and therefore did not accelerate the series expansion of U.

In contrast, the biframe change intrinsically produces an acceleration of this expansion. To see why, let us go back to the standard linear algebra. Suppose for this discussion that some matrix M has a spectral radius $\rho(M) < 1$ so that the Neumann series $\sum_{k=0}^{\infty} M^k$ converges. This series is a natural perturbative expansion scheme (with a small parameter $\rho(M)$) for the resolvent $R := (I - M)^{-1}$. Since $M^{k+1} = M \cdot M^k$, then at order $m > 1$ the standard approximation $R_{\text{std}}^{[m]} := \sum_{k=0}^{m} M^k$ can be calculated with $m - 1$ matrix products. At the same time, it is equally true that

$$ R = (I + M) \cdot (I - M^2)^{-1}, $$

so that another expansion scheme for the resolvent is $R^{[m]} := (I + M) \cdot \sum_{k=0}^{m} N^k$, with $N = M^2$. Calculating this requires $(m - 1) + 1 + 1 = m + 1$ matrix products, $m - 1$ of which come from $\sum_{k=0}^{m} N^k$, one more from calculating $N = M^2$ and another one for the final multiplication by $(I + M)$. Yet, expanding $R^{[m]}$ immediately shows that it reproduces up to order $2m + 1$ of the original Neumann series $R^{[m]} = R_{\text{std}}^{[2m+1]}$, roughly doubling the approximation order at the same computational cost.

Eq. (8b) is a peculiar instance of the above construction.

Observe that for $k \ge 0$, the expression $R_0(M_1 R_1 M_0 R_0)^k R_1$ comprises $2k + 1$ exchanges between the 0 and 1 indices,

precisely as much as $(M^2)^n = (M_0 + M_1)^k$. Since furthermore $R_i := (I - M_i)^{-1}$, $i = 0, 1$, then the truncated series

$$ R^{[m]} := R_0 \cdot \sum_{k=0}^{m} (M_1 \cdot R_1 \cdot M_0 \cdot R_0)^k \cdot R_1, $$

is verified to reproduce the first $2m + 1$ terms of the standard Neumann series while necessitating only $m + 1$ matrix products.

Given that solving N-ODEs is nothing but doing ★-linear algebra, these results hold true for Green's functions. Now the biframe operator $\mathcal{B}$ of Eq. (11a) involves the two parts $A_1$ and $A_0$ so that the Dyson series of its time-ordered exponential, if truncated at order $m$, necessitates $m + 1$ ★-products (i.e. iterated integrals with matrix products) while yielding an expression that is exact up to order $2m+1$ of the Dyson series expansion of U after a standard frame change. Denoting by $\epsilon^n$ the error at order $n$ of the Dyson expansion without frame change (or equally after a standard frame change), the error at the same order of expansion in the biframe is thus $\epsilon^{2n+1}$, quadratically better at nearly fixed computational cost. From an information point of view, this is because we assume that we exactly know the Green's functions of the isolated parts $G_0$ and $G_1$, the biframe change leveraging this knowledge into an acceleration of the original series.

### D. Further frame changes

Given the new perspective on frame changes as ★-linear algebraic statements on resolvents (like Eqs. 5, 8), there exists as many frame changes as there are ways to express the resolvent R of a matrix M := $\sum_i M_i$ in terms of its $M_i$ parts and their the resolvents. Many of these exotic frame changes can produce further accelerations depending on the known quantities or symmetries that one wishes to preserve. For instance, let us devise a trframe change where the systems is seen in three different frames at once, each generated by one part of an operator A(t) = $A_0(t) + A_1(t) + A_2(t)$. This is equivalent to asking for the ordinary resolvent R of a matrix M = $M_0 + M_1 + M_2$ in terms of the resolvents $R_i$, $i = 0, 1, 2$ in a way that these play equivalent roles. Standard linear algebra indicates that,

$$ R = R_0 \cdot (I - M_1 \cdot R_1 \cdot M_0 \cdot R_0)^{-1} \cdot R_1 \cdot (I - M_2 \cdot R_2 \cdot (M_0 \cdot R_0 \cdot (I - M_1 \cdot R_1 \cdot M_0 \cdot R_0)^{-1} \cdot R_1 + M_1 \cdot R_1 \cdot (I - M_0 \cdot R_0 \cdot M_1 \cdot R_1)^{-1} \cdot R_0))^{-1} \cdot R_2, $$

For evolution operators, this is therefore

$$ U = U_0 \star (I_\star - \dot{U}_1 \star \dot{U}_0)^{-1} \star G_1 \star (I_\star - \dot{U}_2 \star (\dot{U}_0 \star (I_\star - \dot{U}_1 \star \dot{U}_0)^{*-1} \star G_1 + \dot{U}_1 \star (I_\star - \dot{U}_0 \star \dot{U}_1)^{*-1} \star G_0))^{*-1} \star G_2 $$

where $\dot{U}_i := A_i(t) \Theta \star G_i$ and $G_i$ is the Green's function associated to Hamiltonian $H_i(t)$. Once expanded as a series
the triframe generates the expansion

$$ \dot{U} = \dot{U}_0 + \dot{U}_1 + \dot{U}_2 + \dot{U}_0 \star \dot{U}_1 + \dot{U}_1 \star \dot{U}_0 + \dot{U}_0 \star \dot{U}_2 + \dot{U}_2 \star \dot{U}_0 \\ + \dot{U}_1 \star \dot{U}_2 + \dot{U}_2 \star \dot{U}_1 + \dot{U}_0 \star \dot{U}_1 \star \dot{U}_2 + \dot{U}_0 \star \dot{U}_2 \star \dot{U}_1 + \dots $$

which is invariant under exchange of indices 0 ↔ 1 ↔ 2.

To determine the acceleration of series expansions in the triframe we follow the analysis of section II C. Let $L := M^3$, then observe that

$$ R = (I + M + M^2) \cdot (I - L)^{-1}, $$

which implies that $(I + M + M^2) \cdot \sum_{k=0}^{m} L^k$ is a series expansion of $R$ necessitating only $m+2$ matrix products to reach order $3m+2$ of the original Neumann series expansion of $R$ in powers of $M$. This implies that cutting the original matrix $M$ into three parts, each of which, if taken in isolation, has a known Green's function, is leveraged into a cubic acceleration of the Dyson series for $U$ in the triframe. Note that the triframe is equivalent to applying a biframe change twice, which is not particularly creative: many more ways of doing exist, e.g. if only some but not all isolated Green's functions $G_i$ are known or if one wants to single emphasize the role of one ordinary frame over the others. For example one could choose to mix a standard frame change with a biframe, leading to smaller acceleration of convergence as compared to the triframe but possibly easier to wield operators.

Clearly, frame changes tailored to any one problem can be devised by leveraging known quantities such as Green's functions of isolated parts into a useful expression with profitable acceleration for perturbative expansions.

### III. EXAMPLE

#### A. Acceleration of Dyson series for a time-dependent Hamiltonian

Consider the Schrödinger equation with Hamiltonian

$$ H(t) = \frac{\omega_0}{2} \sigma_z + 2\beta \cos(\omega t) \sigma_x. $$

Denote $H_0 = (\omega_0/2)\sigma_z$, $H_1 = 2\beta \cos(\omega t)\sigma_x$. In consequence, there are two natural standard frame changes for this Hamiltonian. The first is,

$$ H \rightarrow H_{\text{std},0}(t) := U_0^\dagger(t) \cdot H \cdot U_0(t), \\ = \begin{pmatrix} \frac{\omega_0}{2} & 2\beta e^{i\omega_0 t} \cos(\omega t) \\ 2\beta e^{-i\omega_0 t} \cos(\omega t) & -\frac{\omega_0}{2} \end{pmatrix}. $$

This frame change is mostly used to justify the rotating wave-approximation in near-resonant cases $\omega_0 \approx \omega$ with both large. The second frame change leads to

$$ H \rightarrow H_{\text{std},1}(t) := U_1^\dagger(t) \cdot H \cdot U_1(t), \\ = 2\beta \cos(\omega t) \sigma_x + \frac{\omega_0}{2} \sin\left(\frac{4\beta}{\omega} \sin(\omega t)\right) \sigma_y \\ + \frac{\omega_0}{2} \cos\left(\frac{4\beta}{\omega} \sin(\omega t)\right) \sigma_z $$

This choice is advantageous in the situation where $\omega \gg 1$, the time-average of $H_{\text{std},1}$ being the first order of the high-frequency expansion of $U(t)$ [8, 13]. Subsequent orders correspond to higher moments of $H_{\text{std},1}$. We now turn to the biframe, the driving operator is

$$ \mathcal{B}(t, s) = \beta\omega_0 \cos(\omega t) \cos\left(\frac{2\beta}{\omega} \sin(\omega t)\right) \begin{pmatrix} -iS(t, s) & \bar{C}(t, s) \\ -C(t, s) & i\bar{S}(t, s) \end{pmatrix} \\ + \beta\omega_0 \cos(\omega t) \sin\left(\frac{2\beta}{\omega} \sin(\omega t)\right) \begin{pmatrix} iC(t, s) & \bar{S}(t, s) \\ -S(t, s) & -i\bar{C}(t, s) \end{pmatrix}, $$

with

$$ S(t, s) := e^{\frac{1}{2}i\omega_0 s} \int_s^t e^{-\frac{1}{2}i\omega_0 \tau} \sin\left(\frac{2\beta}{\omega} \sin(\omega\tau)\right) d\tau, $$

$$ C(t, s) := e^{\frac{1}{2}i\omega_0 s} \int_s^t e^{-\frac{1}{2}i\omega_0 \tau} \cos\left(\frac{2\beta}{\omega} \sin(\omega\tau)\right) d\tau, $$

and $\bar{S}$, $\bar{C}$ are the conjugates of $S$ and of $C$, respectively. In order to demonstrate the acceleration of perturbative expansions that is intrinsic to the biframe, we present a numerical comparison of the Dyson expansion in the laboratory frame (no frame change), after the two standard frame changes and in the biframe. To that end, define the $m$th-order Dyson approximations

$$ \begin{align*} U_{\text{lab}}^{[m]} &:= \Theta \star \sum_{k=0}^{m} H^{*k}, \\ U_{\text{std-frame}}^{[m]} &:= U_1 \star \sum_{k=0}^{m} H_{\text{std}}^{*k}, \\ U_{\text{bframe}}^{[m]} &:= U_0 \star \sum_{k=0}^{m} (\mathcal{B}\Theta)^{*k} \star G_1, \end{align*} $$

where $H_{\text{std}}$ is the Hamiltonian after a standard frame change, i.e. either $H_{\text{std},0}$ or $H_{\text{std},1}$. Recall that for any matrix comprising smooth functions of time, $(A\Theta)^{*k}$ reduces to the $k$th iterated integral of $A(t)$ so that $\sum_k (A\Theta)^{*k}$ is the Dyson series $\mathcal{I}\delta(t-s) + A(t) + \int_s^t A(t)A(\tau)d\tau + \dots$ as claimed.

We quantify the error associated with each order $m$ approximation as follows. Let $U_r$ be the reference evolution operator as evaluated within machine precision (relative and absolute tolerances set to $10^{-16}$ for each entry) by a standard numerical solver (in this case Mathematica's NDSolve). We evaluate the accuracy of $U_{\text{lab}}^{[m]}$, $U_{\text{std-frame}}^{[m]}$ and $U_{\text{bframe}}^{[m]}$ by evaluating the deviation from 1 of their normalised Frobenius scalar products with $U_r$ over a total evolution time $T$,

$$ \epsilon := \frac{1}{T} \int_{0}^{T} \left( 1 - \frac{\text{Tr}(U_{r}^{\dagger}(\tau)U^{[m]}(\tau))}{\sqrt{\|U_{r}(\tau)\|_{F} \|U^{[m]}(\tau)\|_{F}}} \right) d\tau, \quad (13) $$

which is the relative error on $U^{[m]}$ with respect to the reference solution. Here $\|A\|_F := \text{Tr}(A^\dagger A)$ designates the Frobenius norm of matrix $A$. As constructed above, the relative error $\epsilon$ evaluates to 0 if and only if $U^{[m]}(t) = U(t)$ exactly at all times between $t=0$ and $t=T$. In Fig. 1 we show the
<figure><img src="image_2.png" /><figcaption>FIG. 1: Logarithm of the relative error $\epsilon$ of Eq. (13) for the first 12 orders of the Dyson series expansion in the laboratory frame (red circles and line), after a standard frame change (black squares and line) and in the biframe (blue stars and line). Scaled parameters: $\omega_0/\omega \simeq 0.67$, $\beta/\omega \simeq 0.53$, total scaled simulation time $\omega T = 6$. The dashed straight line at the bottom represents machine precision. The poor performance of the laboratory frame is due to the inability of its low orders to correctly fit the behavior at long times $t \sim T$, which heavily degrades the overall error measure $\epsilon$.</figcaption></figure>

relative errors in all three frames and up to order $m = 12$ of the Dyson series. Numerical computations of $\star$-products and $\star$-resolvents rely on trapezoidal quadrature as detailed in [1].

## IV. CONCLUSION

We show that changing frame in the context of autonomous and non-autonomous systems of coupled ordinary linear differential equations such a Schrödinger's, is an instance of a simple linear algebraic statement with respect to a generalized convolution-like product, known as the $\star$-product. Exploiting this observation, we showed that infinitely many novel frames exist and that these intrinsically accelerate perturbative expansions of evolution operators and partition functions.

The results presented here are certainly not limited to quantum mechanics nor to first order linear differential equations. Indeed, since any $n$th order non-autonomous linear differential equation is equivalent to a first order non-autonomous $n \times n$ system of coupled linear differential equations, then it is amenable to all the frame changes presented in this work. We shall illustrate this in a coming work with novel perturbative expansion of general Heun functions around hypergeometric solutions of Heun equations.

## Acknowledgements

This work is supported by the French National Research Agency ANR-20-CE29-0007 project MAGICA.

# Appendix A: Properties of the $\star$ product

## 1. Relation to Volterra composition

For bivariate functions $f$ and $g$ that are smooth in both of their arguments $t$ and $s$, the $\star$-product between $f(t, s)\Theta(t-s)$ and $g(t, s)\Theta(t-s)$ is the Volterra composition of the first kind $\star_v$ of the smooth functions $f$ and $g$,

$$ \begin{aligned} (f\Theta\star g\Theta)(t, s) &= \int_{-\infty}^{\infty} f(t, \tau) \Theta(t-\tau) g(\tau, s) \Theta(\tau-s) d\tau, \\ &= \int_{s}^{t} f(t, \tau) g(\tau, s) d\tau \Theta(t-s) = (f \star_v g)(t, s) \Theta. \end{aligned} $$

The Volterra composition, which naturally appears in the Picard iterations is mathematically ill-behaved in general, lacking a proper unit and inverses for example [15]. All of these issues are lifted by the $\star$-product which provides the correct distributional context for non-autonomous linear differential systems and beyond (especially fractional and non-linear ordinary differential systems).

## 2. Relation to the convolution

In general both the $\star$-product and the Volterra composition differ from convolutions because $f, g$ in Eq. (2) and A, B in Eq. (3) may not depend solely on the difference between their arguments. If we explicitly consider two elements $k, l \in \mathcal{D}$ that are invariant under time translations, $k(t, s) = k(t-s, 0) \equiv k(t-s)$ and $l(t, s) = l(t-s, 0) \equiv l(t-s)$; then they both effectively depend on the single variable $t-s$ and their $\star$-product reduces to a convolution

$$ \begin{aligned} (k \star l)(t, s) &= \int_{-\infty}^{\infty} k(t, \sigma)l(\sigma, s)d\sigma, \\ &\equiv \int_{-\infty}^{\infty} k(t - \sigma)l(\sigma - s)d\sigma = (k * l)(t, s) \end{aligned} $$

In the context of ODEs this happens if and only if the coefficient matrix A commutes with itself at all times. Indeed, this is tied to the behavior of the evolution operator U that solves the N-ODE $\dot{U} = AU$. In general, A does not commute with itself at all times if and only if U is not invariant under time translations: $U(t, s) \neq U(t-s, 0)$. This confirms the necessity of working with bivariate objects in the most general theory.

Note that in all cases, the usual one-variable evolution operator is related to the bivariate one by the semi-group property, as follows $U(t, s) = U(t)U^{-1}(s)$. In an abuse of notation we will denote the bivariate and one-variable evolution operators by the same letter U.

## 3. Solution of non-autonomous systems with the $\star$-formalism

Let $A(t)$ be a time-dependent matrix all of whose entries are smooth functions of time and consider the system
of linear Non-autonomous Ordinary Differential Equations (N-ODE)

$$ \dot{U}(t, s) = A(t)U(t, s), \quad U(s, s) = \mathrm{Id}, \qquad (\mathrm{A}1) $$

for $t \ge s$ two real variables. Here $U(t,s)$ designates the matrix solution of the system, termed its evolution operator. The $\star$-product turns such N-ODEs and their solutions into purely linear $\star$-algebra. To see this let $U(t,s) := U(t,s)\Theta(t-s) \in \mathcal{D}$ and introduce the Green's function $G := (\delta^\prime\mathrm{Id}) \star U = \mathrm{Id}_\star + \dot{U}\Theta$. Then the N-ODE Eq. (A1) becomes [4]

$$ G - \mathrm{Id}_\star = A\Theta \star G. \qquad (\mathrm{A}2) $$

The solution of Eq. (A2) is

$$ G = (\mathrm{Id}_\star - A\Theta)^{-1}, \qquad (\mathrm{A}3) $$

and $U = \Theta \star G$ (because $\Theta \star \delta^\prime = \mathrm{Id}_\star$). Here the existence of the $\star$-inverse in $G$ is guaranteed by standard results from the analysis of Picard iterations [4] since

$$ \begin{aligned} G &= \sum_{n \ge 0} (A\Theta)^{*n}, \\ &= \mathrm{Id}_\star + A\Theta + \int_s^t A(t)A(\tau)d\tau \Theta + \\ & \qquad \int_s^t \int_s^\tau A(t)A(\tau)A(\tau')d\tau' d\tau \Theta + \dots \end{aligned} $$

which converges provided $\|A\|$ is finite over a time interval of interest. In addition, the above series shows that $\star$-resolvents are nothing but time-ordered exponentials

$$ G = \mathcal{T}e^{A(t)}, \qquad (\mathrm{A4}) $$

while $U = \Theta \star G = \mathcal{T}e^{\int_s^t A(\tau)d\tau} \Theta$.

The advantage of the $\star$-formalism's outlook is that all methods and results from linear algebra are immediately applicable to the autonomous and non-autonomous differential settings in the distributional sense. This for example implies the existence of formal, exact and explicit $\star$-continued fraction representations of $G$, known as pathsums [2, 4]. This also implies that standard procedure from linear algebra exist in the differential setting, such as Lanczos-triadiagonalization [5–7] and undoubtedly more. In this work we exploit the $\star$-product formalism to produce novel frame changes for ODEs and N-ODEs.

## Appendix B: Standard frame change

As presented in Section II A, the evolution operator $U$ corresponding to matrix $A(t) = A_0(t) + A_1(t)$ can be expressed as

$$ U = U_1 \star (\mathrm{Id}_\star - A_0 \star G_1)^{-1}. \qquad (\mathrm{B1}) $$

Let us now evaluate this explicitly. Firstly, concerning $A_0 \star G_1$ we have

$$ \begin{aligned} (A_0 \star G_1)(t,s) &= \int_{-\infty}^{\infty} A_0(t)\Theta(t-\tau)G_1(\tau,s)d\tau, \\ &= A_0(t)\int_{-\infty}^{\infty} \Theta(t-\tau)G_1(\tau,s)d\tau, \\ &= A_0(t)(\Theta \star G_1)(t,s) = A_0(t)U_1(t,s), \end{aligned} \qquad (\mathrm{B2}) $$

and since $t \ge s$ this is equivalent to $A_0 \star G_1 = A_0 U_1$. Since $U_1(t,s) = U_1(t)U_1^{-1}(s)$ (and if $A$ is Hermitian, $U^{-1} = U^\dagger$), it also follows that

$$ (A_0 \star G_1)^{*n} = A_0(t)U_1(t) \left( \int_s^t U_1^{-1}(\tau)A_0(\tau)U_1(\tau)d\tau \right)^{*n-1} U_1^{-1}(s), $$

so that,

$$ (\mathrm{Id}_\star - A_0 \star G_1)^{-1}(t,s) = \mathrm{Id}_\star + A_0(t)U_1(t) \mathcal{T}e^{\int_s^t U_1^{-1}(\tau)A_0(\tau)U_1(\tau)d\tau} U_1^{-1}(s). $$

Then

$$ \begin{aligned} U(t,s) &= (U_1 \star (\mathrm{Id}_\star - A_0 \star G_1)^{-1})(t,s), \\ &= U_1(t,s) + U_1(t) \int_s^t \{ U_1(\tau')^{-1}A_0(\tau')U_1(\tau') \\ & \qquad \mathcal{T}e^{\int_s^{\tau'} U_1^{-1}(\tau)A_0(\tau)U_1(\tau)d\tau} U_1^{-1}(s) \} d\tau', \end{aligned} $$

which evaluates to

$$ U(t,s) = U_1(t) \mathcal{T}e^{\int_s^t U_1^{-1}(\tau)A_0(\tau)U_1(\tau)d\tau} U_1^{-1}(s). $$

## Appendix C: Biframe presentation for $\dot{U}$

For $j = 0, 1$, by definition $G_j := \delta^\prime \star (U_j) = \mathrm{Id}_\star + \dot{U}_j = (\mathrm{Id}_\star - A_j)^{-1}$ so we have $\dot{U}_j = A_j \star G_j$. Then, Eq. (9) yields

$$ \begin{aligned} G &= G_0 \star (\mathrm{Id}_\star - \dot{U}_1 \star \dot{U}_0)^{-1} \star G_1, \\ &= G_0 \star \sum_{n=0}^{\infty} (\dot{U}_1 \star \dot{U}_0)^{*n} \star G_1. \end{aligned} $$

The series on the second line is guaranteed to converge for the same reasons as Picard iterations: bounding $\|\dot{u}_1 \star \dot{u}_0\| \le C$, the norms of the $n$th $\star$-powers of $\dot{U}_1 \star \dot{U}_0$ are $O(C^n/n!)$.
Given that $\dot{U} = G - \mathrm{Id}_\star$, this result is equivalent to the following unconditionally convergent series representation for $\dot{U}$ that reveals the symmetric roles played by $\dot{U}_0$ and $\dot{U}_1$,

$$ \begin{aligned} \dot{U} &= \dot{U}_0 + \dot{U}_1 + \dot{U}_0 \star \dot{U}_1 + \dot{U}_1 \star \dot{U}_0 \\ & \quad + \dot{U}_0 \star \dot{U}_1 \star \dot{U}_0 + \dot{U}_1 \star \dot{U}_0 \star \dot{U}_1 \\ & \quad + \dot{U}_0 \star \dot{U}_1 \star \dot{U}_0 \star \dot{U}_1 + \dot{U}_1 \star \dot{U}_0 \star \dot{U}_1 \star \dot{U}_0 \\ & \quad + \dots \end{aligned} $$

and from there $U = \mathrm{Id} + \Theta \star \dot{U}$. The above series is just the sum over all alternating products of $\dot{U}_0$ with $\dot{U}_1$, manifestly invariant under exchange of indices $0 \leftrightarrow 1$.
## Appendix D: Proof of Eq. (11b)

To arrive at Eq. (11b) from Eq. (10) we need only evaluating $A_1 \star G_1 \star A_0 \star G_0$. By Eq. (B2) we have $A_j \star G_j = A_j(t)U_j(t)U_j^{-1}(s)$ for $j = 0, 1$. Then

$$ 
\begin{aligned}
A_1 \star G_1 \star A_0 \star G_0 &= \int_s^t A_1(t) U_1(t) U_1^{-1}(\sigma) A_0(\sigma) U_0(\sigma) U_0^{-1}(s) d\sigma \\
&= A_1(t) U_1(t) \int_s^t U_1^{-1}(\sigma) A_0(\sigma) U_0(\sigma) d\sigma U_0^{-1}(s).
\end{aligned}
$$

It follows that the $\star$-resolvent of $A_1\Theta \star G_1 \star A_0\Theta \star G_0$ is the time-ordered exponential of the above, which is Eq. (11b).

## Appendix E: Alternative form for the biframe

To obtain the biframe expression, instead of Eq., (8b) we can equally start from the ordinary identity between matrix resolvents

$$ R = R_0 (I - R_1 \cdot M_1 \cdot R_0 \cdot M_0)^{-1} \cdot R_1, $$

which implies the alternative formula for Green's functions

$$ G = G_0 \star (I_{\star} - G_1 \star A_1 \Theta \star G_0 \star A_0 \Theta)^{-1} \star G_1. \quad (E1) $$

Now let us consider the content of the $\star$-resolvent in more details. Using Eq. (B2) we get

$$ G_1 \star A_1 \Theta \star G_0 \star A_0 = G_1 \star (A_1 U_0 \Theta) \star A_0 \Theta $$

and so

$$ 
\begin{aligned}
(G_1 \star A_1 \star G_0 \star A_0)^{*2}, \\
&= (G_1 \star (A_1 U_0) \star A_0)^{*2}, \\
&= G_1 \star (A_1 U_0) \star (A_0 U_1) \star (A_1 U_0) \star A_0,
\end{aligned}
$$

since $(A_0 \star G_1)(t, s) = A_0(t)U_1(t, s)$, again by Eq. (B2). Continuing in this fashion we verify by induction that for $n \ge 1$,

$$ 
\begin{aligned}
(G_1 \star A_1 \star G_0 \star A_0)^{*n} \\
&= G_1 \star ((A_1 U_0) \star (A_0 U_1))^{*n-1} \star (A_1 U_0) \star A_0
\end{aligned}
$$

which yields, for $n \ge 1$,

$$ 
\begin{aligned}
(G_1 \star A_1 \star G_0 \star A_0)^{*n} \star G_1 \\
&= G_1 \star ((A_1 U_0) \star (A_0 U_1))^{*n-1} \star (A_1 U_0) \star (A_0 U_1) \\
&= G_1 \star ((A_1 U_0) \star (A_0 U_1))^{*n}.
\end{aligned}
$$

Finally, Eq. (E1) now reads

$$ 
\begin{aligned}
G &= G_0 \star G_1 \sum_{n=0}^{\infty} ((A_1 U_0) \star (A_0 U_1))^{*n}, \\
&= G_0 \star G_1 \star (I_{\star} - (A_1 U_0) \star (A_0 U_1))^{*-1},
\end{aligned}
$$

which immediately implies

$$ U = U_0 \star G_1 \star (I_{\star} - (A_1 U_0) \star (A_0 U_1))^{*-1}. \quad (E2) $$

Since $U_0(t, s) = U_0(t)U_0(s)$ and similarly for $U_1$, we have

$$ (A_1 U_0) \star (A_0 U_1) = A_1(t) U_0(t) \int_s^t U_0^{-1}(\tau) A_0(\tau) U_1(\tau) d\tau U_1^{-1}(s). $$

Putting everything together this leads to an alternative presentation of the biframe,

$$ U = U_0 \star G_1 \star \mathcal{T} e^{A_1(t)U_0(t) \int_s^t U_0^{-1}(\tau) A_0(\tau) U_1(\tau) d\tau} U_1^{-1}(s). \quad (E3) $$

Compare with Eq. (11b), in particular the content of the time-ordered exponential.

## Appendix F: Biframe change for constant Hamiltonians A

In the situation where the coefficient matrix A is time-independent or, more generally, commutes with itself at different times, much simplifications occur with respect to the general formula provided in the main text. First of all, in such situation the biframe operator is invariant under time-translations $\mathcal{B}(t-s, 0) = \mathcal{B}(t, s)$. This means that $\mathcal{B}$ effectively depends on a single variable (the difference $t-s$); and that all $\star$-products simplify to convolutions. Consequently the Laplace/Fourier domain may be employed, as for an autonomous ODE, and in this domain Eq. (11b) becomes

$$ \tilde{\mathcal{B}}(z) = \tilde{A}_1(z) (I - \tilde{A}_1(z))^{-1} \tilde{A}_0(z) (I - \tilde{A}_0(z))^{-1}, \quad (F1a) $$

while the Laplace transform of the time-ordered exponential of $\mathcal{B}$ is $(I - \tilde{\mathcal{B}}(z))^{-1}$ and so

$$ \tilde{U}(z) = z^{-1} (I - \tilde{A}_0(z))^{-1} (I - \tilde{\mathcal{B}}(z))^{-1} (I - \tilde{A}_1(z))^{-1}. \quad (F1b) $$

In these expressions $z$ is the Laplace/Fourier domain variable and $\tilde{A}_i$, $\tilde{\mathcal{B}}$ are the Laplace/Fourier transforms of $A_i$ and $\mathcal{B}$, respectively. Similarly $\tilde{G}_1(z) = (I - \tilde{A}_1(z))^{-1}$ and $\tilde{U}_0(z) = z^{-1}(I - \tilde{A}_0(z))^{-1}$.

Although Eqs. (F1) are only analytically correct in the Laplace/Fourier domain when A self-commutes at different times, they are generally valid for numerical calculations in the time domain. Indeed, once time is discretized, $\star$-products become ordinary matrix products, and $\star$-resolvents turn into ordinary matrix resolvents; see [1, 3, 9–11]. Then Eqs. (11) directly take the form of Eqs. (F1) with no Laplace/Fourier transformations involved, no matter how $A(t)$ depends on $t$.
[1] T. Birkandan, P.-L. Giscard, and A. Tamar. Computations of general Heun functions from their integral series representations. In *2021 Days on Diffraction (DD)*, pages 12–18, 2021.

[2] P.-L. Giscard and C. Bonhomme. Dynamics of quantum systems driven by time-varying Hamiltonians: Solution for the Bloch-Siegert Hamiltonian and applications to NMR. *Physical Review Research*, 2:023081, 2020.

[3] P.-L. Giscard and M. Foroozandeh. Exact solutions for the time-evolution of quantum spin systems under arbitrary waveforms using algebraic graph theory. *Computer Physics Communications*, 282:108561, 2023.

[4] P.-L. Giscard, K. Lui, S. J. Thwaite, and D. Jaksch. An exact formulation of the time-ordered exponential using pathsums. *Journal of Mathematical Physics*, 56(5):053503, 2015.

[5] P.-L. Giscard and S. Pozza. Lanczos-like algorithm for the time-ordered exponential: The $\star$-inverse problem. *Applications of Mathematics*, 65:807–827, 2020.

[6] P.-L. Giscard and S. Pozza. Tridiagonalization of systems of coupled linear differential equations with variable coefficients by a Lanczos-like method. *Linear Algebra and its Applications*, 624:153–173, 2021.

[7] P.-L. Giscard and S. Pozza. A Lanczos-like method for non-autonomous linear ordinary differential equations. *Bolletino dell'Unione Matematica Italiana*, 16:81–102, 2023.

[8] T. Mikami, S. Kitamura, K. Yasuda, N. Tsuji, T. Oka, and H. Aoki. Brillouin-Wigner theory for high-frequency expansion in periodically driven systems: Application to Floquet

topological insulators. *Phys. Rev. B*, 93:144307, Apr 2016.

[9] S. Pozza. A new closed-form expression for the solution of odes in a ring of distributions and its connection with the matrix algebra. *Linear and Multilinear Algebra*, 0(0):1–11, 2024.

[10] S. Pozza and N. Van Buggenhout. A new matrix equation expression for the solution of non-autonomous linear systems of odes. *PAMM*, 22(1):e202200117, 2023.

[11] S. Pozza and N. Van Buggenhout. The $\star$-product approach for linear ODEs: A numerical study of the scalar case. In *Programs and Algorithms of Numerical Mathematics*, pages 187–198. Institute of Mathematics CAS, 2023.

[12] M. Ryckebusch, A. Bouhamidi, and P.-L. Giscard. A Fréchet Lie group on distributions. *J. Math. Anal. Appl.*, 546(1):129195, 2025.

[13] J. Venkatraman, X. Xiao, R. G. Cortiñas, A. Eickbusch, and M. H. Devoret. Static Effective Hamiltonian of a Rapidly Driven Nonlinear System. *Phys. Rev. Lett.*, 129:100601, Aug 2022.

[14] V. Volterra. Sopra una proprietà generale delle equazioni integrali ed integrodifferenziali. *Atti della Reale Accademia dei Lincei, Rendiconti - Classe di Scienze Fisiche Matematiche e Naturali*, 20:79–88, 1911.

[15] V. Volterra and J. Pérès. *Leçons sur la composition et les fonctions permutables*. Gauthier-Villars, Paris, 1924. Recent edition by J. Gabay Ed., 2008, ISBN 9782876473232.