We study the propagation of wavepackets along curved interfaces between topological, magnetic materials. Our Hamiltonian is a massive Dirac operator with a magnetic potential. We construct semiclassical wavepackets propagating along the curved interface as adiabatic modulations of straight edge states under constant magnetic fields. While in the magnetic-free case, the wavepackets propagate coherently at speed one, here they experience slowdown, dispersion, and Aharonov–Bohm effects. Several numerical simulations illustrate our results.

§ INTRODUCTION

This paper analyzes wavepackets propagating along an interface between two topologically distinct materials, in the presence of an external magnetic field. It extends constructions carried out in <cit.> for magnetic-free models. We represent here the electron dynamics via a two-dimensional Dirac equation:

\begin{equation}\label{eq:D1}
(\eps D_t + \slashed D)\Psi(t,x) = 0, \quad (t,x) \in \mathbb R \times \mathbb R^2.
\end{equation}

In (<ref>),  $\slashed D$ denotes a Dirac operator with sign-changing mass and magnetic field:

\begin{equation}\label{eq:Dirac}
\slashed D = (\eps D_1-A_1(x))\sigma_1+(\eps D_2-A_2(x))\sigma_2 + \kappa(x)\sigma_3, \qquad \text{where:}
\end{equation}

* $\epsi > 0$ is a small semiclassical parameter and $\epsi D_t=-i \epsi \partial_t$, $\epsi D_j=-i \epsi\partial_j$ denote the self-adjoint semiclassical derivatives;

* $A=(A_1,A_2)^t \in C^\infty(\R^2,\R^2)$ is a magnetic potential with $\nabla A \in C^\infty_b$ (i.e. smooth with all derivatives uniformly bounded), inducing the magnetic field $B=\partial_1A_2-\partial_2A_1$;

* $\kappa \in C^\infty(\R^2,\R)$ has varying sign and satisfies $\nabla\kappa\in C^\infty_b$;

* $\sigma_1, \sigma_2, \sigma_3$ are the standard $2\times2$ Pauli matrices.

The sign of the domain wall $\kappa$ characterizes the topological phase of the material. Specifically, under the transversality condition

\begin{equation}\label{eq-3a}
\kappa(y) = 0 \qquad \Rightarrow \qquad \nabla \kappa(y) \neq 0,
\end{equation}

the interface $\Gamma = \kappa^{-1}(0)$ separates regions of distinct local topology B19b,bal2021edge. For analytic reasons, we consider here a uniform version of (<ref>):

\begin{equation}
\inf\{ \big|\nabla \kappa(x)\big| : x \in \Gamma\} >0.
\end{equation}

The bulk-edge correspondence for Dirac operators <cit.> predicts that a tubular neighborhood of $\Gamma = \kappa^{-1}(0)$ supports asymmetric currents, hence (some analogue of) edge states for small $\epsi$. For vanishing magnetic potentials, we constructed in <cit.> long-lived solutions to (<ref>). These were confined and propagating at speed one along $\Gamma$. We referred to them as dynamical edge states, since we derived them as time-dependent adiabatic modulations of straight edge states. This paper extends the construction to the magnetic case. A non-zero magnetic field induces several new phenomena for dynamical edge states:

(i) A systematic slowdown, see Figure <ref>;

(ii) A large phase-shift, generating a Aharonov–Bohm effect when $\Gamma$ is a loop;

(iii) In general, a mesoscopic dispersion along $\Gamma$.

We set a few notations in <ref>, state a simplified main result in <ref>, and detail the effects (i)–(iii) in <ref>.

Snapshots of the evolution of Gaussian wavepackets propagating along a straight interface $\kappa(x)=x_2$, under a constant magnetic field $B=0,0.5,1,1.5,2$ (from right to left), computed numerically. The packets slow down in stronger fields.

at (0,0) ;
[blue!70] at (.7,-2.5) $B=0$;
[blue!70] at (-2.5,-1.8) $B=2$;

§.§ Notations
We first define the normal and tangent vector fields to the level sets of $\kappa$:

\begin{equation}
\label{eq:tang_vec}
n(x) = \frac{\nabla \kappa(x)}{|\nabla \kappa(x)|}, \qquad \tau(x) = J n(x), \qquad J=\matrice{ 0&-1\\ 1 & 0}.
\end{equation}

Given $y_0 \in \Gamma$, we let $y_t$ be the solution of the ODE

\begin{equation}\label{eq:yt}
\dot y_t = c(y_t) \tau(y_t), \qquad c(y) := \dfrac{|\nabla \kappa(y)|}{\sqrt{|\nabla \kappa(y)|^2 + B(y)^2}},
\end{equation}

where $B = \nabla \times A$ is the magnetic field. Note that $\p_t \kappa(y_t) = 0$ since $\nabla \kappa(x) \cdot \tau(x) = 0$, hence $y_t \in \Gamma$ for any $t \in \R$. We then define the quantities

\begin{equation}
B_t=B(y_t), \quad n_t = \frac{\nabla \kappa(y_t)}{|\nabla \kappa(y_t)|}, \quad r_t = |\nabla \kappa(y_t)|, \qquad \rho_t = \sqrt{r_t^2+B_t^2}, \quad \gamma_t = \frac{B_t}{\rho_t^2}.
\end{equation}

To capture the local geometry of the interface near $y_t$, we define two smoothly varying angles $\varphi_t$ and $\te_t$ (and the corresponding clockwise rotation $R_{\te_t}$) by

\begin{equation}\label{eq:coefst}
\cos \varphi_t= \frac{r_t}{\rho_t} := c_t, \quad \sin \varphi_t= \frac{B_t}{\rho_t}:= s_t, %\ \ \text{ and } \quad R_{\theta_t} n_t=e_2,
\quad R_{\theta_t} = \matrice{ \cos\theta_t & \sin\theta_t \\ -\sin\theta_t & \cos\theta_t },
\end{equation}

see Figure <ref>; we assume $\varphi_0$ and $\theta_0$ belong to $[0,2\pi)$ for concreteness.
With these notations, $R_{\theta_t} n_t=e_2$ and $\dot{y_t} = c_t \tau_t$ with $\tau_t=Jn_t$.
Differentiating the equation defining $\theta_t$, we observe that $\dot n_t+\dot\theta_t\tau_t=0$. Since $\dot y_t=c_t\tau_t$, we deduce that $\dot\theta_t=-c_t (\tau\cdot\partial_\tau n)(y_t)$. In particular, $|\dot\theta_t|=c_t K_t$ where $K_t$ is the curvature of the curve $\Gamma$ at $y_t$.

(-5,-3) rectangle (5,2);

[red!10, domain=-5:5, variable=]
(-5, 3)
– plot (, )
– (5, 3)
– cycle;

[blue!10, domain=-5:5, variable=]
(-5, -4)
– plot (, )
– (5, -4)
– cycle;

[domain=-5:5 , smooth, variable=, ultra thick] plot (, );

[domain=-1:-2.5 , smooth, variable=, ultra thick,red,->] plot (, .3+/2.5+(1+*/6.25)^(1/2)-1.5 +  tanh(-3) +1 );

[red] at (-1.8, -.3) current;

[dotted, thin] (-5,3) grid (5,-4);
[thick,->] (-3,-4) – (-3,3);
[thick,->] (-5,-2) – (5,-2);
[red] at (-4, 1.5) $\kappa > 0$;
[dblue] at (3, -2.5) $\kappa < 0$;
at (-4, -.8) $\kappa^{-1}(0)$;

[domain=0:2 , smooth, variable=, ultra thick, red, dashed] plot (, 1.21+0.601*(-2));

[red, ultra thick, ->] (0,0) – (-1,-0.601);
[red, ultra thick, ->] (0,0) – (-0.601,1);
[thick] (0,0) – (2,0);
[red, domain=0:30.9 , smooth, variable=, ultra thick, ->] plot (cos(), sin());

[red] at (1.3, .4) $\theta_t$;

[red] at (1.8, -.3) $y_t$;
[red,fill=red] (1.8, 0) circle (1mm);
[red] at (1.7, 1.3) $n_t$;
[red] at (1.2, -1) $\tau_t$;

Geometry of curved interface.

We also introduce the unitary pullback operator $\RR_{\te_t}$ by $R_{\te_t}$ as $\RR_{\te_t} g(z) = g(R_{\te_t} z)$. Associated to $\te_t$ and $\varphi_t$ are two spinorial rotations,

\begin{equation}
U_{2,\varphi_t} = e^{-i\frac{\varphi_t}2\sigma_2}, \qquad U_{3,\theta_t} = e^{-i\frac{\theta_t}2\sigma_3}.
\end{equation}

§.§ Simplified main result
For the sake of simplicity, we consider a restricted setup: the domain wall satisfies a geometric condition and the magnetic field is constant. Our main result, <ref>, takes then a clearer form.

The aforementioned assumptions are:

* $\kappa \in C^\infty(\R^2,\R)$ with $\nabla \kappa \in C^\infty_b$

\begin{equation}\label{eq:00g}
y \in \kappa^{-1}(0) \quad \Rightarrow \quad \big| \nabla \kappa(y) \big| = 1, \quad \Delta \kappa(y) = 0;
\end{equation}

* The magnetic field $B$ is constant.

While (<ref>) is analytically restrictive, it is not geometrically restrictive: any one-dimensional submanifold of $\R^2$ is the zero set of a function $\kappa$ satisfying (<ref>).
For instance, for the straight edge $\Gamma = \R e_2$, we can choose $\kappa(y)=y_2$ while for the unit circle $\Gamma = \mathbb{S}^1$ we can choose $\kappa(y)=\ln |y|$ (near $|y|=1$). When the  two above conditions hold, the quantities $r_t, B_t, \rho_t, \gamma_t, \varphi_t, c_t$ and $s_t$ do not depend on $t$, and we omit the subscript $t$.

The most elementary setup with these two conditions consists of $\kappa(x) = x_2$ and $A = -B x_2 e_1$. The corresponding Dirac operator is:

\begin{equation}\label{eq-2a}
\slashed{D}_{0,B} := ( \epsi D_1+Bx_2) \sigma_1 + \epsi D_2 \sigma_2 + x_2 \sigma_3.
\end{equation}

We remark that the magnetic potential vanishes along the interface $\R e_1$ and is parallel to it. The equation $(\epsi D_t + \slashed{D}_{0,B}) \Psi = 0$ admits an explicit family of non-dispersive wavepacket solutions that propagate along $\Gamma$ at speed $c = (1+B^2)^{-1/2}$:

\begin{equation}\label{eq:00i}
\Psi_{0,B}(x) := \frac{1}{\sqrt\eps} \psi_{0,B}\left( \dfrac{x-ct e_1}{\sqrt{\epsi}} \right), \quad \psi_{0,B}(z) := \int_\R e^{i\xione z_1-\frac{1}{2}\rho ( z_2+\gamma \xione)^2} \fhat(\xione) d\xione \cdot U_{2,\varphi}\matrice{1  \\ -1},
\end{equation}

where $\fhat$ is any Schwartz-class function, $\rho = \sqrt{1+B^2}$, $\gamma = \frac{B}{1+B^2}$, and $\varphi=\arctan B$.

We extend this statement to interfaces tilted by an angle $\te \in \R$. Let $\slashed{D}_{\te,B}$ be the Dirac operator with domain wall and magnetic potential

\begin{equation}
\kappa(x) = n \cdot x, \quad A(x) = B \kappa(x) \tau, \quad n = \matrice{-\sin \te \\ \cos \te}, \quad \tau = -\matrice{\cos \te \\ \sin \te}.
\end{equation}

The equation $(\epsi D_t + \slashed{D}_{\te,B}) \Psi = 0$ is unitarily equivalent to the case $\te=0$ and admits a family of solutions constructed from $\Psi_{0,B}$:

\begin{equation}
\Psi_{\te,B}(x) := \frac{1}{\sqrt\eps} \psi_{\te,B}\left( \dfrac{x-ct \tau}{\sqrt{\epsi}} \right), \qquad  \psi_{\te,B}(z) := U_{3,\te} \RR_{\te} \psi_{0,B}(z).
\end{equation}

Our main result produces approximate solutions to $(\epsi D_t + \slashed{D}) \Psi = 0$ as modulations of $\psi_{\te,B}$. To state it, we need the distribution $g_t$ on $\R^2$ defined by:

\begin{equation}\label{eq:00j}
g_t(z) := \int_{\R^2} e^{i\zeta z + i \gamma (\te_t-\te_0) \zeta_1^2} d\zeta = \dfrac{e^{-i\pi/4}}{|2\pi\gamma (\te_t-\te_0)|^{1/2}} e^{i \frac{-z_1^2}{4 \gamma (\te_t-\te_0)}} \cdot \delta_0(z_2).
\end{equation}

The second equality is valid for $\te_t \neq \te_0$ and should be replaced by a Dirac mass when $\te_t = \te_0$. Our simplified theorem for constant magnetic fields reads then as follows.

Assume that $\kappa$ satisfies the condition (<ref>) and that $B$ is constant. For any $\rT >0$, the equation $(\epsi D_t + \slashed{D} ) \Psi = 0$
admits solutions that satisfy, uniformly for $t \in [-\rT,\rT]$ and $x \in \R^2$:

\begin{equation}\label{eq:00ii}
\Psi(t,x) = \dfrac{e^{\frac{i}{\epsi}\chi(t,x)}}{\sqrt{\epsi}}  \psi\left( t,\dfrac{x-y_t}{\sqrt{\epsi}} \right) + O_{L^2}(\epsi^{1/2}), \qquad \text{where} \quad  \psi(t,z) := \big(\RR_{\te_t} g_t\big) * \psi_{\te_t,B}(z),
\end{equation}

\begin{equation}
\text{and} \quad
\chi(t,x) = \int_0^t \dot{y_s} \cdot A(y_s) ds + A(y_t) \cdot (x-y_t) + (x-y_t) \cdot \left( \nabla A(y_t)^\top - B n_t \tau_t^\top\right) (x-y_t).
\end{equation}

The leading order term in (<ref>):

(i) Propagates at speed $c = (1+B^2)^{-1/2}$ along $\Gamma$, in the prescribed direction $\tau$.

(ii) Is semiclassically localized at the phase-space point $(y_t,A_t)$ with $A_t=A(y_t)$. This is where the eigenvalues of the symbol of $\Di$ are degenerate.
While the rapid oscillations generated by $A_t$ can be locally gauged away, they cannot be globally neglected when $\Gamma$ is a closed loop: this is the Aharonov–Bohm effect.

(iii) Disperses along $\Gamma$ at rate prescribed by the difference $\te_t-\te_0$:

\begin{equation}\label{eq:00k}
\big| \psi(t,z) \big| \leq \dfrac{C}{1+|\te_t-\te_0|^{1/2}}\sup_{z \in \R^2} \big| \psi(0,z) \big|.
\end{equation}

In particular, the leading part in (<ref>) is controlled by $\epsi^{-1/2} |\te_t-\te_0|^{-1/2}$.
This is a relatively weak dispersion as it comes from dispersion of the wave envelop  $\R_{\te_t g_t} \star \psi_{\te_t,B}$ rather than dispersive relations of plane waves. It produces effects for times of order one, in contrast with dispersion in e.g. the semiclassical Schrödinger equation which arises at time $\epsi$.

We detail these three effects in <ref>. Theorem <ref> will extend the result of Theorem <ref> to cover varying magnetic fields, general domain walls and longer times of validity. The corresponding wavepackets will more generally have a variable speed, given by the ODE (<ref>); a phase with properties identical to (ii) above; and a more complicated rate of dispersion.

This is a write-up for me. We first construct $\psi_{\te,B}$, we first observe that $\beta(x) = B$, $\tilde{A} = A$, hence $\chi = 0$. We have the formula

\begin{equation}
R_\te = \qquad,  U_{3,\te} = , \qquad U_{2,\varphi} = \qquad
\end{equation}

for the spatial and spinorial rotations. The operators $\mathcal{V}$ and $\mathcal{V}^*$ are

\begin{equation}
\end{equation}

After reduction, our approach produces the solution

\begin{equation}
f(\xi_1) e^{-\frac{\rho \zeta_2^2}{2}} \matrice{1  \\ -1}.
\end{equation}

Acting with $\mathcal{V}$ produces

\begin{equation}
\psi_{0, B} = \int_\R e^{i\xi_1z_1} f(\xi_1) e^{-\frac{\rho ( z_2+\gamma \xi_1)^2}{2}} d\xi_1 \cdot U_{2,\varphi}\matrice{1  \\ -1}.
\end{equation}

If this is radial, then $\psi_{\te,B}$ comes with essentially no change. However I am not sure that this is useful, as we cannot expect that the solution given by the theorem will be radially symmetric (it disperses along the edge only, so its level sets balloon in the direction parallel to the edge only: they cannot be radial).

Now we look for $\psi$. We find $\psi = \mathfrak{U}_t a$ with

\begin{equation}
a(t,\xi_1,\zeta_2) = e^{i\nu_t \frac{\xi_1^2}{2}} f(\xi_1) e^{-\frac{\rho \zeta_2^2}{2}} \matrice{1  \\ -1}, \qquad \nu_t = \gamma(\te_0-\te_t).
\end{equation}

Now act with $\mathcal{V}$ (which does not depend on $t$). We get

\begin{equation}
\mathcal{V}a(t,z) = \mathcal{V} e^{i\nu_t \frac{\xi_1^2}{2}} f(\xi_1) e^{-\frac{\rho \zeta_2^2}{2}} \matrice{1  \\ -1} = \mathcal{V} e^{i\nu_t \frac{\xi_1^2}{2}} \mathcal{V}^* \psi_{0,B}.
\end{equation}

We observe that

\begin{align}
\mathcal{V} e^{i\nu_t \frac{\xi_1^2}{2}} \mathcal{V}^* f(z)
% &  = \dfrac{1}{\sqrt{2\pi}} \int_\R e^{i\xi_1z_1 + i\nu_t \frac{\xi_1^2}{2}} (\mathcal{V}^* f)(\xi_1,z_2+\gamma \xi_1)      \\     & = \dfrac{1}{2\pi} \int_\R e^{i\xi_1 z_1 + i\nu_t \frac{\xi_1^2}{2}} \int_\R e^{-iy_1\xi_1} f(y_1,z_2) dy_1 \\     &=  \dfrac{1}{2\pi}\int_{\R^2}  e^{i\xi_1(z_1 - y_1) + i\nu_t \frac{\xi_1^2}{2}} f(y_1,z_2) d\xi_1   dy_1
= \dfrac{1}{|\nu_t|^{1/2}} e^{i \frac{z_1^2}{2\nu_t}} *_1 f(z),
\end{align}

where $*{}_1$ denotes convolution with respect to the first variable.

§.§ Effects of the magnetic field
We comment here on the structure of the wavepacket (<ref>), when $B$ is constant and $\kappa$ satisfies (<ref>); and detail how things change when these conditions are relaxed (see Theorem <ref>).

The left figure illustrates the Aharonov-Bohm effect showing the phase of one spinor component as the wavepacket performs one full revolution around the unit circle with $A(x)=\vert x \vert^{-1}$. The figure in the center shows a sequence of snapshots of the state undergoing one full revolution on a circular edge profile $\kappa(y) :=\log(\vert y \vert)$ with constant magnetic field $B=1/2$ and on the right for $B=1$ corresponding to a maximal dispersion rate $\gamma=1/2$.

Aharonov–Bohm effect for the domain wall $\kappa(y) = \ln |y|$ (corresponding to a circle interface). The plot shows the rapid evolution of the phase of the first spinor component as the wavepacket performs a single revolution around the circle. See <ref> for details.

at (0,0) ;

Snapshots of the evolution of a wavepacket around the unit disc

(with $\kappa(y) = \ln |y|$ and $\dot\theta=c$) in magnetic fields $B = 1/2,1/\sqrt{2},1$ (top row) and $B=2,3,4$ (bottom row).
We observe gradual dispersion consistent with our predicted dispersion rate $\gamma |\te_t-\te_0| = \gamma c t = B (1+B^2)^{-3/2} t$, see (<ref>). For equal distance of propagation the dispersion is maximal when $B=1$ while for equal time of propagation the dispersion is strongest when $B=1/\sqrt{2}$.

Turning on a magnetic field systematically slows down the propagation. In the general setup of Theorem <ref>, the wavepackets move at speed $(1+B_t^2)^{-1/2}$, which is smaller than $1$ whenever the magnetic field does not vanish.

The wavepacket is semiclassically localized at $(y_t,A_t)$. This point lies in the crossing set of the semiclassical symbol

\begin{equation}
\slashed{D}(x,\xi) = \big(\xi_1-A_1(x)\big) \sigma_1 + \big(\xi_2-A_2(x)\big) \sigma_2  + \kappa(x) \sigma_3;
\end{equation}

that is, the eigenvalues of $\slashed{D}(y_t,A_t)$ are repeated. Hence,  $(y_t,A_t)$ is an exotic semiclassical trajectory, in the sense that it is not among those predicted by standard propagation of singularity, such as <cit.>.

The semiclassical action $\dot{y_t} \cdot A_t$ generates a large phase-shift $e^{\frac{i}{\epsi}\int_0^t \dot{y_s} \cdot A_s ds}$ that can be locally – but not globally – gauged away. When $\Gamma$ is a loop, after a full revolution the phase shift relates to the magnetic flux $\Phi = \int_\Gamma A$: it is $e^{\frac{i}{\epsi} \Phi}$, see Figure <ref> for the case of the circle.
This is the Aharonov–Bohm effect.

The envelop of the wavepacket typically disperses along $\Gamma$. Under the geometric condition (<ref>), the rate of change of dispersion is $|\te_t-\te_0|^{-1/2}$. Indeed, writing $\psi$ in terms of $g_t$ and $\psi_{0,B}$, and using the formula (<ref>), we have for $|\te_t-\te_0| \geq 1$:

\begin{align}
\sup_{z \in \R^2} \big| \psi(t,z) \big| = \sup_{z \in \R^2} \big|g_t * \psi_{0,B}(z)\big| & \leq \dfrac{1}{\sqrt{ 2\pi \gamma |\te_t-\te_0|}} \sup_{z_2 \in \R} \int_\R |\psi(y_1,z_2)| dy_1
\\
& \leq \dfrac{C}{\sqrt{\gamma |\te_t-\te_0|}} = \dfrac{C'}{\sqrt{ \gamma |\te_t-\te_0|}} \sup_{z \in \R^2} \big| \psi(0,z) \big|. \label{eq-7m}
\end{align}

This yields (<ref>). While the rate of change of dispersion $|\te_t-\te_0|^{-1/2}$ is bounded above if $\Gamma$ is asymptotically flat, it can be as small as $t^{-1/2}$ when $\Gamma$ is a loop or a spiral. In these cases, the wavepacket loses coherence over long times as displayed in Figure <ref>; it should be noted however that larger magnetic fields do not necessarily give rise to larger dispersion. In the general setup of Theorem <ref>, the rate of dispersion $\nu_t$ takes a more complicated form: see Lemma <ref>, (<ref>) and (<ref>). We comment that dispersion can have long-time effects. When the dispersion is strongest,  $\nu_t \sim t$ and our construction holds up to times $\rT\ll\eps^{-1/8}$; when it is weakest, $\nu_t = O(1)$ and we recover the time of validity  $\rT\ll\eps^{-1/2}$ of <cit.>.

When $B$ or $|\nabla \kappa|$ vary along $\Gamma$ (that is, outside the setup of Theorem <ref>), an additional effect emerges: time-dependent anisotropic compression/stretching in the normal and tangent directions of $\Gamma$. We refer the reader to the formula (<ref>) and Theorem <ref>.
The compression factors, $\sqrt{\rho_t}$ (in the normal direction) and $c_0^2/c_t^2$ (in the tangent direction), remain bounded above and below in the limit $t \rightarrow \infty$. This contrasts with the parameter $\nu_t$ that controls the dispersion, which can grow like $t$.

§.§ Strategy of proof
Conceptually speaking, our approach consists in constructing successive transformations of $\Di$ that bring us closer to the flat Dirac operator $\Di_{0,B}$ of (<ref>). The main steps are as follows:

1. In <ref>, we conjugate $\Di$ by a gauge transform $e^{i\chi(t,x)/\epsi}$, cooked up so that the resulting magnetic potential $A-\nabla \chi$ vanishes along $\Gamma$ and is tangent to the level sets of $\kappa$. These are shared features with the flat Dirac operator $\Di_{0,B}$.  The expansion of the gauge term $e^{i\chi(t,x)/\epsi}$ near $y_t$ produces the large oscillatory phase of (<ref>).

2. In <ref> we look for solutions to the gauge-modified Dirac equation $(\epsi D_t + \tilde{\slashed{D}}) \tilde{\Psi} = 0$ as semiclassical wavepackets localized at $(y_t,0)$. A formal Taylor expansion produces a hierarchy of equations for the envelops.

3. In <ref>, we perform a series of spatial and spinorial rotations on the leading equation. The first two rotations (which already appear in <cit.>) flatten the interface; the last is magnetically induced and is among the new ingredients. The result is the leading equation that one would get starting from $\slashed{D}_{0,B}$.

4. Up to a partial shifted Fourier transform, the leading operator takes the same form as in the absence of magnetic fields <cit.>. In <ref>, we compute explicitly its kernel and prove stability estimates – later needed for the subleading transport equation.

5. In <ref> we explicitly integrate the transport equation. In Fourier variables, the solutions have a quadratic phase with Hessian $\te_t-\te_0$ (in the setup of <ref>). In physical space, this transfers to dispersion at rate $|\te_t-\te_0|^{-1/2}$.

6. Higher-order approximations are constructed iteratively in <ref>.  Combining them with the unitarity of $e^{-it\Di}$ and a Duhamel argument, we obtain our main result, Theorem <ref>. It gives an approximate solution to $(\epsi D_t + \Di) \Psi = 0$ that propagates along $\Gamma$, slowed down by the magnetic field, and explicitly expressed through the above transformations.

We then illustrate our findings with a series of numerical simulations in <ref>.

§.§ Related literature
For systems of semiclassical PDEs, the symbols that govern the macroscopic transport are the eigenvalues of the (matrix-valued) symbol. A symbolic diagonalization argument shows that their Hamiltonian flow governs the leading-order dynamics.
When the eigenvalues of the symbol are degenerate, the classical equations of motion break down. The situation studied here is among the simplest such cases. Our analysis show that the phase-space crossing set $\{ \big(x, A(x)\big) , x\in \Gamma\}$ support wavepackets with the dynamics (<ref>). As mentioned above, (<ref>) is an exotic semiclassical trajectory not predicted by the standard results on propagation of singularities <cit.>.

In other setups, wavepackets may e.g. start away from the crossing set, reach it, and undergo a Landau–Zener transition; see e.g. Hag94,HJ98,FG2,Col04.

From a physical point of a view, our motivation stems from the ubiquity of $\Di$ in the field of topological phases of matter WI,moessner2021topological, and in particular one-particle models of topological insulators and topological superconductors VO,Be13,PSB16,B19b,B19a, which generically come with conical points <cit.>. The domain wall $\kappa(x)$ models the interface between two topologically distinct insulating phases FLW16,Drouot:19,DW20. This in turn generates an asymmetric transport along the interface $\Gamma$ by a principle called the bulk-interface correspondence; see e.g. EG02,GP,PSB16,Drouot:19b,B20,bal2021topological,Drouot2020microlocal,bal2021. The wavepackets analyzed here encode this asymmetry; see <cit.> for a discussion when $A=0$.  The operator $\Di$ also emerges in the effective analysis of graphene  and its pseudomagnetic (strained) analogues, see for instance <cit.>.

Note that the magnetic field is essential in the integer quantum Hall effect TKNN,avron1994,BES94, which was the first observed example of topologically non-trivial state of matter.
There, the insulating gaps are obtained from the degenerate Landau levels associated to a constant magnetic field. This contrasts with the situation considered here: the non-trivial topology imposed by the domain wall is stable against magnetic contributions; see <cit.>. This means that the magnetic field does not influence the existence of edge states. It however affect their quantitative features, see <ref>.

Our results in Theorems <ref> and <ref> concern weakly dispersive wavepackets with a macroscopic center $y_t$ propagating along $\Gamma = \kappa^{-1}(0)$. Changing the metric (i.e. replacing $D$ by $\big(D_j+\Gamma_j(x))g^{jk}(x)\big)\sigma_k$) would likely preserve this structure. However, adding an electric potential $V(x)$ to $\Di$ modifies the energy landscape, hence the interface $\Gamma = \kappa^{-1}(0)$ between the two topological media cannot properly support wavepackets. These are likely to propagate instead within a thicker strip close to $\Gamma$, with splitting according to Landau–Zener transition rules.

§.§ Acknowledgments
This work is a sequel to <cit.>, which started during a 2020 AIM workshop, Mathematics of topological insulators. The authors thank the organizers: Daniel Freed, Gian Michele Graf, Rafe Mazzeo and Michael Weinstein. In addition, the authors are very grateful to Clotilde Fermanian Kammerer, Jianfeng Lu, and Alexander Watson for initial discussions. The authors acknowledge support form the NSF grants DMS-2118608 and DMS-2054589 (AD), DMS-1908736 and EFMA-1641100 (GB).

§ GAUGE TRANSFORMATION

This section constructs a gauge function that reduces $\Di$ to an operator with magnetic properties closer to those of the model operator (<ref>).

§.§ Equivalent tangent magnetic potential
We construct first a magnetic potential $\tilde A$ on $\R^2$ such that $\nabla \times \tilde{A} = B$ near $\Gamma$, and $\tilde A$ vanishes along $\Gamma=\kappa^{-1}(0)$ and is carried by the vector field $\tau$. In particular, the operator

\begin{equation}\label{eq-3f}
\tDi = \big( \epsi D_1 - \tilde A_1(x) \big) \sigma_1 + \big( \epsi D_2 - \tilde A_2(x) \big) \sigma_2 + \kappa(x) \sigma_3
\end{equation}

will share many of the characteristics of the model (<ref>): the magnetic potential is tangent to $\Gamma$ and vanishes along $\Gamma$. Moreover, from magnetic equivalence between $A$ and $\tilde A$, $\Di$ and $\tDi$ are locally conjugated. We use below the notation (with Euclidean distance $d$)

\begin{equation}
U_\eta = \big\{ x \in \R^2 : d(x,\Gamma) < \eta \big\}.
\end{equation}

There exist $\eta > 0$ and a function $\beta \in C^\infty(\R^2,\R)$ vanishing on $\Rm^2\backslash U_{2\eta}$ such that, defining $\tilde A:=\beta\kappa\tau$, we have

\begin{equation}\label{eq-0h}
\nabla \times (A - \tilde A) = \nabla \times (A - \beta\kappa\tau) = 0 \quad \text{ on } U_\eta; \qquad \beta = \dfrac{B}{|\nabla \kappa|} \quad \text{ on } \Gamma.
\end{equation}

We first note that it suffices to construct $\beta$ on $U_\eta$ and extend it to $\R^2$ as a smooth function with support in $U_{2\eta}$. For $f$ a smooth function, we observe that

\begin{equation}\label{eq-0w}
\nabla \times (f \tau) = \nabla f \times \tau +f \nabla \times \tau = \p_n f + f \nabla \times \tau,
\end{equation}

where we used $\nabla f \times \tau = \p_1 f \tau_2 - \p_2 f \tau_1 = \p_n f$. Let $B = \nabla \times A$. To find $\beta$ such that $B = \nabla \times (\beta \kappa \tau)$, we first solve

\begin{equation}\label{eq-0y}
\p_n f + f  \nabla \times \tau = B, \ \ \ \ f|_\Gamma = 0.
\end{equation}

We note that the coefficients of $\p_n$ are in $C^\infty_b$. In particular, the flow $e^{s\p_n}(x)$ is defined for all times. We now define the map $\Phi : \Gamma \times \R \rightarrow \R^2$ by

\begin{equation}
\Phi(x,s) = e^{s\p_n}(x).
\end{equation}

If $f$ solves (<ref>), then $\tf = f \circ \Phi$ solves

\begin{equation}
\p_s \tf + \tf \cdot (\nabla \times \tau) \circ \Phi = B \circ \Phi, \ \ \ \ \tf|_\Gamma = 0.
\end{equation}

This equation clearly admits a solution $\tf \in C^\infty(\Gamma \times \R,\R)$. Moreover, note that $\p_n$ is transverse to $\Gamma$; that $|\nabla \kappa|$ is uniformly bounded above and below; and that $\kappa \in C^\infty_b$, $\Phi$ is a diffeomorphism from $\Gamma \times (-\delta, \delta)$ to its range, which contains a neighborhood of the form $U_\eta$. Therefore, $\tf$ induces the solution $f = \tf \circ \Phi^{-1} \in C^\infty(U_\eta,\R)$ of (<ref>).

Because $f$ vanishes on $\Gamma$ and $\kappa$ vanishes transversely on $\Gamma$ (with $|\nabla \kappa|$ bounded below), after potentially reducing $\eta$ we can write $f = \beta \kappa$ for some smooth function $\beta$ on $U_\eta$. From the equation (<ref>) and the identity (<ref>), we conclude that

\begin{equation}
\nabla \times A = B = \p_n f + f \nabla \times \tau = \nabla \times (f \tau) = \nabla \times (\beta \kappa \tau).
\end{equation}

Moreover, using that $f$ and $\kappa$ vanish on $\Gamma$, we have along $\Gamma$:

\begin{equation}
B = \p_n f = \p_n (\beta \kappa) = \beta \p_n \kappa = \beta |\nabla \kappa|.
\end{equation}

This completes the proof.

Thanks to (<ref>), $A$ and $\tilde A$ give rise to the same magnetic field on $U_\eta$, hence the difference $A-\tilde A$ is locally a gradient field:

Let $\eta > 0$ given by Lemma <ref> and $U \subset U_\eta$ be a simply connected set. For any $(y_0,\eta_0) \in U \times \R$, there exists a unique $\chi_0 \in C^\infty(U,\R)$ such that

\begin{equation}\label{eq-2p}
\tilde{A} = A -\nabla \chi_0 \quad \text{ in } U; \quad \text{ and } \quad \chi_0(y_0) = \eta_0.
\end{equation}

Thanks to (<ref>), we have $\nabla \times (A -\beta \kappa \tau) = 0$ on $U_\eta$, hence on $U$. Since $U$ is simply connected, by Poincaré's lemma there exists a unique $\chi_0 \in C^\infty(U,\R)$ such that $\nabla \chi_0 = A - \beta \kappa \tau$ and $\chi(y_0)=\eta_0$.

It follows from Lemma <ref> that the operators $\Di$ and $\tDi$ defined in (<ref>) are locally conjugate: on open sets $U$ produced by Lemma <ref>, we have

\begin{equation}\label{eq-3d}
e^{-i\chi_0/\epsi} \Di e^{i\chi_0/\epsi} = \tDi.
\end{equation}

§.§ Global gauge

When $\Gamma$ is simply connected, we can pick $U = U_\eta$ in Lemma <ref>, and the conjugation relation (<ref>) holds on a full neighborhood of $\Gamma$. This however fails when $\Gamma$ is a loop: $\chi_0$ is only defined on part of $\Gamma$. We circumvent this obstacle by using instead a time-dependent gauge that follows the center of mass $y_t$ of our wavepacket, i.e. defined on a set of the form

\begin{equation}
\Omega_\delta = \big\{(t,x), \  t\in \R, \ |x-y_t| < \delta \big\}.
\end{equation}

Let $\tilde A$ be as defined in Lemma <ref>. There exist $\delta > 0$ and $\chi \in C^\infty(\R \times \R^2, \R)$ with support in $\Omega_{2\delta}$ with $\chi(0,y_0) = 0$ and such that for $(t,x) \in \Omega_\delta$:

\begin{equation}\label{eq-2h}
\tilde{A}(x) = A(x) - \nabla \chi(t,x), \qquad \p_t \chi(t,x) = 0.
\end{equation}

As in the proof of Lemma <ref>, it suffices to construct $\chi$ on $\Omega_\delta$ and to extend it to $\R^2$ as a smooth function with support in $\Omega_{2\delta}$. Without loss of generality, we can assume that $\Gamma$ is connected. If $\Gamma$ is also simply connected, then we simply take $\delta = \eta$ and $\chi(t,x) = \chi_0(x)$. If $\Gamma$ is not simply connected, then it is a loop; in particular it is compact.

Fix $t \in \R$. According to Lemma <ref>, there exists $\delta_t > 0$ and a smooth function $\chi(t,\cdot)$ defined on the ball $B(y_t,\delta_t)$ such that

\begin{equation}\label{eq-2k}
\tilde{A}(x) = A(x) - \nabla \chi(t,x), \quad x \in B(y_t,\delta_t); \qquad \chi(t,y_t) = \int_0^t \dot{y_s} \cdot A(y_s) ds.
\end{equation}

Since $\Gamma$ is compact, we can pick $\delta_t$ independent of $t$; we write below $\delta = \delta_t$.  Varying $t$, we obtain a uniquely defined function $\chi$ on $\Omega_\delta$. We take time-derivative of both identities in (<ref>). The first one yields $\nabla \p_t \chi(t,x) = 0$ for $x \in B(y_t,\delta)$. The second one gives

\begin{equation}\label{eq-2n}
\p_t \chi(t,y_t) = \p_t \big(\chi(t,y_t) \big) - \dot{y_t} \cdot \nabla \chi(t,y_t)  %=\dfrac{d}{dt} \int_0^t \dot{y_s} \cdot A(y_s) ds + \dot{y_t} \cdot \nabla \chi(t,y_t)
= \dot{y_t} \cdot A(y_t) - \dot{y_t} \cdot A(y_t) = 0,
\end{equation}

where we used $\nabla \chi = \beta \kappa \tau - A$ hence $\nabla \chi(t,y_t) = A(y_t)$. From these identities, we deduce that $\p_t \chi(t,x) = 0$ for $x \in B(y_t,\delta)$; in particular $\p_t \chi$ is smooth. Moreover, $\nabla \chi = A-\tilde A$ is also smooth. This implies that $\chi$ is smooth, which completes the proof.

Thanks to (<ref>), we have the relation

\begin{equation}\label{eq-2o}
e^{-i\chi/\epsi} \big(\epsi D_t + \Di\big) e^{i\chi/\epsi} = \epsi D_t + \tDi + R, \qquad \text{where:}
\end{equation}

\begin{equation} \label{eq:R}
R= \p_t\chi + (\p_1\chi-A_1+\tilde A_1)\sigma_1 + (\p_2\chi-A_2+\tilde A_2)\sigma_2.
\end{equation}

It follows from (<ref>) that $\Psi$ solves (<ref>) if and only if $\tilde \Psi = e^{-i\chi/\epsi}\Psi$ solves

\begin{equation}\label{eq:D2}
\big( \epsi D_t + \tilde {\slashed D} + R\big)  \tilde \Psi = 0.
\end{equation}

While $R$ looks like a leading-order term, it will effectively be of order $\eps^\infty$ because it vanishes on $\Omega_\delta$ – a domain where our wavepacket is concentrated. We will eventually treat it as a small source term in <ref>.

§.§ Gauge expansion
We conclude this section with an expansion of $\chi$ near $(t,y_t)$, which will serve to prove Theorem <ref>.

We have:

\begin{equation}\label{eq-2g}
\chi(t,y_t) = \int_0^t \dot{y_s} \cdot  A_s ds; \qquad \nabla \chi(t,y_t) = A_t; \qquad \nabla^2 \chi (t,y_t) = \nabla A^\top(y_t) - B_t n_t \tau_t^\top.
\end{equation}

The first formula of (<ref>) comes from the equation (<ref>) defining $\chi$. The second one follows from

\begin{equation}
\nabla \chi(t,y_t) = A(y_t)-\beta(y_t) \kappa(y_t) \tau(y_t) = A_t.
\end{equation}

In order to see the last one, we notice that

\begin{align}
\nabla^2 \chi (t,y_t) & = \nabla \cdot \big(\nabla^\top \chi  \big)(y_t) = \nabla \big( A^\top - \beta \kappa \tau^\top \big)(y_t)
\\
& = \nabla A^\top(y_t) - \beta(y_t) \nabla \kappa(y_t) \tau^\top (y_t) = \nabla A^\top(y_t) - B_t n_t \tau_t^\top,
\end{align}

where we used $\beta_t =B_t/r_t$ on $\Gamma$, see (<ref>). This completes the proof.

§ GAUGE TRANSFORMATION

Our first lemma constructs a magnetic potential $\tilde A$ equivalent to $A$ (in the sense that $\nabla \times \tilde A = B$) that (a) vanishes along $\Gamma=\kappa^{-1}(0)$; (b) is carried by the vector field $\tau$. We use below the notation

\begin{equation}
U_\eta = \big\{ x \in \R^2 : \ d(x,\Gamma) < \eta\big\} .
\end{equation}

There exist $\eta > 0$ and a function $\beta \in C^\infty(U_\eta,\R)$ such that

\begin{equation}\label{eq-0h}
\nabla \times (A - \beta \kappa \tau) = 0 \quad \text{ on } U_\eta; \qquad \beta = \dfrac{B}{|\nabla \kappa|} \quad \text{ on } \Gamma.
\end{equation}

For $f$ a smooth function, we observe that

\begin{equation}\label{eq-0w}
\nabla \times (f \tau) = \nabla f \times \tau +f \nabla \times \tau = \p_n f + f \nabla \times \tau,
\end{equation}

where we used $\nabla f \times \tau = \p_1 f \tau_2 - \p_2 f \tau_1 = \p_n f$. Let $B = \nabla \times A$. To find $\beta$ such that $B = \nabla \times (\beta \kappa \tau)$, we first solve

\begin{equation}\label{eq-0y}
\p_n f + f  \nabla \times \tau = B, \ \ \ \ f|_\Gamma = 0.
\end{equation}

A unique solution $f \in C^\infty(U,\R)$ exists on a neighborhood $U$ of $\Gamma$ in $\R^2$ because the vector field $\p_n$ is transverse to $\Gamma$. Since $f$ vanishes on $\Gamma$ and $\kappa$ vanishes transversely on $\Gamma$, after potentially shrinking $U$ we can write $f = \beta \kappa$ for some smooth function $\beta$ on $U$. From the equation (<ref>) and the identity (<ref>), we conclude that

\begin{equation}
\nabla \times A = B = \p_n f + f \nabla \times \tau = \nabla \times (f \tau) = \nabla \times (\beta \kappa \tau).
\end{equation}

Moreover, using that $f$ and $\kappa$ vanish on $\Gamma$, we have along $\Gamma$:

\begin{equation}
B = \p_n f = \p_n (\beta \kappa) = \beta \p_n \kappa = \beta |\nabla \kappa|.
\end{equation}

This completes the proof.

We now pick $\tilde A \in C^\infty(\R^2,\R^2)$ such that $\tilde A  = \beta \kappa \tau$ on $U$, and we define the operator

\begin{equation}\label{eq-0j}
\tilde{\slashed D} = (\eps D_1-\tilde A_1(x))\sigma_1+(\eps D_2-\tilde A_2(x))\sigma_2 + \kappa(x)\sigma_3
\end{equation}

acting on functions on $\R^2$. Thanks to (<ref>), $A$ and $\tilde A$ give rise to the same magnetic field on $U$: $\nabla \times A = \nabla \times \tilde{A} = B$. This implies the operators $\slashed D$ and $\tilde{\slashed D}$ are locally conjugate via a gauge function:

Let $U$ be given by Lemma <ref> and $U_0 \subset U$ be a simply connected set. For any $(y_0,\eta_0) \in U_0 \times \R$, there exists $\chi_0 \in C^\infty(U_0,\R)$ such that

\begin{equation}\label{eq-2p}
\tilde{A} = A -\nabla \chi_0 \quad \text{ in } U_0; \quad \text{ and } \chi_0(y_0) = \eta_0.
\end{equation}

Thanks to (<ref>) and to $\tilde A  = \beta \kappa \tau$ on $U$, we have $\nabla \times (A -\tilde{A}) = 0$ on $U$, hence on $U_0$. Since $U_0$ is simply connected, by Poincaré's lemma there exists $\chi_0 \in C^\infty(U_0,\R)$ such that $\nabla \chi_0 = \tilde{A} - A$. This is (<ref>) up to a shift of $\chi_0$ by a constant.

§.§ Global time-dependent gauge.

When $\Gamma$ is simply connected, we can pick $U_0$ to be a fully simply connected, open neighborhood of $\Gamma$. Lemma <ref> produces then $\chi_0$ on a full neighborhood of $\Gamma$. However, when $\Gamma$ is a loop, $\chi_0$ is only defined on part of $\Gamma$ – in contrast to $\tilde A$ and $\beta$. Since we aim to construct solutions to $(\epsi D_t + \Di) \psi = 0$ concentrated at time $t$ near $y_t$, we go around this obstruction by constructing a time-dependent gauge that (a) is defined on a neighborhood of $\{(t,y_t) : t \in \R\}$; (b) is locally time-independent; and (c) transforms $A$ to $\tilde{A}$.

Let $y_0\in\Gamma$. There exists an open neighborhood $\Omega$ of $\{(t,y_t), t\in \R\}$ in $\R \times U$ and a smooth function $\chi \in C^\infty(\Omega, \R)$ with $\chi(0,y_0) = 0$ and such that for $(t,x) \in \Omega$:

\begin{equation}\label{eq-2h}
\tilde{A}(x) = A(x) - \nabla \chi(t,x), \qquad \p_t \chi(t,x) = 0.
\end{equation}

Fix $t \in \R$. According to Lemma <ref>, there exists $\delta_t > 0$ and a smooth function $\chi(t,\cdot)$ defined on the ball $B(y_t,\delta_t)$ such that

\begin{equation}\label{eq-2k}
\tilde{A}(x) = A(x) - \nabla \chi(t,x), \quad x \in B(y_t,\delta_t); \qquad \chi(t,y_t) = \int_0^t \dot{y_s} \cdot A(y_s) ds.
\end{equation}

Now varying $t$, we obtain a function $\chi$ on $\Omega = \big\{ (t,x) : |x-y_t| < \delta_t \big\}$, which is clearly a neighborhood of $\Gamma$. We take time-derivative of both identities in (<ref>). The first one yields $\nabla \p_t \chi(t,x) = 0$ for $x \in B(y_t,\delta_t)$. The second one gives

\begin{equation}\label{eq-2n}
\p_t \chi(t,y_t) = \p_t \big(\chi(t,y_t) \big) - \dot{y_t} \cdot \nabla \chi(t,y_t)  %=\dfrac{d}{dt} \int_0^t \dot{y_s} \cdot A(y_s) ds + \dot{y_t} \cdot \nabla \chi(t,y_t)
= \dot{y_t} \cdot A(y_t) - \dot{y_t} \cdot A(y_t) = 0,
\end{equation}

where we used $\nabla \chi = \beta \kappa \tau - A$ hence $\nabla \chi(t,y_t) = A(y_t)$. From these identities, we deduce that $\p_t \chi(t,x) = 0$ for $x \in B(y_t,\delta_t)$; in particular $\p_t \chi$ is smooth. Moreover, $\nabla \chi = A-\tilde A$ is also smooth. This implies that $\chi$ is smooth, which completes the proof.

So far, the constructions of $\beta(x)$, $\tilde A(x)=(\beta\kappa\tau)(x)$ and $\chi(t,x)$ are local in the vicinity of $\Gamma$. Since the wavepacket $\Psi(t,x)$ typically extends over a larger region, we need to extend the above constructions to $\Rm^2$.

We choose $\eta>0$ sufficiently small such that $U_{2\eta}$, a $2\eta-$neighborhood of $\Gamma$, is a tubular neighborhood of $\Gamma$ (so that each point in $U_{2\eta}$ is uniquely parametrized by a coordinate along $\Gamma$ and one across $\Gamma$). The existence of $\eta>0$ is guaranteed by the transversality condition $\inf\{ \big|\nabla \kappa(x)\big| : x \in \Gamma\} >0$. Let $U_\eta$ be the $\eta-$neighborhood of $\Gamma$. Then we have the following result:

Let $y_0\in\Gamma$. There exists a function $\chi(t,x)\in C^\infty(\Rm\times \Rm^2)$ with $\chi(0,y_0)=0$, a vector field $\tilde A(x)\in C^\infty(\Rm^2,\Rm^2)$ such that $A(x)-\tilde A(x)$ is compactly supported, and a function $\beta\in C^\infty(\Rm^2)$ such that for all $(t,x)\in \Rm\times U_\eta \cap \{ |x-y_t|<\eta\}$, we have

\begin{equation}\label{eq:globalgauge}
(\beta\kappa\tau)(x)=  \tilde A(x)=A(x)-\nabla\chi(t,x),\qquad \partial_t\chi(t,x)=0.
\end{equation}

I do not think we have $\partial_t\chi(t,x)=0$ on $(t,x)\in\Rm\times U_\eta$? For instance, if $\Gamma$ is a circle then integrating over $t$ we would get $\chi(0,x) = \chi(t,x)$ for all $t$, so in particular if $y_T = y_0$ we would have $\chi(0,y_0) = \chi(T, y_T)$, which we know is wrong? The $(t,x)$ domain needs to be restricted to the sliding domain of course, I agree.

We may choose in the proof of the preceding lemma $\min_{t\in \Rm} \delta_t=2\eta>0$.   Let $\beta_\eta(x)$, $\tilde A_\eta=\beta_\eta\kappa\tau$ and $\chi_\eta(t,x)$ be the corresponding functions constructed in Lemmas <ref>, <ref>, and <ref> on $U_{2\eta}$.

Let $\phi_\eta(x)\in C^\infty(\Rm^2)$ be equal to $1$ in $U_\eta$, the $\eta$-neighborhood of $\Gamma$, and equal to $0$  on $\Rm^2\backslash U_{2\eta}$.

We then define $\beta(x)=\phi_\eta(x)\beta_\eta(x)$, $\tilde A(x)=\phi_\eta(x)\tilde A_\eta(x)+(1-\phi_\eta(x))A(x)$, and $\chi(t,x)=\phi_\eta(x)\chi_\eta(t,x)$. The above construction ensures that $\chi$ and $A-\tilde A$ are compactly supported with support in the vicinity of $\Gamma$.

We observe that $\tilde A(x)=A(x)-\nabla \chi(t,x)$ on the domain $\{|x-y_t|<\eta\}$ but not (necessarily) on $\Rm^2$. We still define the conjugate Dirac operator $\tilde\Di$ by (<ref>), but now for the globally defined magnetic potential $\tilde A$.

We then observe that

\begin{equation}\label{eq-2o}
e^{-i\chi(t,x)/\epsi} (\epsi D_t + \Di) e^{i\chi(t,x)/\epsi}  = \epsi D_t + \tilde{\Di} + \delta_\eta,
\end{equation}

where the error term $\delta_\eta$ is the multiplication operator by the following matrix-valued function

\begin{equation}\label{eq:deltaeta}
\delta_\eta = (iD_t\chi) + (iD_1\chi-A_1+\tilde A_1)\sigma_1 + (iD_2\chi-A_2+\tilde A_2)\sigma_2.
\end{equation}

The construction of $\tilde A(x)$ and $\chi(t,x)$ in (the proof of) Proposition <ref> ensures that $\delta_\eta$ is smooth and compactly supported away from $\Gamma$ in space.

It follows from (<ref>) that $\Psi$ is a solution of (<ref>) if and only if $\tilde \Psi(t,x) = e^{-i\chi(t,x)/\epsi}\Psi(t,x)$ is a solution of

\begin{equation}\label{eq:D2}
\big( \epsi D_t + \tilde {\slashed D} + \delta_\eta\big)  \tilde \Psi = 0.
\end{equation}

Whereas $\delta_\eta(t,x)$ in (<ref>) appears to be a leading-order term, it vanishes on $U_\eta$ at all times and hence on a domain where most of our wavepacket is concentrated. It is therefore effectively of order $\eps^\infty$. It will be neglected in our subsequent asymptotic expansions and ultimately treated as a small source term in the error analysis presented in <ref>. Sections <ref> to <ref> construct wavepackets $\tilde\Psi$ that solve $(\epsi D_t + \tilde {\slashed D})\tilde\Psi=0$.

§.§ Gauge expansion
Since our wavepacket will be concentrated within $\sqrt{\epsi}$ of $y_t$, eventually only the values of $\chi$ and its derivatives at $(t,y_t)$ will matter. We compute them at order up to two (which is enough to derive the leading order expansion of the wavepacket):

We have:

\begin{equation}\label{eq-2g}
\chi(t,y_t) = \int_0^t \dot{y_s} \cdot \dot A_s ds; \qquad \nabla \chi(t,y_t) = A_t; \qquad \nabla^2 \chi (t,y_t) = \nabla A^\top(y_t) - B_t n_t \tau_t^\top.
\end{equation}

The first formula of (<ref>) comes from the equation (<ref>) defining $\chi$. The second one follows from

\begin{equation}
\nabla \chi(t,y_t) = A(y_t)-\beta(y_t) \kappa(y_t) \tau(y_t) = A_t.
\end{equation}

In order to see the last one, we notice that

\begin{align}
\nabla^2 \chi (t,y_t) & = \nabla \cdot \big(\nabla^\top \chi  \big)(y_t) = \nabla \big( A^\top - \beta \kappa \tau^\top \big)(y_t)
\\
& = \nabla A^\top(y_t) - \beta(y_t) \nabla \kappa(y_t) \tau^\top (y_t) = \nabla A^\top(y_t) - B_t n_t \tau_t^\top,
\end{align}

where we used $\beta_t =B_t/r_t$ on $\Gamma$, see (<ref>). This completes the proof.

§ LOCAL NORMAL FORM

We now start our construction of approximate solutions to $(\epsi D_t + \tilde {\slashed D})\tilde\Psi=0$. This section describes a series of transformations that allows us to explicitly invert an operator that governs the leading dynamics: a rescaling of $\tilde \Psi$ in natural coordinates, spatial and spinorial rotations, and finally a shifted partial Fourier transform.

§.§ Spatial rescaling
We first write the wavepacket in natural coordinates:

\begin{equation}\label{eq:psiz}
\psi(t,z)%=\eps^{\frac12}(Se^{-i\frac1\eps\chi}\psi_M)(t,z)
= \eps^{\frac12} S \tilde \Psi (t,z),
\end{equation}

with the scaling transformation defined as

\begin{equation}\label{eq:S}
Sf(z)=f(y_t+\sqrt\eps z),\qquad S^{-1}f(x)=f\left(\frac{x-y_t}{\sqrt\eps}\right).
\end{equation}

We observe that $\eps^{\frac12}S$ is an isometry on $L^2(\Rm^2)$, while formally,

\begin{equation}\label{eq:Sj}
Sf(z) = \sum_{\alpha} \frac{\eps^{\frac{|\alpha|}2}z^\alpha}{\alpha!} \partial^\alpha f(y_t) = \dsum_{j\geq0} \eps^{\frac j2} S_jf(z),\quad S_jf(z):= \dsum_{|\alpha|=j} \frac{z^\alpha}{\alpha!} \partial^\alpha f(y_t).
\end{equation}

This scaling emerges from the following consideration. For $\kappa(x) = x_2$, the wavepacket is confined by a harmonic-like oscillator $\epsi D_2 \sigma_2 + x_2 \sigma_3$ and hence comes with a natural scale $\sqrt\eps x_2$. This scale is then imposed to all directions when $\kappa$ admits more sophisticated variations, i.e. when $\Gamma$ is curved.

We then verify that $\psi(t,z)$ solves

\begin{equation}\label{eq:L}
L\psi=0 ;\qquad L= \eps^{-\frac12} S(\eps D_t+\tilde{\slashed D})S^{-1}  = \eps^{\frac12} D_t -\dot y_t \cdot D + (D+\eps^{-\frac12} Sh(z))\cdot \sigma
\end{equation}

with $D=(D_1,D_2,0)^t$, $\sigma=(\sigma_1,\sigma_2,\sigma_3)^t$, and $h=(-\tilde A_1,-\tilde A_2,\kappa)^t=(-\beta\kappa\tau,\kappa)^t$. We still use $\dot y_t$ for $(\dot y_t,0)^t\in\Rm^3$.

Using (<ref>), we find formally $L=\sum_{j\geq0} \eps^{\frac j2}L_j$, where we observe that $L_{-1}=0$ since $\tilde A$ and hence $h$ vanish on $\Gamma = \kappa^{-1}(0)$, and

\begin{eqnarray}
L_0 &=& -\dot y_t\cdot D + D\cdot\sigma + S_1h \cdot\sigma \label{eq:L0} \\
L_1 &=& D_t + S_2h\cdot\sigma \label{eq:L1}\\
L_j&=& S_{j+1}h\cdot\sigma  = \dsum_{|\alpha|=j+1} \frac{1}{\alpha!} z^\alpha \partial^\alpha h(y_t)\cdot\sigma,\qquad j\geq2 \label{eq:Lj}.
\end{eqnarray}

Our next objective is to transform $L_0$ in an appropriate basis so that its infinite dimensional kernel and its inverse on the orthogonal complement may be written explicitly.

§.§ Rotations
We now introduce spatial and spinorial rotations already mentioned in the introduction. We define, for $1\leq j\leq 3$ and $\theta\in\Rm$, the spinor rotation acting on $\Cm^2$:

\[
U_{j,\theta} = e^{-i\frac\theta2 \sigma_j} = \cos \left(\tfrac\theta2\right) -i \sin\left(\tfrac \theta2\right) \sigma_j.
\]

Associated to it is a standard three-dimensional rotation acting on $\Rm^3$ of angle $\theta$ about the $j$-axis (with usual orientation) $\tilde R_{j,\theta}$ such that for $c \in \mathbb R^3$:

\begin{equation}\label{eq:rotations}
\tilde R_{j,\theta} c \cdot \sigma = c\cdot U_{j,\theta}^*\sigma U_{j,\theta},\qquad U_{j,\theta}^*\sigma U_{j,\theta} =: \tilde R_{j,\theta}^*\sigma.
\end{equation}

Note that

\begin{equation}\label{eq:spinrot}
U_{k,\theta}^*\sigma_j U_{k,\theta} = \delta_{j\not=k}\cos\theta \sigma_j+ \epsilon_{jkl}\sin\theta \sigma_l+ \delta_{jk}\sigma_j.
\end{equation}

From this, we deduce two formulas that we will use later:

\begin{equation}\label{eq:rotations23}
\tilde R_{2,\varphi_t} = \matrice{ \cos\varphi_t &0& -\sin\varphi_t \\ 0&1&0 \\ \sin\varphi_t &0&\cos\varphi_t} = \matrice{ c_t & 0 & -s_t \\ 0 & 1 & 0 \\ s_t & 0 & c_t}, \qquad
\tilde R_{3,\theta_t} = %\matrice{ \cos\theta_t & \sin\theta_t&0 \\ -\sin\theta_t & \cos\theta_t &0\\0&0&1 } =
\matrice{R_{\te_t} & 0 \\ 0 & 1},
\end{equation}

where $R_{\te_t}$ is the spatial rotation defined in (<ref>).
We finally define the $L^2$-unitary transform related to the spatial rotations

\begin{equation}\label{eq:mR}
\mR_{\theta_t} f=f\circ R_{\theta_t}
\end{equation}

and the operator

\begin{equation}\label{eq:bUt}
\bU_t = U_{2,\varphi_t} U_{3,\theta_t} \mR_{\te_t}.
\end{equation}

When $\varphi_t\equiv0$ (i.e., when $B=0$), this transformation already appears in <cit.>.
We now compute the operators $\bU_t^* L_j \bU_t$. We start with preliminary relations:

We find:

\begin{equation}\label{eq-0m}
\bU_t^* ( -\dot  y_t \cdot D) \bU_t = c_t D_1,
\quad
\bU_t^* D\cdot\sigma \bU_t  = (c_t\sigma_1+s_t\sigma_3) D_1 + \sigma_2 D_2,
\quad
\end{equation}

\begin{equation}\label{eq-0q}
\bU_t^* D_t \bU_t
= D_t-\frac{\dot\varphi_t}{2} \sigma_2 - \frac{\dot\te_t}{2} (-s_t \sigma_1 + c_t \sigma_3) - \dot{\theta_t} (z_1 D_2 - z_2 D_1).
\end{equation}

Moreover $\bU_t^* S_j h \cdot \sigma \bU_t$ is a multiplication operator by a time-dependent homogeneous polynomial of degree $j$ in $z$; and for $j=1,2$:

\begin{equation}\label{eq-0o}
\bU_t^* S_1 h \cdot \sigma \bU_t = \rho_t z_2 \sigma_3, \quad \bU_t^* S_2 h \cdot \sigma \bU_t = -r_t c_t z_2 \lr{z, \nabla \beta(y_t)} + q_t(z) \cdot \sigma
\end{equation}

where $q_t\cdot\sigma$ is a matrix-valued quadratic form in $z$ carried by $\sigma_2, \sigma_3$.

1. Space differentiation. We first observe that

\begin{equation}\label{eq-0n}
\mR_{\theta_t}^*D \mR_{\theta_t} = \tilde R_{3,{\theta_t}}^*D,\quad \mR_{\theta_t}^*D\cdot v \mR_{\theta_t} = D\cdot \tilde R_{3,{\theta_t}} v.
\end{equation}

Applied to $v=-\dot y_t$ with $R_{\theta_t} \dot y_t=-c_te_1$ by construction, we find $-\tilde R_{3,{\theta_t}} \dot y_t \cdot D=c_t D_1$. Since spinorial rotations commute with scalars, we obtain the first relation in (<ref>).

We now apply (<ref>) to $v=\sigma$ and find

\[
(\mR_{\theta_t} U_{3,{\theta_t}})^* D\cdot\sigma \mR_{\theta_t} U_{3,{\theta_t}} = D\cdot\sigma.
\]

Therefore, using (<ref>), we find

\[
\bU_t^* D\cdot\sigma \bU_t = U_{2,\varphi_t}^* D\cdot\sigma U_{2,\varphi_t} = D \cdot \tilde{R}_{2,\varphi_t}^* \sigma = (c_t\sigma_1+s_t\sigma_3) D_1 + \sigma_2 D_2.
\]

This proves the second identity in (<ref>).

2. Time differentiation. We note that

\begin{equation}
\mR_{\theta_t}^*D_t\mR_{\theta_t}=D_t-\dot{\theta_t} Jz \cdot D_z,
\quad
U_{3,\te_t}^*D_t U_{3,\te_t} = D_t-\frac{\dot\te_t}{2} \sigma_3,
\quad
U_{2,\varphi_t}^*D_t U_{2,\varphi_t} = D_t-\frac{\dot\varphi_t}{2} \sigma_2.
\end{equation}

We deduce that

\begin{align}
\bU_t^* D_t \bU_t & = U_{2,\varphi_t}^* \Big( D_t-\frac{\dot\te_t}{2} \sigma_3 \Big) U_{2,\varphi_t} - \dot{\theta_t}  (z_1 D_2 - z_2 D_1)
\\
& = D_t-\frac{\dot\varphi_t}{2} \sigma_2 - \frac{\dot\te_t}{2} (-s_t \sigma_1 + c_t \sigma_3) - \dot{\theta_t}  (z_1 D_2 - z_2 D_1).
\end{align}

This proves (<ref>).

3. Multiplicative operators.  We observe that $\mR^*_{\theta_t} Sh \mR_{\theta_t}= h(y_t+\sqrt\eps R_{\theta_t}z)$ as a multiplication operator so that

\[
\bU_t^* Sh\cdot\sigma \bU_t= \tilde R_{2,\varphi_t} \tilde R_{3,\theta_t} h (y_t+\sqrt\eps R_{\theta_t}z) \cdot\sigma =\sum_{j\geq1} \eps^{\frac j2} \dsum_{|\alpha|=j} z^\alpha \nu_{\alpha} \cdot \sigma,
\]

where we have defined the $3-$vectors $\nu_\alpha = \frac{1}{\alpha!} \tilde R_{2,\varphi_t}\tilde R_{3,\theta_t} (R_{-\theta_t}\nabla)^\alpha h(y_t)$. Indeed,

\[
h(y_t+\sqrt \eps R_\theta z) %= \dsum_{j\geq1} \frac{1}{j!} (\sqrt\eps R_\theta z\cdot\nabla)^j h(y_t)
= \dsum_{j\geq1} \frac{1}{j!} (\sqrt\eps z\cdot R_{-\theta} \nabla)^j h(y_t) = \dsum_{j\geq1} \eps^{\frac j2} \dsum_{|\alpha|=j} \frac{z^\alpha}{\alpha!} (R_{-\theta}\nabla)^\alpha h(y_t)
\]

using multinomial coefficients.  Therefore, $S_j h \cdot \sigma$ is a multiplicative operator by a homogeneous polynomial of degree $j$ depending on $t$;  composing with rotations preserves this feature: $\bU_t^* S_j h \cdot \sigma \bU_t$ is also a homogeneous polynomial of degree $j$.

We now focus on $j=1$. Define $\tilde h = (-\beta \tau,1)^t$ so that $h = \kappa \tilde h$.
Since $\kappa(y_t) =0$, we have $S_1 h = \tilde{h}(y_t) \cdot S_1 \kappa$. Since spinorial rotations commute with the scalar $S_1 \kappa$,

\begin{equation}\label{eq-0p}
\bU_t^* S_1 h \cdot \sigma \bU_t = \RR_{\te_t}^* S_1 \kappa \RR_{\te_t} \  \bU_t^*\tilde{h}(y_t) \cdot  \sigma   \bU_t.
\end{equation}

We now compute separately $\mR_{\theta_t}^* S_1 \kappa \mR_{\te_t}$ and $\bU_t^*\tilde{h}(y_t) \cdot  \sigma   \bU_t$ (note that $\RR_{\te_t}$ does not technically affect that last term since $\tilde{h}(y_t)$ does not depend on $z$). We have:

\begin{equation}\label{eq-1g}
\mR_{\theta_t}^* S_1 \kappa \cdot \sigma \mR_{\te_t} = \mR_{\theta_t}^* z \cdot \nabla \kappa(y_t) \mR_{\te_t} = z \cdot R_{\te_t}\nabla \kappa(y_t) = r_t z_2.
\end{equation}

Moreover, since $B_t = \beta_t r_t$, we have $\tilde h (y_t) = (-r_t^{-1} B_t \tau_t, 1)^t$ hence $R_{3,\te_t} \tilde{h} = (r_t^{-1}B_t,0,1)$. We deduce

\begin{equation}
U_{3,\te_t}^* \tilde{h}(y_t) \cdot  \sigma   U_{3,\te_t} =  R_{3,\te_t} \tilde h(y_t) \cdot \sigma  = \dfrac{B_t \sigma_1}{r_t} + \sigma_3 = \dfrac{\rho_t}{r_t} (s_t \sigma_1 + c_t \sigma_3).
\end{equation}

By (<ref>), $U_{2,\varphi_t}^* (s_t \sigma_1 + c_t \sigma_3) U_{2,\varphi_t} = \sigma_3$, which is the main motivation for the definition of $U_{2,\varphi_t}$. Hence,

\begin{equation}\label{eq-1h}
\bU_t^*\tilde{h}(y_t) \cdot  \sigma  \bU_t = U_{2,\varphi_t}^*  U_{3,\te_t}^* \tilde{h}(y_t) \cdot  \sigma   U_{3,\te_t}U_{2,\varphi_t} = \dfrac{\rho_t}{r_t} \sigma_3.
\end{equation}

Multiplying (<ref>) and (<ref>) and going back to (<ref>), we conclude that:

\begin{equation}\label{eq-0c}
\bU_t^* S_1 h \cdot \sigma \bU_t = \rho_t z_2 \sigma_3.
\end{equation}

We go on with the case $j=2$. We compute $\bU_t^* S_2 h \cdot \sigma \bU_t$. We recall that $h = \kappa \tilde h$, therefore (using again $\kappa(y_t) = 0$):

\begin{equation}\label{eq-0u}
S_2 h = 2 (S_1 \kappa) (S_1 \tilde{h})  + (S_2 \kappa) \tilde{h}(y_t).
\end{equation}

We start with the contribution of $(S_2 \kappa) \tilde{h}(y_t)$, i.e. the term $\bU_t^* (S_2 \kappa) \tilde{h}(y_t)  \cdot \sigma \bU_t$. Using that $S_2 \kappa$ is a scalar and $\tilde h \cdot \sigma$ does not depend on $z$ we have

\begin{equation}\label{eq-0v}
\bU_t^* (S_2 \kappa) \tilde{h}(y_t) \bU_t = \RR_{\te_t}^* (S_2 \kappa) \RR_{\te_t} \ \bU_t^* \tilde{h}(y_t) \cdot  \sigma \bU_t = \RR_{\te_t}^* (S_2 \kappa) \RR_{\te_t} \dfrac{\rho_t}{r_t} \sigma_3,
\end{equation}

where we used (<ref>) in the last equality. This term is carried by $\sigma_3$.

We now focus on $\bU_t^* (S_1 \kappa) (S_1 \tilde{h}(y_t) \cdot \sigma) \bU_t$.
Using (<ref>) and that spinorial rotations commute with scalars, we obtain:

\begin{equation}\label{eq-1i}
\bU_t^* (S_1 \kappa) (S_1 \tilde{h}) \cdot \sigma \bU_t = \mR_{\te_t}^* (S_1 \kappa) \mR_{\te_t}\, \bU_t (S_1 \tilde{h}) \cdot \sigma \bU_t = r_t z_2 \, \bU_t (S_1 \tilde{h}) \cdot \sigma \bU_t.
\end{equation}

It remains to compute $\bU_t (S_1 \tilde{h}) \cdot \sigma \bU_t$. We recall that $\tilde h = (-\beta \tau,1)^t$. Moreover, $\tau$ has unit norm: $\lr{\tau, \tau} = 1$. Taking derivatives of this expression shows that $\p_j \tau$ is normal to $\tau$; we write $\p_j \tau = \alpha_j n$ below. We deduce that

\begin{align}
(S_1 \tilde{h}) \cdot \sigma & = -\sum_{j,k=1}^2 z_j \p_j (\beta \tau_k) \sigma_k = -\sum_{j,k=1}^2 z_j (\p_j \beta) \tau_k \sigma_k - \beta \sum_{j,k=1}^2 z_j \p_j \tau_k \sigma_k
\\ & = -(S_1\beta) \tau \cdot \sigma - \beta \sum_{j,k=1}^2 z_j \alpha_j n_k \sigma_k = -(S_1\beta) \tau \cdot \sigma - \lr{z,\alpha} n \cdot \sigma.
\end{align}

We observe that:

\begin{equation}
U_{3,\theta_t}^*  \tau \cdot \sigma  U_{3,\theta_t} =  R_{3,\theta_t} \tau \cdot \sigma = -\sigma_1, \qquad U_{2,\varphi_t}^*  \sigma_1  U_{2,\varphi_t} = R_{2,\varphi_t} e_1 \cdot \sigma = c_t \sigma_1 - s_t \sigma_3;
\end{equation}

\begin{equation}
U_{3,\theta_t}^*  n \cdot \sigma  U_{3,\theta_t} =  R_{3,\theta_t} n \cdot \sigma = \sigma_2, \qquad U_{2,\varphi_t}^*  \sigma_2  U_{2,\varphi_t} = \sigma_2.
\end{equation}

Therefore, the following equality is valid modulo terms carried by $\sigma_2, \sigma_3$:

\begin{equation}
\bU_t (S_1 \tilde{h}) \cdot \sigma \bU_t = \RR_{\te_t}^* S_1\beta \RR_{\te_t} c_t \sigma_1 = c_t \lr{z, R_{\te_t} \nabla \beta(y_t)} \sigma_1.
\end{equation}

Going back to  (<ref>), (<ref>) and (<ref>), we conclude that  – again with equality valid  modulo terms carried by $\sigma_2, \sigma_3$:

\begin{align}
\bU_t^* S_2 h \cdot \sigma \bU_t = \bU_t^* (S_1 \kappa) (S_1 \tilde{h}) \cdot \sigma \bU_t = -c_t r_t z_2 \lr{z, R_{\te_t} \nabla \beta(y_t)} \sigma_1.
\end{align}

This completes the proof of (<ref>), hence of the lemma.

Thanks to Lemma <ref>, we can describe  $\bU_t^* L_j \bU_t$. For $j \geq 2$, $\bU_t^* L_j \bU_t$ is a multiplication operator by a homogeneous polynomial of degree $j$. For $j = 0, 1$, we have explicit expressions:

\begin{align}
\label{eq-1s}  \bU_t^* L_0 \bU_t & = c_t (1+\sigma_1) D_1 + \sigma_2 D_2 + (\rho_t z_2 + s_t D_1) \sigma_3;
\\
\label{eq-1r}   \bU_t^* L_1 \bU_t & = D_t - \dot{\te_t} (z_1 D_2 - z_2 D_1) + \left( \dfrac{\dot{\te_t}s_t}{2} + j_t z_1z_2 + k_t z_2^2\right) \sigma_1 + E_1,
\end{align}

where $E_1$ is a multiplication operator by a polynomial of degree two and carried by $\sigma_2$ and $\sigma_3$ (we will see later that it does not contribute to leading order). The coefficients $j_t$ and $k_t$ are explicitly given in terms of $\beta$:

\begin{equation}\label{eq-1c}
\matrice{j_t \\ k_t}
= r_t c_t R_{\te_t} \nabla \beta(y_t) = r_t c_t \matrice{-\p_\tau \beta(y_t) \\ \p_n \beta(y_t)}.
\end{equation}

We end this section with a few  important relations:

\begin{equation}\label{eq-1d}
j_t = -r_t \dot{\beta_t}; \quad
j_t \gamma_t =  \dfrac{d \ln c_t}{dt}; \quad k_t = \dfrac{c_t}{2} \left( \p_n B(y_t) - B_t \dfrac{\Delta \kappa(y_t)}{r_t}\right).
\end{equation}

We start with the first identity in (<ref>). We recall that $\dot{y_t} = c_t\tau_t$. Since $\p_\tau$ is tangent to $\Gamma$ and $\beta = B/|\nabla \kappa|$ along $\Gamma$:

\begin{equation}\label{eq-1m}
j_t = -r_t c_t \p_\tau \beta(y_t) = -r_t c_t \tau_t \cdot \nabla \left( \dfrac{B}{|\nabla \kappa|} \right) (y_t) = -r_t \dot{y_t} \cdot \nabla \left( \dfrac{B}{|\nabla \kappa|} \right) (y_t) = -r_t \dot{\beta_t}.
\end{equation}

Regarding the second identity: using (<ref>) and $\beta_t=B_t/r_t$ on $\Gamma$, we obtain:

\begin{equation}
j_t \gamma_t = - \dfrac{r_t B_t}{r_t^2+B_t^2} \dot{\beta_t} = - \dfrac{\beta_t \dot{\beta_t}}{1+\beta_t^2}  = -\dfrac{1}{2} \dfrac{d}{dt}  \ln\big(1+\beta_t^2\big) = -\dfrac{1}{2} \dfrac{d}{dt}  \ln\frac{1}{c_t^2}= \dfrac{d \ln c_t}{dt}.
\end{equation}

For the last identity, we first note that $|\nabla \kappa| = \p_n \kappa$. Since $f = \beta \kappa$ solves (<ref>):

\begin{equation}
B = \p_n f + f \nabla \times \tau  = \p_n(\beta \kappa) + \beta \kappa \nabla \times \tau  = \kappa\p_n \beta  + \beta \p_n \kappa + \beta \kappa \nabla \times \tau.
\end{equation}

Therefore, since $\kappa$ vanishes on $\Gamma$, we deduce that on $\Gamma$:

\begin{equation}\label{eq-1b}
\p_n B = 2\p_n \kappa \p_n \beta + \beta \big(\p_n^2 \kappa + \p_n \kappa \nabla \times \tau \big) = 2|\nabla \kappa| \p_n \beta + B\left( \dfrac{\p_n^2 \kappa }{|\nabla \kappa|}+\nabla \times \tau\right).
\end{equation}

Moreover, using $|\nabla \kappa| = \p_n \tau$ and $\nabla f\times J\nabla g=\nabla f\cdot\nabla g$:

\begin{equation}\label{eq-1l}
\nabla \times \tau =  \dfrac{\nabla \times J\nabla \kappa}{|\nabla \kappa|} + \nabla \dfrac{1}{|\nabla \kappa|} \cdot J\nabla \kappa = \dfrac{\Delta \kappa}{|\nabla \kappa|} - \dfrac{\p_n^2 \kappa}{|\nabla \kappa|}.
\end{equation}

Plugging (<ref>) into (<ref>) and using the defining equation $\kappa_t=r_tc_t\partial_n\beta(y_t)$  yield (<ref>).

§.§ Shifted Fourier transform
We finally introduce the $L^2-$ unitary transformation

\begin{equation}\label{eq:FS}
\mV_ta(z) =\dfrac{1}{\sqrt{2\pi}} \dint_{\Rm} e^{iz_1\xione} a(\xione,z_2+\gamma_t\xione) d\xione.
\end{equation}

Its inverse is given explicitly by

\begin{equation}\label{eq:FSinv}
\mV_t^* a(\Xitot) = \dfrac{1}{\sqrt{2\pi}} \dint_{\Rm} e^{-iz_1\xione} a(z_1,\zetatwo-\gamma_t\xi) dz_1.
%, \qquad \Xitot=(\xione,\zetatwo).
\end{equation}

Above, $\xione$ is the one-dimensional dual Fourier variable to $z_1$ while the one-dimensional variable $\zetatwo$ is a shift of $z_2$ given by $\zetatwo=z_2+\gamma_t\xione$.

This shifted partial Fourier transform, which reduces to a partial Fourier transform in the first variable $z_1\to\xione$ when $B=0$ and introduces an additional shift $\gamma_t\xione$ to $z_2$ otherwise is reminiscent of the parametrization of Landau level eigenfunctions in a Landau gauge when $B\not=0$.
The construction of our wavepackets involves many functions of the form $a(t,\Xitot)=f(t,\xione)\phi(\zetatwo)$,

with $f$ describing the wavepacket profile along $\Gamma$ while $\phi$ describes the wavepacket profile across $\Gamma$.

The global transformation we need to consider is the composition of $\bU_t$ introduced in the preceding section and $\mV_t$ given above. We thus set $\fU_t = \bU_t \mV_t$. The operator in this new set of variables is $T:= \fU_t^* L \fU_t$ with $L$ defined in (<ref>). We similarly define $T_j = \fU_t^* L_j \fU_t$ for $j\geq0$. We will also be using $D_\xi=-i\partial_\xi$ and $D_\zeta=-i\partial_\zeta$.

In $T_1$, it will turn out that the terms carried by $\sigma_2$ or $\sigma_3$ and the terms that are odd in $(D_\zeta, \zetatwo)$ do not appear to leading order in our expansion. Therefore, we use the notation

\begin{equation}
P_1  \equiv P_2
\end{equation}

if $P_1-P_2$ is a differential operator made of odd terms in $(\zetatwo,D_\zeta)$ and terms carried by $\sigma_2, \sigma_3$.

With  $\fU_t = \bU_t \VV_t$ and $T_j = \fU_t^* L_j \fU_t$, we have:

\begin{align}\label{eq:T0}
T_0 & = c_t (1+\sigma_1) \xione + D_\zeta \sigma_2  + \rho_t\zetatwo \sigma_3,
\\
\label{eq-1t} T_1 & \equiv  D_t + \dot{\te_t} \gamma_t \big(D_\zeta^2-\xione^2\big)  + \Big( \dfrac{\dot{\te_t}s_t}{2}  +  j_t \gamma_t (D_\xi \xione - D_\zeta \zetatwo) + k_t \big(\zetatwo^2 + \gamma_t^2 \xione^2 \big)\Big) \sigma_1.
\end{align}

Moreover $T_j$ is a linear combination  of differential operators of the form $\xi^k \zeta^\ell D_\xi^m D_\zeta^{j+1-k-\ell-m}$, with $k,l\geq0$ and coefficients bounded with respect to $t$.

We note that $\mV_t$ satisfies the canonical relations

\begin{equation}\label{eq:relmVt}
\VV_t^* (z_2 + \gamma_t D_1) \VV_t =  \zetatwo, \quad \VV_t^* D_2 \VV_t = D_\zeta, \quad \VV_t^* D_1 \VV_t = \xione, \quad \mV_t^* z_1 \mV_t = -(D_\xi+\gamma_t D_\zeta).
\end{equation}

We deduce (<ref>) for $T_0 = \fU_t^* L_0 \fU_t$ from the expression  (<ref>). Thanks to (<ref>) we obtain $\VV_t^* z_2 \VV_t = \zetatwo - \gamma_t \xione$, hence using the notation $\equiv$:

\begin{align}
\VV_t^* z_1 z_2 \VV_t & = -(D_\xi+\gamma_t D_\zeta) (\zetatwo - \gamma_t \xione) \equiv  \gamma_t (D_\xi \xione - D_\zeta \zetatwo)
\\
\VV_t^* z_2^2 \VV_t & = (\zetatwo - \gamma_t \xione)^2 \equiv \zetatwo^2 + \gamma_t^2 \xi^2
\\
\VV_t^*  (z_1 D_2 - z_2 D_1) \VV_t & = -(D_\xi+\gamma_t D_\zeta) D_\zeta - (\zetatwo - \gamma_t \xione) \xione \equiv - \gamma_t (D_\zeta^2-\xione^2).
\end{align}

Moreover, we observe that

\begin{equation}
\mV_t^*D_t \mV_t = D_t + \dot \gamma_t \xione D_\zeta \equiv D_t.
\end{equation}

We deduce from (<ref>) the expression for $T_1$:

\begin{align}
T_1 \equiv   D_t + \dot{\te_t} \gamma_t (D_\zeta^2-\xione^2)  + \Big( \dfrac{\dot{\te_t}s_t}{2}  +  j_t \gamma_t (D_\xi \xione - D_\zeta \zetatwo) + k_t \big(\zetatwo^2 + \gamma_t^2 \xione^2 \big) \Big) \sigma_1.
\end{align}

Finally, since $L_j$ is multiplication operator by a homogeneous polynomial of degree $j+1$, $T_j$ is a linear combination  of differential operators of the form $\xi^k \zeta^\ell D_\xi^m D_\zeta^{j+1-k-\ell-m}$, with coefficients bounded with respect to $t$. This completes the proof.

§.§ Conjugation and asymptotic expansion

Introducing $a=\fU_t^*\psi$ and recalling that $T=\fU_t^*L\fU_t$, equation (<ref>) is equivalent to $Ta=0$. With $L_j$ given in (<ref>)–(<ref>) and $T_j=\fU_t^*L_j\fU_t$, we decompose $a$ as $a=\sum_{j\geq0}\eps^{\frac j2}a_j$ to reduce the equation $Ta=0$ to the triangular system

\begin{equation}\label{eq:seqT}
T_0 a_0 =0 ,\quad T_1a_0+T_0a_1=0,\quad \dsum_{k=0}^{j} T_{j-k}a_k =0,\ j\geq2.
\end{equation}

In <ref> we identify the kernel of $T_0$ and we produce explicit formula for its inverse on the orthogonal complement. We describe how to solve the subleading equation $T_1a_0+T_0a_1=0$ in <ref>. We finally present the higher-order asymptotic expansions and corresponding error estimates in <ref>.

§ INVERSION OF LEADING OPERATOR

To study the kernel of the operator

\begin{equation}
T_0 = c_t (1+\sigma_1) \xione + D_\zeta \sigma_2  + \rho_t\zetatwo \sigma_3,
\end{equation}

and invert it on the orthogonal complement, we first bring it to a normal form. Thanks to the change of variables $(\xi,\zeta) \mapsto (\rho_t^{1/2}c_t^{-1}\xinorm,\rho_t^{-1/2}\ynorm)$, we can assume that $\rho_t = c_t = 1$. Moreover, with $Q = Q^*= \frac{1}{\sqrt{2}}(\sigma_1+\sigma_3)$, we have

\begin{equation}
Q \sigma_1 Q=\sigma_3, \quad  Q \sigma_2 Q=-\sigma_2, \quad Q \sigma_3 Q=\sigma_1.
\end{equation}

It suffices then to work with the model operator

\begin{equation}\label{eq:H}
H=\xinorm (1+\sigma_3) + \sigma_1 \ynorm  - D_\ynorm  \sigma_2 = \begin{pmatrix} 2\xi & \fa_\ynorm  \\ \fa_\ynorm ^* & 0\end{pmatrix} ,\quad \fa_\ynorm  = \partial_\ynorm +\ynorm ,
\end{equation}

instead of $T_0$. Multiplication by $2\xi$ may also be written as $\fa_\xinorm+\fa_\xinorm^*$ with $\fa_\xinorm = \partial_{\xinorm}+\xinorm$.

§.§ Functional setting.

To study the operator $H$, we use the standard basis of Hermite functions on $L^2(\Rm)$ given by

\begin{equation}
h_0(\zeta)=\pi^{-\frac14}e^{-\frac12 \zeta^2}, \qquad h_n(\zeta) = \dfrac{1}{2^{\frac n2} \sqrt{n!}} (\fa_\zeta^*)^n h_0(\zeta),\ n\geq1.
\end{equation}

They satisfy the relations

\begin{equation}\label{eq:aop}
\fa_\zeta h_0=0;\qquad \fa_\zeta h_n=\sqrt{2n} h_{n-1},\  n\geq1; \qquad \fa_\zeta^* h_n=\sqrt{2n+2} h_{n+1}, \ n\geq0.
\end{equation}

With the decomposition $L^2(\Rm)\ni \psi(\zeta)=\sum_{n\geq0} \psi_nh_n(\zeta)$, we may then define

\begin{align}\label{eq:inva}
\fa_\zeta^{-1} \psi(\zeta) = \sum_{n\geq1} \frac{1}{\sqrt{2n}} \psi_{n-1} h_{n}(\zeta), \quad & \psi \in L^2(\Rm) \\
(\fa_\zeta^*)^{-1} \psi(\zeta) = \sum_{n\geq0} \frac{1}{\sqrt{2n+2}} \psi_{n+1} h_{n}(\zeta), \quad & \psi \in ({\rm Ker}\,  \fa_\zeta)^\perp = \overline{{\rm Ran}\, \fa_\zeta^*} = \{\psi\in L^2(\Rm); \ \psi_0=0\}.
\end{align}

A natural scale of Hilbert spaces for $p\in\Nm$ associated to this decomposition is $\mS_p(\Rm)=\{\psi \in L^2(\Rm); (\fa_\zeta^*)^p\psi \in  L^2(\Rm)\}$ endowed with the norm

\begin{align}\label{eq:normp1}
\|\psi\|^2_p = \|(\fa_\zeta^*)^p \psi \|^2_0 \ \cong\  \|\psi\|^2_0 + \|D_\zeta^p \psi\|_0^2 + \|\zeta^p \psi\|^2_0,
\end{align}

where $\|\cdot\|_0$ is the usual $L^2(\Rm)-$norm and where $a\cong b$ when there exists $C_p>0$ such that $C_p^{-1}a\leq b \leq C_p a$. That the two above expressions for the norm are equivalent may be obtained by induction from the classical result

\begin{align}
\|\fa_\zeta^* \psi \|^2_0 = \|D_\zeta \psi\|^2_0 + \|\zeta \psi\|^2_0 + \|\psi\|^2_0
\end{align}

for $\psi\in \mS_1(\Rm)$. Note that $\mS_0(\Rm)\equiv L^2(\Rm)$.

With the decomposition  $L^2(\Rm^2)\ni \psi(\xi,\zeta)=\sum_{m,n\geq0}\psi_{mn}h_m(\xi)h_n(\zeta)$ in two-dimensional spaces, we similarly define $\mS_p(\Rm^2)= \{ \psi\in L^2(\Rm^2); (\fa_\xi^*)^p\psi\in L^2(\Rm^2) \mbox{ and } (\fa_\zeta^*)^p\psi\in L^2(\Rm^2)\}$ endowed with the norms

\begin{align}\label{eq:normp2}
\|\psi\|^2_p = \dsum_{j=0}^p \|(\fa_\xi^*)^j(\fa_\zeta^*)^{p-j} \psi \|^2_0 \ \cong\   \|(1+|\xi|^p+|\zeta|^p) \psi\|^2_0+ \|D_\xi^p \psi\|_0^2 + \|D_\zeta^p \psi\|_0^2.
\end{align}

For vector-valued functions, we also define the spaces $\mS_p(\Rm^d,\Cm^q)$ component-wise. Below we simply write $\mS_p$ when the domain $\Rm^d$ and range $\Cm^q$ are clear and $\|\cdot\|_p$ the associated norm.

§.§ Inversion of model operator.

As the lemma below demonstrates, the kernel $N\subset L^2(\Rm^2,\Cm^2)$ of $H$ and its orthogonal complement are given by

\begin{align}
\label{eq:N}
N=&\left\{\psi(\xinorm,\ynorm)=f(\xinorm)h_0(\ynorm) \matrice{ 0 \\ 1},\ f\in L^2(\Rm)\right\},\\  N^\perp =& \left\{\psi \in L^2(\Rm^2,\Cm^2),\ \big(\psi_2(\xi,\cdot),h_0(\cdot)\big)_2=0,\  \forall \xi\in\Rm \right\},
\end{align}

where $(\cdot,\cdot)_2$ is the standard $L^2(\Rm)-$inner product in the second variable $\zeta$ and $\psi=(\psi_1,\psi_2)^t$.

For $p\in \Nm$, we define the bounded linear operator $H^{-1}$ from $N^\perp\cap\mS_{p+1}$ to $N^\perp\cap \mS_{p}$ as

\begin{align}
H^{-1} = \matrice{ 0  & (\fa_\zeta^*)^{-1} \\ \fa_\zeta^{-1} & - \fa_\zeta^{-1} 2\xi (\fa_\zeta^*)^{-1} } .
\end{align}

All solutions $\psi\in \mS_{p}$ of the equation $H\psi=g$ for $g \in N^\perp\cap\mS_{p+1}$ are of the form $\psi=H^{-1}g+\psi_0$ with $\psi_0$ arbitrary in $N\cap\mS_{p}$. In particular ${\rm Ker}\, H=N$ on $\mS_0$.

The system $H\psi=g$ is equivalent to $\fa^*_\zeta \psi_1=g_2$ and $\fa_{\zeta}\psi_2 + 2\xi \psi_1 = g_1$. Using (<ref>), we find for $(g_2,h_0(\zeta))_2=0$, i.e., $g\in N^\perp$, that $\psi_1=(\fa^*_\zeta)^{-1}g_2$ and then $\psi_2=\fa_\zeta^{-1} (g_1-2\xi\psi_1)$ and hence the above result. We then note that $\fa_\zeta^{-1}$ has range in $N^\perp$.

That $\psi\in\mS_{p}$ when $g\in\mS_{p+1}$ follows directly from the fact that the following operators are bounded: $\fa_\zeta^{-1}:\mS_{p+1}\to\mS_{p+1}\cap N^\perp$, $(\fa_\zeta^*)^{-1}:\mS_{p+1}\cap N^\perp\to\mS_{p+1}$, as well as $\xi:\mS_{p+1}\to \mS_p$.

Any function $\psi\in\mS_{p}$ solution of $H\psi=g$ may be decomposed as $\psi_0+\psi_1$ with $\psi_0\in N$ and $\psi_1\in N^\perp$. Since $H\psi_0=0$ and hence $\psi_1=H^{-1}g$, $\psi_0$ is arbitrary in ${\rm Ker}\,H=N$.

§.§ Microscopic balance.

Our geometric assumptions impose that $C^{-1}\geq\rho_t\geq C>0$ and that $0<c_0\leq c_t<1$ for constants $C$ and $c_0$ independent of $t\geq0$.

For a fixed time $t$, we now solve the equation $T_0 a=b$.
We introduce

\begin{equation}\label{eq:phi}
\phi_t(\zetatwo) := \left(\frac {\rho_t}{4\pi}\right)^{\frac14} e^{-\frac{\rho_t}{2} \zetatwo^2} \matrice{ 1\\-1 },\qquad N_t = {\rm Ker}\ T_0 = \big\{ f(\xione) \phi_t(\zetatwo); \ f\in \mS_0(\Rm;\Cm) \big\}.
\end{equation}

The normalization implies that $\|\phi_t\|_{L^2(\Rm,\Cm^2)}=1$. That $N_t$ is the kernel of $T_0$ comes from Lemma <ref> and the change of variables given in the last section.
We still denote by $\mS_p(\Rm^2,\Cm^2)$ the spaces of functions in the $(\xione,\zetatwo)$ variables, which are equivalent to the corresponding spaces in the variables $(\rho_t^{-1}c_t\xione,\rho_t^{1/2}\zetatwo)$ by assumption on $(c_t,\rho_t)$.

The leading order equation is then solved as follows.

Let $b\in N_t^\perp \cap \mS_{p+1}$. The equation $T_0a=b$ admits a unique solution $T_0^{-1}b:=a\in N_t^\perp \cap \mS_{p}$ with inverse operator given explicitly by

\begin{align}
T_0^{-1} = Q \matrice{ 0 & (\rho_t\zeta-\partial_\zeta)^{-1} \\ (\rho_t\zeta+\partial_\zeta)^{-1} & - 2\xi c_t (\rho_t\zeta+\partial_\zeta)^{-1}(\rho_t\zeta-\partial_\zeta)^{-1} } Q ,\quad Q=\frac{1}{\sqrt 2} (\sigma_1+\sigma_3).
\end{align}

All solutions of that equation in $\mS_{p}$ are of the form $a=T_0^{-1}b+f(\xione)\phi(\zetatwo)$ for arbitrary $f(\xione)\in \mS_p(\Rm,\Cm)$.

This is a corollary of Lemma <ref> and the invertible transform from $T_0$ to $H$.
The operator $T_0$ is thus invertible on the orthogonal complement of $\phi_t(\zetatwo)$ with $T_0a=b$ solvable if and only if $(b,\phi_t)_2=0$ in which case all solutions are given by $a=T_0^{-1}b+f(\xione)\phi_t$ for $f$ arbitrary in $\mS_p$.

The solution to the leading equation $T_0a_0=0$ is therefore given by

\begin{equation}\label{eq:a0}
a_0(t,\Xitot) = f_0(t,\xione) \phi_t(\zetatwo)
\end{equation}

with $f_0(t,\xione)$ arbitrary at this level.

§.§ Functional setting.

To study the operator $H$, we use the basis of Hermite functions $h_{mn}(\xi,\zeta) = h_m(\xi) h_n(\zeta)$ on $L^2(\R^2)$, where $h_n$ denote the one-dimensional Hermite functions:

\begin{equation}
h_0(\zeta)=\pi^{-1/4}e^{-\frac12 \zeta^2}, \qquad h_n(\zeta) = \dfrac{(\fa_\zeta^*)^n}{2^{n/2} \sqrt{n!}} h_0(\zeta).
\end{equation}

These satisfy the relations

\begin{equation}\label{eq:aop}
\fa_\zeta h_0=0;\qquad \fa_\zeta h_n=\sqrt{2n} h_{n-1},\  n\geq1; \qquad \fa_\zeta^* h_n=\sqrt{2n+2} h_{n+1}, \ n\geq0.
\end{equation}

We decompose below functions $\psi\in L^2(\Rm^2,\C^2)$ as

\begin{equation}\label{eq:decpsimn}
\psi(\xinorm,\ynorm)=
\dsum_{n\geq0} \psi_n(\xi) h_n(\zeta) =  \sum_{m,n\geq0}\psi_{mn}h_{mn}(\xinorm,\ynorm), \qquad \psi_{mn} \in \C^2.
\end{equation}

We introduce the Hilbert spaces $\mS_p(\Rm^2,\C^2)\subset L^2(\Rm^2,\C^q)$, $q=1,2$ equipped with the  norm

\begin{equation}\label{eq:normp}
\|\psi\|^2_p = \dsum_{j=1}^2 \dsum_{m,n\geq0} \aver{m,n}^{p} |\psi_{mnj}|^2, \quad \aver{m,n}=\sqrt{1+m^2+n^2}.
\end{equation}

We note that the operators $D_\xi, D_\zeta$, as well as multiplication by $\xi$ and $\zeta$, are bounded from $\SSS_p(\R^2, \C^2)$. We define likewise spaces $\SSS_p(\R^d,\C^q)$; below we simply write $\mS_p$ when the domain $\Rm^d$ and range $\Cm^q$ are clear.

§.§ Inversion of model operator.

As the lemma below demonstrates, the kernel $N\subset L^2(\Rm^2,\Cm^2)$ of $H$ and its orthogonal complement are given by

\begin{equation}
\label{eq:N}
N=\left\{\psi(\xinorm,\ynorm)=f(\xinorm)h_0(\ynorm) \matrice{ 0 \\ 1},\ f\in L^2(\Rm)\right\},\ \  N^\perp=\left\{\psi \in L^2(\Rm^2,\Cm^2),\ \psi_{02}(\xi)=0 \ \forall \xi\in\Rm\right\}.
\end{equation}

For $p,q\geq0$, we define $H^{-1}$ as the bounded linear operator mapping $g \in N^\perp\cap\mS_{p+1}$ to $\psi\in N^\perp\cap \mS_{p}$ given in the partial Hermite basis by

\[
\psi_{n1}(\xi)=\frac{1}{\sqrt{2n+2}} g_{n+1,2}(\xi),\ \ n\geq0,\qquad \psi_{n2}(\xi) = \frac{1}{\sqrt{2n}} g_{n-1,1}(\xi) - \frac \xi n g_{n2}(\xi), \ \ n\geq1.
\]

and $\psi_{02}(\xi)=0$.

All solutions $\psi\in \mS_{p}$ of the equation $H\psi=g$ for $g \in N^\perp\cap\mS_{p+1}$ are of the form $\psi=H^{-1}g+\psi_0$ with $\psi_0$ arbitrary in $N\cap\mS_{p}$. In particular ${\rm Ker}\, H=N$ on $\mS_0$.

Decomposing the spinors $\psi$ and $g$ in the partial basis of Hermite functions in $\zeta$ and using the expressions (<ref>), the equation $H\psi=g$ is equivalent to

\begin{align}
2\xi \psi_{n1}(\xi) + \sqrt{2n+2} \psi_{n_1,2}(\xi) &=g_{n1}(\xi) , \\
\sqrt{2n}\psi_{n-1,1}(\xi) &= g_{n2}(\xi),
\end{align}

for $n\geq0$ and $\xi\in\Rm$ with the convention that $\psi_{(-1)j}=0$ for $j=1,2$ (since $\fa_\zeta h_0(\zeta)=0$). The second set of equations implies the compatibility condition $g_{02}(\xi)=0$ for all $\xi\in\Rm$, i.e., $g\in N^\perp$.

This is a triangular system of equations for $\psi_{nj}(\xi)$ whose solution is given in the lemma. Note that $\psi_{02}(\xi)$ is not constrained by the above system and is therefore set to $0$. This defines $H^{-1}$ with range in $N^\perp$. We may equivalently define

\begin{align}
\psi_1=(\fa_\zeta^*)^{-1} g_2,\qquad \psi_2= \fa_\zeta^{-1} (g_1-2\xi (\fa_\zeta^*)^{-1} g_2).
\end{align}

It remains to verify that $\psi\in\mS_{p}$ when $g\in\mS_{p+1}$. The map $g_{n}(\xi)\to g_{n-1}(\xi)$ is bounded from $\mS_{p}$ to itself  so that both $(\fa_\zeta^*)^{-1}$ and $\fa_\zeta^{-1}$ are bounded on $N^\perp\cap \mS_p(\Rm^2,\Cm^2)$. Multiplication by $\xi$ is bounded from $N^\perp\cap\mS_{p+1}$ to $N^\perp\cap \mS_{p}$, and hence so is $H^{-1}$.

Any function $\psi\in\mS_{p}$ solution of $H\psi=g$ may be decomposed as $\psi_0+\psi_1$ with $\psi_0\in N$ and $\psi_1\in N^\perp$. Since $H\psi_0=0$ and hence $\psi_1=H^{-1}g$, $\psi_0$ is arbitrary in ${\rm Ker}\,H=N$.

§.§ Microscopic balance.

Our geometric assumptions impose that $C^{-1}\geq\rho_t\geq C>0$ and that $0<c_0\leq c_t<1$ for constants $C$ and $c_0$ independent of $t\geq0$.

For a fixed time $t$, we now solve the equation $T_0 a=b$.
We introduce

\begin{equation}\label{eq:phi}
\phi_t(\zetatwo) := \left(\frac {\rho_t}{4\pi}\right)^{1/4} e^{-\frac{\rho_t}{2} \zetatwo^2} \matrice{ 1\\-1 },\qquad N_t = {\rm Ker}\ T_0 = \big\{ f(\xione) \phi_t(\zetatwo); \ f\in \mS_0(\Rm;\Cm) \big\}.
\end{equation}

The normalization implies that $\|\phi_t\|_{L^2(\Rm,\Cm^2)}=1$. That $N_t$ is the kernel of $T_0$ comes from Lemma <ref> and the change of variables given in the last section.

We still denote by $\mS_p(\Rm^2,\Cm^2)$ the spaces of functions in the $(\xione,\zetatwo)$ variables, which are equivalent to the corresponding spaces in the variables $(\rho_t^{-1}c_t\xione,\rho_t^{1/2}\zetatwo)$

by assumption on $(c_t,\rho_t)$.

The leading order equation is then solved as follows.

Let $b\in N_t^\perp \cap \mS_{p+1}$. The equation $T_0a=b$ admits a unique solution $T_0^{-1}b:=a\in N_t^\perp \cap \mS_{p}$. The inverse operator $T_0^{-1}$ is bounded from $N_t^\perp \cap \mS_{p+1}$ to $N_t^\perp \cap \mS_{p}$. All solutions of that equation in $\mS_{p}$ are of the form $a=T_0^{-1}b+f(\xione)\phi(\zetatwo)$ for arbitrary $f(\xione)\in \mS_p(\Rm,\Cm)$.

This is a corollary of the previous result knowing that the time dependent transform from $T_0$ to $H$ preserves the above functional spaces.
The operator $T_0$ is thus invertible on the orthogonal complement of $\phi_t(\zetatwo)$ with $T_0a=b$ solvable if and only if $(b,\phi_t)_2=0$ (standard inner product in $L^2(\Rm,\Cm^2)$ in the $\zetatwo$ variable) in which case all solutions are given by $a=T_0^{-1}b+f(\xione)\phi_t$ for $f$ arbitrary in $\mS_p$.

The solution to the leading equation $T_0a_0=0$ is therefore given by

\begin{equation}\label{eq:a0}
a_0(t,\Xitot) = f_0(t,\xione) \phi_t(\zetatwo)
\end{equation}

with $f_0(t,\xione)$ arbitrary at this level.

§ TRANSPORT EQUATION

Lemma <ref> states that the subleading equation in (<ref>), $T_0a_1=-T_1a_0$, admits a solution if and only if $T_1a_0(t,\cdot) \in N_t^\perp$ for every $t$. From the expressions (<ref>) of $a_0$ and (<ref>) of $N_t$ and $\phi_t(\zeta)$, we deduce that $T_0a_1=-T_1a_0$ is solvable if and only if $\mT f_0 = 0$, where

\begin{equation}\label{eq:mT}
\mT f(t,\xione):=  \int_{\mathbb{R}} \phi_t(\zetatwo)\cdot T_1 [f(t,\xione) \phi_t(\zetatwo)] \ d\zetatwo.%, \qquad \phi(\zetatwo) = \left(\frac {\rho_t}{4\pi}\right)^{1/4} e^{-\frac{\rho_t}{2} \zetatwo^2} \matrice{ 1\\-1 }.
\end{equation}

In this section we provide an explicit expression for $\mT$ and for the solutions to $\mT f_0 = 0$ and more generally $\mT f_0=g$. We also provide a functional setting to analyze the map $g\to f_0$ when the dispersion is strongest.

§.§ Derivation of the transport operator

We have the identity

\begin{equation}\label{eq:transport}
\mT = D_t  - \dfrac{k_t}{2\rho_t}  - j_t \gamma_t \dfrac{\xione D_\xi +  D_\xi \xione}{2} -  \left(\dot{\te_t} \gamma_t + k_t\gamma_t^2  \right)\xione^2.
\end{equation}

We first observe that

\begin{equation}
\left \langle \matrice{1 \\ -1}, \sigma_j \matrice{1 \\ -1} \right \rangle = -2 \delta_{1j}.
\end{equation}

Therefore, the terms in $T_1$ that are carried by $\sigma_2, \sigma_3$ do not contribute to $\mT$. Likewise, we observe that

\begin{equation}
\int_{\mathbb{R}} \phi(\zetatwo) \cdot \zetatwo \phi(\zetatwo) d\zetatwo = \int_{\mathbb{R}} \phi(\zetatwo) \cdot D_\zeta \phi(\zetatwo) d\zetatwo = 0.
\end{equation}

Hence, for the purpose of computing $\mT$, we can ignore terms in $T_1$ carried by $\sigma_2$ or $\sigma_3$, and terms linear in $(\zetatwo,D_\zeta)$. In other words, we can replace $T_1$ in (<ref>) by the right hand side of (<ref>):

\begin{equation}
D_t + \dot{\te_t} \gamma_t \big(D_\zeta^2-\xione^2\big)  + \left( \dfrac{\dot{\te_t}s_t}{2}  +  j_t \gamma_t (D_\xi \xione - D_\zeta \zetatwo) + k_t \big(\zetatwo^2 + \gamma_t^2 \xione^2 \big)\right) \sigma_1.
\end{equation}

We note moreover that

\begin{equation}\label{eq-0x}
D_\zeta \zetatwo\phi =  i(\rho_t \zetatwo^2-1) \phi, \quad D_\zeta^2 \phi(\zetatwo) = (\rho_t - \rho_t^2 \zetatwo^2) \phi(\zetatwo), \ \ \ \ D_t \phi(\zetatwo) = i \dfrac{\dot{\rho_t}}{2} \zetatwo^2 \phi(\zetatwo) - i \dfrac{\dot{\rho_t}}{4\rho_t} \phi(\zetatwo).
\end{equation}

Therefore, we deduce that

\begin{equation}
\begin{array}{l}
\mT =\Big(\dfrac{\rho_t}{\pi}\Big)^{\frac12} \dint_{\mathbb{R}}    e^{-\rho_t \zetatwo^2}
\Big( D_t + i \dfrac{\dot{\rho_t}}{2} \zetatwo^2 -i\frac1{4} \frac{\dot\rho_t}{\rho_t}
+ \dot{\te_t} \gamma_t \big(\rho_t   - \rho_t^2 \zetatwo^2-\xione^2\big) \\[3mm]  \qquad - \Big( \dfrac{\dot{\te_t}s_t}{2}  +  j_t \gamma_t (D_\xi \xione - i \rho_t \zetatwo^2 + i) + k_t \big(\zetatwo^2 + \gamma_t^2 \xione^2 \big)\Big) \Big)
\end{array}
\end{equation}

Furthermore, we have

\begin{equation}
\left(\frac {\rho_t}{\pi}\right)^{\frac12} \int_{\mathbb{R}}
e^{-\rho_t \zetatwo^2}  d\zetatwo = 1, \ \ \ \ \left(\frac {\rho_t}{\pi}\right)^{\frac12} \int_{\mathbb{R}} \zetatwo^2
e^{-\rho_t \zetatwo^2}  d\zetatwo = \dfrac{1}{2\rho_t}.
\end{equation}

Hence, after performing the integration and realizing that the coefficients involving $\dot\rho_t$ cancel out, we obtain the formula

\begin{align}
\label{eq-0r} \mT & =
D_t %+ i \dfrac{\dot{\rho_t}}{4\rho_t}
+ \dot{\te_t} \gamma_t \left(\dfrac{\rho_t}{2}-\xione^2\right) - \dfrac{\dot{\te_t}s_t}{2}  -  j_t \gamma_t \left(D_\xi \xione + \dfrac{i}{2}\right) - k_t \left(\dfrac{1}{2\rho_t} + \gamma_t^2 \xione^2 \right)
\\
\label{eq-0s}
& = D_t
%+ i \dfrac{\dot{\rho_t}}{4\rho_t}
- \dot{\te_t} \gamma_t \xione^2 -  j_t \gamma_t \dfrac{\xione D_\xi +  D_\xi \xione}{2} - k_t \left(\dfrac{1}{2\rho_t} + \gamma_t^2 \xione^2 \right)
\\
& = D_t
%+ i \dfrac{\dot{\rho_t}}{4\rho_t}
- \dfrac{k_t}{2\rho_t}  -  j_t \gamma_t \dfrac{\xione D_\xi +  D_\xi \xione}{2} -  \left(\dot{\te_t} \gamma_t + k_t\gamma_t^2  \right)\xione^2,
\end{align}

where in the second line we used $\gamma_t \rho_t = s_t$.

I worked a bit on the solution of the transport equation. There are a few  advantages of writing it as follows:

* The compression factor is in fact just $c_t = \cos\varphi_t$. This rules out exponential compression / stretching. (before we had $e^\lambda_t$ with $\lambda_t = -\int_0^t j_s \gamma_s ds$, so we could only say $e^{-ct} \leq e^{\lambda_t} \leq e^{Ct}$).

* The whole solution is dilated through $c_t$. Benefit is two-fold: simpler formula for the coefficient $\nu_t$ in $\xi_1^2$ in the phase (before $\nu_t$ solved an ODE); simpler scheme when we will have to later reverse the Fourier-like transform.

The compression factor $\mu_t$ disappears eventually by an opposite compression in $z_2$ once we are back in spatial variables. Where is this written? It is very surprising to me: we should still see compression/stretching, just not of infinite stength. For instance see Figure 1 in magnalexis (reproduced here as Figure <ref>, but the stretching is a bit harder to see).
That's because the stretching comes from somewhere else; see numerical section. The thing to look at, if correct of course, is (6). For gaussian profiles, this gives (7). The latter should explain all the effects we see in the numerics. In particular, there is a factor $Q_t$. It is the one responsible for amplitude changes via $\nu_t$. We should also observe an amplifying factor coming from $\rho_t$ but I do not think we have numerics for that.
These are all bounded compression (in fact I would argue that they are rather dispersion factors). But you have dilated variables using (your) $\mu_t$. How do you see that $\mu_t$ is finite? I agree that was an assumption; it's nice that you relate it to $c_t$; I totally agree using (58). In some ways, that also gives an expression for $j_t$ in terms of $\gamma_t$ and $c_t$.
Yes, the one from $\rho_t$ should be finite. The only one that is increasing with time come from (my) $\nu_t$ and from the coefficient $k_t$. Incidentally, it is also the one that is most painful to write. $j_t$ has no such role.

Ok, well why not use (<ref>) as the solution of the transport equation then, since it is clear from it (and without asumptions) that there is no exponential compression?

There is infinite compression coming from $Q_t$, that's what I'm saying. It's not exponential but power law as observed by Simon; in fact $T^{-\frac12}$ it seems.  I would call that effect dispersion.

The whole wavepacket has the same macroscopic speed  and so we do not have microscopic dispersion. I agree this effect is mesoscopic dispersion. This has to be verified but I was surprised to see that $\mu_t$ has no rescaling effect; just some kind of rotation in the $z$ variables. The main compression/dispersion comes from $Q_t$ and $\nu_t$.

Look at the solution of the Schrodinger equation:

\begin{equation}
\dfrac{1}{\sqrt{t}} \int_\R e^{it \xi^2} \hat{f} (\xi) d\xi.
\end{equation}

It disperses, right? Well, that's the same as (<ref>) with no phase and $\nu_t = t$. ($\nu_t = t$ happens e.g. for the circle and $B$ constant).

OK. So, $\mu_2 = \nu_t$ increases linearly in time I agree. This is what causes the dispersion by $Q_t$. Bottom line, the important effect is hidden in $Q_t$. I could only make sense of it for a gaussian profile in order to get an explicit expression. I did not see the effects of $\rho_t$ or (my) $\mu_t$ in the simulations. ?? Look at the circle: there is clearly dispersion; the dispersion is maximal for $B=1/\sqrt{2}$, also predicted by $\mu_t$. I'm not sure this is predicted in the $z$ variables; only in the $(\xi_1,\zetatwo)$ ones.

OK. $\mu_t$ generates some form of compression. But again, what we gain in $z_1$ we lose in $z_2$. So amplitude effects (which are not dispersive effects) are not to be found in $\mu_t$ imo. Do you see maximal compression for a given value of $B_t$ in (6) or (7)? These are the expressions for the spatial profiles. We do not see them at the level of the $\xi/\zeta$ variables.  Ok, I'll check

Ok, I'll continue working details out of <ref>.

Here is why I hesitate with dispersion. The phase factor $e^{i\nu_t\xi_1^2}$ is indeed generating dispersion when $\nu_t\propto t$, say. But the wavepacket at the microscale still propagates with speed $\dot y_t$. We do not have any microscopic- (Dirac-) level dispersion. However, we indeed observe mesoscopic- $\sqrt\eps-$ level dispersion.

§.§ Solving the transport equation

To produce an explicit solution of the transport equation $\mT f = 0$ we define $(\lambda_t,\mu_t,\nu_t)$ such that

\begin{equation}\label{eq-1a}
\lambda_t =  \int_0^t \dfrac{k_s}{2\rho_s}ds, \qquad \nu_t =  2\int_0^t \dfrac{c_0^2}{c_s^2} \left( \dot{\te_s} \gamma_s + k_s \gamma_s^2\right) ds, \qquad e^{\mu_t} = \frac{c_t}{c_0}.
\end{equation}

The solution to $\mT f = 0$ is given by

\begin{equation}\label{eq-1n}
%f(t,\xione) =  e^{\frac{\mu_t}2} \exp\left( i \lambda_t + i\dfrac{\nu_t}{2} \left( e^{\mu_t}\xione \right)^2 \right) f\left(0, e^{\mu_t}\xione \right).
f(t,\xione) =  \exp\left( i \lambda_t + i\dfrac{\nu_t}{2} \left(e^{\mu_t}\xi \right)^2 \right) e^{\frac{\mu_t}2} f\left(0, e^{\mu_t} \xi \right).
\end{equation}

1. We recall that the self-adjoint operator $\frac{1}{2} (\xione D_\xi +  D_\xi \xione)$ generates the semigroup of ($L^2-$unitary) dilations; that is, for $F$ independent of $t$:

\begin{equation}\label{eq-1k}
\left(D_\mu - \dfrac{\xione D_\xi + D_\xi \xione}{2}\right) U_\mu F = 0, \quad U_\mu F(\xione) = e^{\frac\mu 2} F\big(e^\mu \xione\big).
\end{equation}

We note that $e^{\mu_0} =1$; moreover, thanks to (<ref>), we have $\dot{\mu_t} = \p_t \ln(c_t) = j_t \gamma_t$. Therefore, by the chain rule and (<ref>), we have

\begin{equation}
U_{\mu_t}^{-1} \left(D_t - j_t \gamma_t \dfrac{\xione D_\xi + D_\xi \xione}{2}\right) U_{\mu_t} = D_t.
\end{equation}

2. We deduce, using $U_{\mu_t}^{-1} \xione^2 U_{\mu_t}=e^{-2\mu_t}\xione^2$, that

\begin{align}
U_{\mu_t}^{-1} \mT U_{\mu_t} & =  D_t - \dfrac{k_t}{2\rho_t} - e^{-2\mu_t} \left( \dot{\te_t} \gamma_t + k_t \gamma_t^2\right) \xione^2
%& = D_t + \dot{\lambda_t} + \dfrac{1}{2}\dot{\nu_t} \xione^2  = e^{i\lambda_t + i\nu_t \frac{\xione^2}{2}} D_t e^{-i\lambda_t -i \nu_t \frac{\xione^2}{2}},
= D_t - \dot{\lambda_t} - \dfrac{1}{2}\dot{\nu_t} \xione^2  = e^{i\lambda_t + i\nu_t \frac{\xione^2}{2}} D_t e^{-i\lambda_t - i \nu_t \frac{\xione^2}{2}},
\end{align}

where we used the relations (<ref>) for $\lambda_t$ and $\nu_t$. For the formula (<ref>) for $U_\mu$, we deduce that the solution to $\mT f=0=\mT U_{\mu_t}e^{i\lambda_t + i \frac{\nu_t}{2} \xione^2} g$ is $D_tg=0$ and hence

\begin{equation}
f(t,\cdot) = U_{\mu_t} e^{i\lambda_t + i \frac{\nu_t}{2} \xione^2} f(0,\cdot).
\end{equation}

Since $e^{\mu_t} = c_t/c_0$, we conclude that

\begin{equation}\label{eq-3r}
f(t,\xione) = e^{\frac{\mu_t}2} \exp\left( i \lambda_t + i\frac{\nu_t}{2} \left( e^{\mu_t}\xi \right)^2 \right) f\left(0, e^{\mu_t}\xi \right).
\end{equation}

This completes the proof.

§.§ Dispersive estimate

In this section, we study the $L^\infty$-decay of the leading order solution of (<ref>), $\psi_0 = \fU_t a_0$, where $a_0$ takes the form prescribed by (<ref>) and (<ref>):

\begin{equation} \label{eq-00m}
a_0 = e^{i\lambda_t} \tilde a_0 \matrice{1\\-1}, \quad  \tilde a_0(t,\xi,\zeta)= \left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}  e^{\frac{\mu_t}2}  \exp\left(i\dfrac{\nu_t}{2} \left( e^{\mu_t}\xione \right)^2 -\frac{\rho_t}{2} \zetatwo^2 \right)  \hf\left(e^{\mu_t}\xione \right),
\end{equation}

where $\hf(\xi)$ denotes the Fourier transform of a function $f(z_1)$.

(i) With $\psi_0 = \fU_t a_0$ and $a_0$ given by (<ref>), we have whenever $B_t \neq 0$ or $\nu_t \neq 0$:

\begin{equation}\label{eq-3i}
\psi_0(t,z) = e^{i\lambda_t} \bU_t \big( G_t *_1 f \big)
%\left(e^{\mu_t}z \right)
(e^{-\mu_t} z_1, e^{\mu_t} z_2)
\matrice{1\\-1},
\end{equation}

where $*_1$ denotes convolution with respect to the first variable and

\begin{equation}\label{eq-3h}
G_t(e^{-\mu_t} z_1,e^{\mu_t}z_2) =
\left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}
\Big( \frac{e^{\mu_t}}{Q_t} \Big)^{\frac12}
e^{-\frac12 \rho_tz_2^2}  e^{-\frac12 Q_t^{-1}(z_1+i s_t z_2)^2}, \quad Q_t = s_t\gamma_t -i e^{2\mu_t}\nu_t.
\end{equation}

We recall that $s_t=\gamma_t\rho_t$.

(ii) In particular, there exists $C > 0$ such that as long as $\nu_t \neq 0$:

\begin{equation}\label{eq-3k}
\sup_{z \in \R^2} \big| \psi_0(t,z) \big| \leq   C  \min \Big ( \dfrac{\| f \|_{L^1}}{|\nu_t|^{1/2}}  \,,\,\|\fhat\|_{L_1} \Big).
\end{equation}

We recall that $\fU_t = \bU_t \mV_t$ and we first write a formula for $\mV_t \tilde a_0$. We have:

\begin{equation}
\mV_t \tilde a_0(z) =  \left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}e^{\frac{\mu_t}2} \int_\R e^{iz_1\xione}  \exp\left( i\dfrac{\nu_t}{2} \left( e^{\mu_t}\xione \right)^2 -\frac{\rho_t}{2} (z_2+\gamma_t \xi)^2 \right) \fhat\left(e^{\mu_t}\xione \right) \dfrac{d\xione}{\sqrt{2\pi}}.
\end{equation}

\begin{align}
\mV_t \tilde a_0(e^{\mu_t}z_1,\frac{z_2}{e^{\mu_t}}) &= \left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}
\int_\R e^{i e^{\mu_t} z_1\xione}
\exp\left( i\dfrac{\nu_t}{2} \left( e^{\mu_t}\xione \right)^2 -\frac{\rho_te^{-2\mu_t}}{2} (z_2+\gamma_t e^{\mu_t}\xi)^2 \right)
\fhat\left(e^{\mu_t}\xione \right) \dfrac{e^{\frac{\mu_t}2} d\xione}{\sqrt{2\pi}}
\\ \label{eq:secondbound}
&=\left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}
\int_\R e^{iz_1\xione}  \exp\left( i\dfrac{\nu_t}{2}
\xione^2 -\frac{\rho_te^{-2\mu_t}}{2} (z_2+\gamma_t \xi)^2 \right)
\fhat\left(\xione \right) \dfrac{e^{-\frac{\mu_t}2}d\xione}{\sqrt{2\pi}}.
\end{align}

The left-hand side involves the inverse Fourier transform in $z_1$ of a product and may therefore be written as the convolution

\begin{align}\label{eq:convta0}
\mV_t \tilde a_0(e^{\mu_t}z_1,z_2) = (f *_1 G_t)(z_1,e^{\mu_t}z_2)
\end{align}

where the Gaussian $G_t$ is given by

\begin{equation}\label{eq-3m}
G_t(z) = \left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}e^{\frac{-\mu_t}2} \int_\R e^{iz_1 \xi }  \exp\left( i\dfrac{\nu_t}{2} \xi^2 -\frac{\rho_te^{-2\mu_t}}{2} (z_2+\gamma_t \xi )^2 \right) \dfrac{d\xione}{\sqrt{2\pi}},
\end{equation}

so that

\begin{align}
G_t(z_1,e^{\mu_t}z_2) &= \left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}e^{\frac{-\mu_t}2} e^{-\frac12 \rho_tz_2^2} \dint_{\Rm} e^{-\frac12 \tilde Q_t\xi^2} e^{i\xi(z_1+i\rho_t\gamma_t e^{-\mu_t} z_2)} \dfrac{d\xione}{\sqrt{2\pi}},\qquad \tilde Q_t = \rho_t\gamma_t^2e^{-2\mu_t} -i\nu_t
\\ \label{eq:Qt}
&= \left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}e^{\frac{-\mu_t}2} e^{-\frac12 \rho_tz_2^2} \tilde Q_t^{-\frac12} e^{-\frac12 \tilde Q_t^{-1}(z_1+i\rho_t\gamma_t e^{-\mu_t} z_2)^2},
\end{align}

as the Fourier transform of a Gaussian function when $\tilde Q_t\not=0$.  Since $Q_t=e^{2\mu_t}\tilde Q_t$, we find (<ref>). Then (<ref>) follows from applying $\bU_t$ to $\mV_t a_0$ using (<ref>) and  (<ref>).

To prove the second part of the lemma, we first take $L^\infty$-norms on both sides of (<ref>) and we apply Young's $L^1-L^\infty$ convolution inequality  (in the variable $z_1$). This produces:

\begin{align}\label{eq-3q}
\sup_{z \in \R} \big| \mV_t \tilde a_0(z) \big| & = \sup_{z_2 \in \R} \sup_{z_1 \in \R}\big| f *_1 G_t(z) \big| \leq \sup_{z_2 \in \R} \| f \|_{L^1} \sup_{z_1 \in \R} \big| G_t(z) \big| = \| f \|_{L^1}  \sup_{z \in \R} \big| G_t(z) \big|.
\end{align}

The Gaussian $|G_t|$ attains its maximum at $z = 0$ because the real part of the quadratic form in $z$ in (<ref>) satisfies

\begin{align}
\Re \Big(\rho_tz_2^2 + \frac{1}{ Q_t}(z_1+is_t z_2)^2\Big) = \frac{\rho_t}{(s_t\gamma_t)^2+\tilde \nu_t^2}  (\gamma_t z_1 - \tilde \nu_t z_2)^2 \geq0, \qquad \tilde\nu_t=e^{2\mu_t}\nu_t,
\end{align}

as we verify by an elementary computation using $\rho_t\gamma_t=s_t$.

Using (<ref>), we find

\begin{equation}
\sup_{z \in \R} \big| \mV_t \tilde a_0(z) \big| \leq \| f \|_{L^1}  \big| G_t(0) \big| =  \left( \dfrac{\rho_t}{4\pi}\right)^{\frac14}\left( \dfrac{e^{\mu_t}}{|Q_t|} \right)^{\frac12} \| f \|_{L^1}.
\end{equation}

From (<ref>), we also obtain that $\sup_{z \in \R} \big| \mV_t \tilde a_0(z) \big|\leq C \|\hat f\|_{L^1}$.

To end up with (<ref>), we observe that $\bU_t$ preserves $L^\infty$-norms, that $|Q_t| \geq |\nu_t|$ and that $c_t$ is bounded above and below for all times $t$.

When $Q_t = 0$, we have $\nu_t = 0$, $c_t = 1$, $\rho_t = r_t$ and the formulas (<ref>) and (<ref>) remain valid. But instead of being a Gaussian, $G_t$ is now a multiple of the Dirac mass, and

\begin{equation}
G_t(z) =  c_0^{1/2} \left( \dfrac{r_t}{4\pi}\right)^{1/4}\exp\left( -\frac{r_t}{2} z_2^2 \right) \delta_0(z_1), %\qquad \mV_t a_0(t,z) =  c_0^{1/2} \left( \dfrac{r_t}{4\pi}\right)^{1/4} e^{-\frac{r_t}{2} z_2^2} f(z_1),
\qquad \psi_0(t,z) = c_0^{1/2} \left( \dfrac{r_t}{4\pi}\right)^{1/4} \bU_t e^{-\frac{r_t}{2} z_2^2} f(z_1).
\end{equation}

§.§ Stability estimates

For every $p$, there exist a constant $C$ independent of $T$ such that if $f, g \in \mC\mS_p$ are such that $\mT f = g$ and $f(0,\cdot) \equiv 0$ then

\begin{equation}
\big\| f \big\|_{\mC\mS_p}  \leq C T \lr{\nu}^p_T \big\| g\big\|_{\mC\mS_p}.
\end{equation}

1. Given $s \in \R$, let $G_s(t,\xi)$ be the solution to the equation $\mT G_s = 0$ starting from $G_s(s,\xi) = g(s,\xi)$. The analysis of <ref> – with initial time replaced by $s$ instead of $0$ – provides an explicit formula for $G_s$:

\begin{equation}
G_s(t,\xi) = \left( \dfrac{c_t}{c_s}\right)^{1/2} \exp\left( i (\lambda_t-\lambda_s) + i\dfrac{\nu_t-\nu_s}{2} \left( \dfrac{c_t \xione}{c_s} \right)^2 \right) g\left(s, \dfrac{c_t\xione}{c_s} \right).
\end{equation}

The norm on the space $\SSS_p$ involves at most $p$ derivatives with respect of $\xi$. Since $c_t$, $c_s$ are uniformly bounded with respect to time, we deduce that there exists a constant $C$ depending on $p$ only such that

\begin{equation}
\big\| G_s(t,\cdot) \big\|_{\mS_p} \leq C \lr{\nu_t-\nu_s}^p  \big\| g(s,\cdot) \big\|_{\mS_p}
\end{equation}

2. Duhamel's principle gives the solution to $\mT f = g$ with vanishing initial data:

\begin{equation}
f(t,\xi) = \int_0^t G_s(t,\xi) ds.
\end{equation}

From (<ref>) we obtain for every $t \in [0,T]$:

\begin{align}
\big\| f(t,\cdot) \big\|_{\mS_p} & \leq \int_0^t \big\| G_s(t,\cdot) \big\|_{\mS_p} ds \leq C \int_0^t\lr{\nu_t-\nu_s}^p  \big\| g(s,\cdot) \big\|_{\mS_p} ds
\\
\leq  C T \sup_{s,t \in [0,T]} \lr{\nu_t-\nu_s}^p \cdot \big\| g\big\|_{\mC\mS_p} \leq2^p C T \lr{\nu}^p_T \big\| g\big\|_{\mC\mS_p}.
\end{align}

It suffices to take the supremum over $t$ in the left hand side to complete the proof.

§.§ Construction of amplitudes

We construct the coefficients $a_j$ by recursion. Assume that $a_0, \dots a_{J-1}, b_J$ are constructed and let's construct $f_J$ and $b_{J+1}$. These must satisfy the equation:

\begin{equation}
T_0 b_{J+1} + T_1 (b_J + f_J \phi_0) + T_2 a_{J-1} + \dots + T_{J+1} a_0 = 0.
\end{equation}

For a solution to exist, the term $T_1 (b_J + f_J \phi_0) + T_2 a_{J-1} + \dots + T_J a_0$ must be orthogonal to the kernel of $T_0$, hence

\begin{equation}
\TT f_J = g_J, \qquad g_J := -\lr{T_1 b_J + T_2 a_{J-1} + \dots + T_J a_0, \phi_0}.
\end{equation}

The operators $T_j$ are bounded from $\SSS_p$ to $\SSS_{p-j-1}$ for all $p$. Therefore, we have the estimate

\begin{equation}
\| g_J \|_p \leq \| T_1 b_J \|_p + \| T_2 a_{J-1} \|_p + \dots +  \| T_{J+1} a_0 \|_p \lesssim \| b_J\|_{p-2} + \|a_{J-1}\|_{p-3} + \dots + \| a_0 \|_{p-J-2}.
\end{equation}

Previous derivation.

§.§ Dispersive estimate.

Let $\psi_0 = \fU_t a_0$, where $a_0(t,\Xitot)=f_0(t,\xione)\phi(\zetatwo)$ is the leading solution and assume that $f_0(0,\xione)=f(\xione)$ with inverse Fourier transform $\check f(z_1)\in L^1(\Rm)\cap L^\infty(\Rm)$. We wish to understand the spatial structure of $\psi_0(t,z)$ and in particular its dispersion along $\Gamma$.

Since $\mT f_0=0$, we obtain from Lemma <ref> that

\begin{equation}\label{eq-00mbis}
a_0(t,\xi,\zeta) = \left( \dfrac{\rho_t}{4\pi}\right)^{1/4}
\exp\left( i \lambda_t + i\frac{\tilde\nu_t}{2}\xione^2   - \dfrac{1}{2} \rho_t \zetatwo^2 \right) \mu_t^{\frac12} f(\mu_t\xione) \matrice{1 \\ -1},
\end{equation}

where we define $\mu_t=\frac{c_t}{c_0}$ We already had $e^\mu_t = c_t/c_0$ so maybe we should use something else? and $\tilde\nu_t=\nu_t\mu_t^2$ to simplify expressions. In the variables $z$, we thus find

\begin{equation}\label{eq:Nua0}
\mV_t  a_0(t,z) = %\dfrac{\rho_t^{1/4}\mu_t^{\frac12}}{2\pi^{\frac34}}
\rho_t^{1/4}\mu_t^{\frac12} \pi^{-1/4}
e^{i\lambda_t} \dint_{\Rm} \frac{1}{2\pi} e^{iz_1\xione}
\exp\left( i\frac{\tilde\nu_t}{2}\xione^2   - \dfrac{1}{2} \rho_t (z_2+\gamma_t\xione)^2 \right) f(\mu_t\xione) d\xione \matrice{1 \\ -1}.
\end{equation}

The inverse Fourier transform of the above product is the convolution of the inverse Fourier transform of the above exponential and of $\mu_t^{-1}\check f(\mu_t^{-1}z_1)$. We find

\begin{equation}\label{eq:Qtbis}
\exp\left( i\frac{\tilde\nu_t}{2}\xione^2   - \dfrac{1}{2} \rho_t (z_2+\gamma_t\xione)^2 \right) = e^{-\frac12 \rho_tz_2^2} e^{i(z_1+is_t z_2)\xione}e^{-\frac12 Q_t\xi^2},\qquad Q_t=s_t\gamma_t-i\tilde\nu_t,
\end{equation}

so that the inverse Fourier transform $\xione\to z_1$ of $\exp( i\frac{\tilde\nu_t}{2}\xione^2   - \frac{1}{2} \rho_t (z_2+\gamma_t\xione)^2)$ is

\[
G(t,z)=(2\pi Q_t)^{-\frac12} e^{-\frac12 \rho_t z_2^2} e^{-\frac12 Q_t^{-1} (z_1+is_tz_2)^2}.
\]

We use here the convention for Fourier transforms that $\mathcal{F}f(\xi)=\int_{\Rm}e^{-ix\xi}f(x)dx$, with thus a different normalization than the unitary operator $\mV_t^*$.

We verify when $Q_t\not=0$ that $Q_t^{-1}$ and $\rho_t-Q_t^{-1}s_t^2$ (since $\rho_t\gamma_t=s_t$) have non-negative real parts so that $|G|$ is maximal at $z=0$.

Define $\|X(t,\cdot)\|_\infty$ for a vector-valued function $X(t,z)$ as the supremum over $z$ at time $t$ of the Euclidean norm $\|X(t,z)\|$. Then we observe from the above bound and the expression of $\mV_t a_0$ as a convolution that

\begin{equation}\label{eq:dispersion}
\| \mV_t a_0(t,\cdot) \|_\infty \leq C \min(|Q_t|^{-\frac12} \|\check f\|_{L^1} ,\| f\|_1)
\end{equation}

for a constant $C$ independent of time. The second inequality is a direct consequence of (<ref>).

The terms $c_t$ (or $\mu_t$) and $\rho_t$ are responsible for a quantitative compression (or dilation) of the wave packet but do not qualitatively affect its structures. Similarly, the term $\lambda_t$ produces a phase shift that does not impact the wavepacket amplitude. The contribution associated to the term $\nu_t$ is more interesting and as mentioned in the introduction causes a potentially large dispersion over long times when $\nu_t$ grows linearly in time (i.e., when the interface remains curved as is the case for closed loop or when the magnetic field varies along the interface). Such an effect is very well captured by the numerical simulations in Figs. <ref>, <ref>, and <ref>.

The following lemma summarizes this dispersive effect.

Let $\psi_0=\fU_t a_0$. There exists a constant $C$ such that for all $t > 0$:

\begin{equation}
\sup_{z\in \R^2} \big\| \psi_0(t,z) \big\| \leq C |\nu_t|^{-\frac 12}  \|\check f\|_{L^1}.
\end{equation}

We observe that
\big\| \psi_0(t,\cdot) \big\|_\infty =\big\| \fU_t a_0(t,\cdot) \big\|_\infty =\big\| \mV_t a_0(t,\cdot) \big\|_\infty
by definition of the $\|\cdot\|_\infty$ norm.

We used above that
$\fU_t = \bU_t\mV_t$ and that rotations preserve amplitudes of functions. The proof then follows from (<ref>) and (<ref>) since $s_t\gamma_t$ is uniformly bounded.

Note that when $s_t=\gamma_t=0$ because $B_t=0$, then $\zetatwo=z_2$ and $a_0(t,\xione,z_2)$ in (<ref>) has a Gaussian profile in $z_2$ and experiences dispersion proportional to $|\nu_t|^{-\frac12}$ in the $z_1$ variable. When $B\equiv0$ so that $\nu_t=0$ and $\mu_t=1$, then $a_0(t,\xione,z_2)$ is the Fourier transform $z_1\to\xione$ of the leading profile considered in <cit.>.

We now use the definition (<ref>) of $\mV_t$:

\begin{align}
& \big\| \mV_t a_0(t,\cdot) \big\|  = \left\| \int_\R e^{iz_1\xione} a_0(\xione,z_2+\gamma_t\xione) d\xione \right\|_{L^\infty}
\\
& = \dfrac{\rho_t^{1/4}}{2^{1/2}\pi^{3/4}} e^{\frac{\mu_t}2} \left\| \int_\R e^{iz_1\xione} \exp\left( - \dfrac{1}{2} \rho_t (z_2+\gamma_t\xione)^2 \right) \exp\left(i\frac{\nu_t}{2} \left( e^{\mu_t}\xione \right)^2\right)  f\left(0, e^{\mu_t}\xione \right) d\xione \right\|_{L^\infty}.
\end{align}

We thus need to analyze the inverse Fourier transform $\xione\to z_1$ of a product of three functions, which is given by two successive convolutions in $z_1$.

Define the Gaussian

\begin{equation}
G_{t,z_2}(z_1) = \exp\left( - \dfrac{1}{2} \rho_t (z_2+\gamma_t\xione)^2 \right).
\end{equation}

When $\gamma_t=0$, the above term is uniformly bounded.
When $\gamma_t\not=0$, multiplication by $G_{t,z_2}$ involves a convolution by $\hat{G}_{t,z_2}$. Like $G_{t,z_2}$, this is a Gaussian with uniformly bounded coefficients since its dependence in $z_2$ is of the form $e^{i\gamma_t^{-1}z_2z_1}$. Hence applying Young's $L^1-L^1$ convolution inequality, and rescaling $\xione \mapsto \frac{c_0}{c_t}\xione$, we deduce that the above is controlled by

\begin{equation}\label{eq-1z}
C \left\| \int_\R e^{iz_1\xione} e^{\frac{i\nu_t}{2}\xione^2}  f\left(0, \xione \right) d\xione \right\|_{L^1}.
\end{equation}

The $L^\infty$-norm of the Fourier transform of the Gaussian $e^{\frac{i\nu_t}{2}\xione^2}$ is $\nu_t^{-1/2}$. Therefore, applying Young's $L^\infty-L^1$ convolution inequality, we obtain that (<ref>) is controlled by $C\nu_t^{-1/2}\|\check f\|_{L_1}$. This concludes the proof.

I think this section should be rewritten; right now the proof of Lemma 5.5 is overly complicated and it introduces multiple notations that we do not use later.

Regarding more specifically the proof of Lemma 5.5: We have $\TT = D_t - L$ for some  operator $L$ and that we know how to solve $(D_t + L)f = 0$ with an initial data $f_0$ – I write the solution $e^{-itL} f_0$. We specifically have (<ref>)
which pretty straightforwardly implies

\begin{equation}
\big\| e^{-itL}f \big\|_{\SSS_p} \leq C \lr{\nu_t}^p |f_0|_{\SSS_p}.
\end{equation}

In particular

\begin{equation}\label{eq-3n}
\big\| e^{-itL}f \big\|_{\mC \SSS_p} \leq C \lr{\nu_t}^p_T |f_0|_{\SSS_p}.
\end{equation}

Now, as you mention, the solution to $(D_t+L) f = g$ with vanishing initial data is provided by Duhamel's formula:

\begin{equation}
f(t,\xi) = \int_0^t e^{-i(t-s) L} g(s) ds.
\end{equation}

In particular, thanks to (<ref>) we have the estimate

\begin{equation}
\big\| f \big\|_{\mC \SSS_p} \leq T \big\| e^{-i(t-s)L} g\big\|_{\mC\SSS_p} \leq C(1+T) \lr{\nu_t}_T^p \|g\|_{\mC\SSS_p}
\end{equation}

This is a full proof of Lemma 5.5, and it does not require new notations besides something for the propagator $e^{itL}$.

Don't we have time-dependent coefficients? The integrating factors on $[s,t]$ are not of the form $(t-s)\cdots$ but rather integrals from $s$ to $t$. How do you want to simplify this? I'm interested. The notation is only used for that proof.

OK, you are right that the coefficients depend on time. But I don't think it changes much? Let me make corrections below:

We have $\TT = D_t - L_t$ for some  operator $L_t$ with time-dependent coefficients and that we know how to solve $(D_t + L_t)f_t = 0$ with given data $f_s$ at time $\tau$ – I write the solution $U(t,\tau) f_\tau$. We have an explicit formula for $U(t,\tau)$, which is obtained by replacing $t=0$ by $\tau$ in the quantities introduced in 5.2: $c_0$ becomes $c_\tau$ and the domain of integration in the formula (<ref>) for $\lambda_t$, $\nu_t$ changes from $[0,t]$ to $[\tau,t]$; in particular these get changed to $\nu_t-\nu_\tau$ and $\lambda_t-\lambda_\tau$, respectively. We then have an analog of the formula (<ref>) for $U(t,\tau)$, where $c_0, \lambda_t, \nu_t$ are changed accordingly. It implies that $U(t,\tau)$ is bounded on $\SSS_p$, with norm controlled by $\lr{\nu_{\tau,t}t-\nu_\tau}^p$.

The solution to $(D_t+L_t) f_t = g_t$ with $f_0 = 0$ is provided by Duhamel's formula:

\begin{equation}
f_t = \int_0^t U(t,\tau) g_\tau d\tau.
\end{equation}

In particular, thanks to (<ref>) we have the estimate

\begin{equation}
\big\| f_t \big\|_{\SSS_p} \leq t \sup_{\tau \in [0,t]} \big\| U(t,\tau) g_\tau \big\|_{\SSS_p} \leq t \sup_{\tau \in [0,t]} \lr{\nu_t-\nu_\tau}^p \|g_\tau\|_{\SSS_p} .
\end{equation}

Take the supremum for $t \in [0,T]$ to obtain

\begin{equation}
\big\| f \big\|_{\mC\SSS_p} \leq T \sup_{t,\tau \in [0,T]} \lr{\nu_{t,\tau}}^p \|g_\tau\|_{\SSS_p} \leq T \|g\|_{\mC\SSS_p} \sup_{t,\tau \in [0,T]} \lr{\nu_t-\nu_\tau}^p .
\end{equation}

Please, let's not introduce a new set of spaces and make things even more confusing... as written above, a reasonable estimate on $\big\| f \big\|_{\mC\SSS_p}$ is fairly easy to get... What you are writing above is roughly what I'm doing, both to invert $\mT$ and to get estimates. I'm just introducing notation that explains what is happening. I do not understand your phobia of notation. Too much is too much. But not enough is definitely not enough. Your estimate on $\|f_t\|_p$ requires justification which is painful; we are not just taking derivatives. I agree it is easy to state it and roughly get it. The $\mCS$ norm is not only taking a supremum in time so getting the final estimate is a bit more involved. This estimate is the reason for the loss of two ranks for each time derivative... Writing all of this is painful and not  the most interesting, I agree. But it needs to be written and the best way I found was to use the explicit expression for the norms as in (5.16), which is an interesting result in itself. Once you accept this (classical) form of the spaces, the proofs are reasonably simple and that should make it into the paper unless we find a reference in the literature. It might be more minor but we also need to bound the other operators. We have to decompose $T_1$ into its part linear in $D_\xi$ and the rest. Same for $T_0^{-1}$. Otherwise, we do not get optimal results. Look at the structure of the complete derivation without going to details before we can decide what stays or not. The scaling with $\vsa$ is very natural in our context. We can probably write things without it, but I'm not sure without the explicit form of the $\|\cdot\|_p$ norms. I'll think about this.

§.§ Dispersion-dependent functional setting

The solution of $\mT f=0$ with given initial condition admits an explicit expression as we saw in Lemma <ref>. In the analysis of the asymptotic expansion of $\psi(t,z)$ in powers of $\eps$, we need to solve transport equations of the form $\mT f=g$ with time-dependent source terms. To quantify the stability of the inverse transport operators and that of other relevant transforms, we introduce the following functional setting.

We recall the functional spaces $\mS_p$ were defined in section <ref>. To handle the time-dependence for $t\in [0,\rT]$ of the wavepackets, we introduce the spaces for $p\in\Nm$ and $k=\lfloor \frac p2 \rfloor \in \Nm$ (i.e., $p=2k$ or $p=2k+1$) defined for an interval ${\rm I}\subset\Rm$ by

\begin{equation}\label{eq:CSp}
\mCS_p({\rm I},\Rm^d,\Cm^q) = \cap_{r=0}^k C^r({\rm I}; \mS_{p-2r}(\Rm^d,\Cm^q))
\end{equation}

with norm given by the sum of the natural norms for the above spaces (see (<ref>) below). The spaces are constructed so that any derivative in time corresponds to a loss of order $2$ in the remaining variable.

To quantify the effects of dispersion, we define

\begin{align} \label{eq:anut}
\anut= 1+\sup_{0\leq t\leq \rT} |\nu_{t}|  \quad \mbox{ and } \quad \vsa=\anut^{-\frac12}.
\end{align}

We saw in Lemma <ref> that dispersion resulted in a multiplication operator of the form $\phi_t(\xi)=e^{i\frac12 \nu(t)\xi^2}$.
The operator $D_t \phi_t = \phi_t (D_t+\frac12 \nu'(t)\xi^2)$. This explains why the spaces $\mCS_p$ are constructed so that both $\partial_t$ and $\xi^2$ map  $\mCS_{p+2}$ to $\mCS_p$.

Similarly, $\partial_\xi \phi_t = \phi_t(\partial_\xi+i\nu(t)\xi)$, so that (the operator of multiplication by) $\phi_t$ is large as an operator on $\mCS_1$ when $\anut$ is.
Yet clearly, $\phi_t^2$ is comparable to $\phi_t$ in the same sense. Since the construction of our wavepackets requires repeated application of operators of the form $\phi_t$, we introduce scaled metrics on $\mS_p$ and $\mCS_p$ so that application of $\phi_t$ results in a bounded operation independent of $\anut$.

This is simply achieved by replacing $\fa_\xi^*=-\partial_\xi+\xi$ by $\fa_{\xi\vsa}^*=-\vsa \partial_\xi + \vsa^{-1}\xi$ and endowing the spaces $\mS_p(\Rm)$ and $\mS_p(\Rm^2)$ respectively with the norms

\begin{align} \label{eq:normsmS}
\|\psi\|^2_{p\vsa} = \|(\fa_{\xi\vsa}^*)^{p} \psi\|^2_0 \quad \mbox{ and } \quad \|\psi\|_{p\vsa}^2 =  \dsum_{j=0}^{p} \|(\fa_{\xi\vsa}^*)^j (\fa_\zeta^*)^{p-j}\psi \|^2_0.
\end{align}

We call $\mS_{p\vsa}$ the spaces $\mS_p$ endowed with these dispersion-scaled norms.

The space $\mCS_p({\rm I},\Rm^d,\Cm)$  are similarly endowed for $d=1,2$ with the norms (and called $\mCS_{p\vsa}$)

\begin{align}\label{eq:normsmCS}
\|\psi\|^2_{p\vsa} = \sup_{t\in {\rm I}} \dsum_{r=0}^{\lfloor \frac p2 \rfloor} \|(\fa_{\xi\vsa}^*)^{p-2r} \partial_t^r \psi\|^2_0 \ \mbox{ and } \ \|\psi\|_{p\vsa}^2 =  \sup_{t\in {\rm I}} \dsum_{r=0}^{\lfloor \frac p2 \rfloor} \dsum_{j=0}^{p-2r} \|(\fa_{\xi\vsa}^*)^j (\fa_\zeta^*)^{p-2r-j} \partial_t^r\psi \|^2_0.
\end{align}

Spaces of vector-valued functions are similarly constructed componentwise.

Note that $\vsa^{p}\|f\|_p \lesssim  \|f\|_{p\vsa} \lesssim  \vsa^{-p} \|f\|_{p}$. Here and below, we use the notation $a\lesssim b$ to mean the existence of a $\vsa-$independent constant $C$ such that $a\leq Cb$.
Here is a number of dispersion-dependent estimates we will be using.

All operator bounds below are meant to be $\vsa$-independent bounds.

Any operator $B\in \{ \vsa D_\xi, \vsa^{-1}\xi, D_\zeta, \zeta\}$ is bounded from $\mS_{p+1,\vsa}$ to $\mS_{p\vsa}$ and from $\mCS_{p+1,\vsa}$ to $\mCS_{p\vsa}$.  The operator $\partial_t$ is bounded from $\mCS_{p+2,\vsa}$ to $\mCS_{p\vsa}$.

Let $\phi(\xi)=e^{i\frac12 \nu \xi^2}$ with $|\nu|\lesssim \anut$. Then the operator of multiplication by $\phi(\xi)$ is bounded from $\mS_{p\vsa}$ to itself.

Let $\phi(t,\xi)=e^{i\frac12 \nu(t) \xi^2}$ with $\sup_{t\in {\rm I}}|\nu(t)|\lesssim \anut$ and $\sup_{t\in {\rm I}} |\nu^{(j)}(t)|\leq C_j$ for $j\geq1$. Then the operator of multiplication by $\phi(t,\xi)$ is bounded from $\mCS_{p\vsa}$ to itself.

The first statement comes from the construction of the spaces since $\fa^*_{\xi\vsa}$ controls $\vsa\partial_\xi$ and $\vsa^{-1}\xi$ in the sense that
\begin{align}
\|\fa^*_{\xi\vsa} \psi\|_0^2 = \vsa^2\|\partial_\xi \psi\|_0^2 + \vsa^{-2}\|\xi \psi\|_0^2 + \|\psi\|^2_0.
\end{align}

One obtains from (<ref>) for $p\geq2$ that $\|\partial_t\psi\|^2_{p\vsa}\lesssim \|\psi\|^2_{p-2,\vsa}$ and hence the bound on $\partial_t$.

Consider now the operator of multiplication by $\phi(\xi)$ in one dimension $d=1$. We wish to show that

\begin{align}
\|\phi\psi\|_{p\vsa}^2 =  \|(\fa_{\xi\vsa}^*)^p (\phi\psi)\|^2_0 \lesssim  \|(\fa_{\xi\vsa}^*)^p \psi\|^2_0 = \|\psi\|_{p\vsa}^2.
\end{align}

This holds when $p=0$. Assume it holds for $p-1\geq0$. Then

\begin{align}
\|\phi\psi\|_{p\vsa}^2 =  \|(\fa_{\xi\vsa}^*)^{p-1} (\fa_{\xi\vsa}^*)(\phi\psi)\|^2_0 =\|(\fa_{\xi\vsa}^*)^{p-1} (\phi\psi_1)\|^2_0 \lesssim \|\psi_1\|_{p-1,\vsa}^2
\end{align}

by induction hypothesis, where $\psi_1(\xi)=(\fa_{\xi\vsa}^*-i\vsa \nu \xi)\psi$. By construction of the functional spaces and the above result for the operator $\vsa^{-1}\xi$ knowing that $\vsa |\nu|\lesssim \vsa^{-1}$, we find $\|\psi_1\|_{p-1,\vsa}^2\lesssim \|\psi\|_{p,\vsa}^2$ and the result is proved. The same proof applies in two dimensions $d=2$ as well using the norm for $\mS_p(\Rm^2)$ in (<ref>) and that the commutator $[\fa_\zeta^*,\phi]=0$.

The proof in the time-dependent setting uses that

\begin{align}
\partial_t^r (\phi\psi)= \phi \big(\partial_t+i\frac12 \nu'(t)\xi^2\big)^r \psi = \phi \psi_r,\qquad \psi_r=\big(\dsum_{j=0}^r \nu_j(t) \xi^{2(r-j)} \partial_t^j \big) \psi
\end{align}

for smooth and bounded functions $\nu_j(t)$ independent of $\anut$ by assumption on $\nu(t)$. Assume dimension $d=1$ as $d=2$ is treated similarly.  Using that $\xi^{2(r-j)} \partial_t^j$ maps $\mCS_{p\vsa}$ to $\mCS_{p-2r,\vsa}$ for each $0\leq j\leq r$ and the bounds proved above in the time-independent setting, we find

\begin{align}
\|(\fa_{\xi\vsa}^*)^{p-2r} \partial_t^r \phi\psi\|_0 = \|(\fa_{\xi\vsa}^*)^{p-2r} \phi \psi_r\|_0 \lesssim \|\psi_r\|_{p-2r,\vsa} \lesssim \|\psi\|_{p\vsa}.
\end{align}

The explicit expression of the norms in (<ref>) and the above estimate conclude the proof of the lemma.

§.§ Stability of the transport and other operators

We then have the following stability result for the transport equation $\mT f=g$. We recall that $\vsa=\anut^{-\frac12}$.

Let $p\in \Nm$. The solution $f_0(t,\xione)$ of $\mT f_0=g$ on $[0,\rT]\times\Rm$ with initial condition $f_0(0,\xione)=\fhat(\xione)\in \mS_{p\vsa}(\Rm,\Cm)$ and source term $g\in \mCS_{p\vsa}(\Rm,\Cm)$ satisfies the estimate

\begin{equation}
%\|f_0\|_{{\mCS}_{pq}} \leq C_{pq} \big( \|f\|_{\mS_{pq}} +  \aver{\rT} \|g\|_{\mCS_{pq}}\big)   + C_{pq} \anut^p \big( \|f\|_{\mS_{0,q+p}} +  \aver{\rT} \|g\|_{\mCS_{0,q+p}}\big),
\|f_0\|_{p\vsa} \leq C_{p} \big(\|\fhat\|_{p\vsa} + \aver{\rT} \|g\|_{p\vsa}\big)
\end{equation}

where $\aver{\rT}=1+\rT$ and $C_{p}$ is independent of $\rT$.

We first adapt Lemma <ref> to handle volume sources and define

\[
\lambda_{s,t} = \int_s^t \frac{k_\tau}{2\rho_\tau}d\tau,\qquad \tilde\nu_{s,t} = 2 \dint_s^t \frac{c_t^2}{c_\tau^2}\big( \dot\theta_\tau\gamma_\tau + k_\tau \gamma_\tau^2)d\tau,\qquad e^{\mu_{s,t}} = \frac{c_t}{c_s},
\]

to obtain that the solution to $\mT f=0$ with initial condition $f(s,\cdot)$ is given by

\begin{equation}\label{eq:mTstot}
f(t,\xione) =  \exp\left( i \lambda_{s,t} + i\frac{\tilde \nu_{s,t}}{2} \xione^2 \right)  e^{\mu_{s,t}/2} f\left(s, e^{\mu_{s,t}}\xione \right).
\end{equation}

By an application of the Duhamel principle, the solution to $\mT f=g$ for $0\leq t\leq T$ with $f(0,\xione)=0$ is thus given explicitly by

\begin{equation}\label{eq:solmTsource}
f(t,\xione)=i\dint_0^t \exp\Big(i\lambda_{s,t}+i\frac12 \tilde \nu_{s,t} \xione^2\Big) e^{\mu_{s,t}/2} g(s, e^{\mu_{s,t}}\xione) ds.
\end{equation}

Consider the solution in (<ref>) at a fixed time $s$ with the explicit time-dependence $\lambda(t)=\lambda_{s,t}$, $e^{\mu(t)}=e^{\mu_{s,t}}$ and $\nu(t)=\tilde \nu_{s,t}$ to simplify notation. We prove the lemma for the operator

\begin{align}
f(t,\xi)  \to  e^{i\lambda(t)} e^{\frac{\mu(t)}2} e^{i\frac12 \nu(t)\xi^2} f(t,e^{\mu(t)}\xi).
\end{align}

The term $e^{i\lambda(t)}$ generates a smooth in time modulation. Since $\mu(t)$ is smooth and bounded above and below by positive constants independent of $t$, $e^{\mu(t)/2}$ is also smooth.  The operator $f(t,\xi)\to e^{\frac{\mu(t)}2}e^{i\lambda(t)} f(t,\xi)$ is therefore bounded in the norms $\|\cdot\|_{p\vsa}$.

All other time dependent coefficients are smooth and uniformly bounded in time independent of $\rT$ except for $\nu(t)$ that may grow linearly with $\nu'(t)$ uniformly bounded. Consider the operator $f(t,\xi)\mapsto h(t,\xi)=f(t,e^{\mu(t)}\xi)$ with $\fhat\in \mCS_{p\vsa}$. Plugging $h(t,\xi)$ into the definition of the norms (<ref>), we directly obtain that $\|h\|_{p\vsa}\lesssim \|f\|_{p\vsa}$.

The final transformation (with $\nu(t)$  replaced by $\nu(t)e^{-2\mu(t)}$)

\begin{align}
f(t,\xi) \mapsto h(t,\xi) = f(t,\xi) e^{i\frac12 \nu(t) \xi^2}
\end{align}

was analyzed in  Lemma <ref>. This concludes the analysis of the map from $\fhat(\xi)$ to $f_0(t,\xi)$.

The volume source term $g(t,\xi)$ is treated similarly by the Duhamel principle, with an additional possible integration in time that provides the extra multiplication by $\aver{\rT}$. This concludes the proof of the lemma.

We conclude this section with a summary of the operators we introduced to construct approximations of the Dirac equation and some estimates they satisfy.

The terms in $T_0$ and $T_1$ that contribute to the construction of the leading term $a_0(t,\Xitot)$ in the formal expansion $a=\sum_{j\geq0} \eps^{\frac j2} a_j$ were given in Lemma <ref>. Constructing higher-order terms $a_j$ and proving convergence results require estimates on the operators $T_j$, which are constructed as in the proof of Lemma <ref> and given explicitly by

\begin{equation}\label{eq:T1explicit}
T_1 = D_t + \mA + \tilde T_1,\qquad T_j=\tilde T_j,\quad j\geq2,
\end{equation}

with $\fU_t^*D_t \fU_t=D_t+\mA$ where

\begin{equation} \label{eq:mA}
\mA:=- \frac 12(\dot\varphi_t\sigma_2+\dot\theta_t(-s_t\sigma_1+c_t\sigma_3)-\dot\theta_t (-D_\xione D_\zetatwo -\zetatwo\xione+\gamma_t(\xione^2-D_\zetatwo^2))) + \dot\gamma_t \xione D_\zetatwo,
\end{equation}

and for $j\geq1$,

\begin{equation}\label{eq:Tjexplicit}
\tilde T_j =\dsum_{|\alpha|=j+1}\mV_t^* z^\alpha \mV_t\nu_\alpha\cdot\sigma, \qquad \Cm^3\ni
%\nu_\alpha := \dfrac{1}{\alpha!} \tilde R_{2,\varphi_t} \tilde R_{3,\theta_t}\partial^\alpha (\mR_{\theta_t}^*h)(R_{\theta_t} y_t).
\nu_\alpha = \frac{1}{\alpha!} \tilde R_{2,\varphi_t}\tilde R_{3,\theta_t} (R_{-\theta_t}\nabla)^\alpha h(y_t).
\end{equation}

We verify for completeness that $\mV_t^* z^\alpha \mV_t =(\zetatwo-\gamma_t\xione)^{\alpha_2}(-1)^{\alpha_1} (D_\zetatwo+\gamma_tD_\xione)^{\alpha_1}$.

We recall that $(j_t,k_t)=(\nu_{111},\nu_{021})$ while $\nu_{201}=0$ was important to obtain an explicit expression for the transport solution since the transport operator involves no term of the form $D^2_{\xione}$. We now summarize the estimates we will be needing.

We may write $T_0^{-1}=T_{01}^{-1}+T_{02}^{-1}$ and $T_1=T_{11}+T_{12}$ such that the following operators are bounded with $\anut-$independent bounds:

\begin{align}
\vsa^{-1}T_{01}^{-1}:& \qquad N_t^\perp\cap\mCS_{p+1,\vsa} \to \mCS_{p\vsa},  &
T_{02}^{-1}:& \qquad N_t^\perp\cap \mCS_{p\vsa} \to \mCS_{p\vsa}, \\ %%%N_t^\perp\cap \mCS_{p\vsa} \\ %true?
\vsa T_{11}:& \qquad  \mCS_{p+1,\vsa} \to \mCS_{p\vsa},  &
T_{12}:& \qquad  \mCS_{p+2,\vsa} \to \mCS_{p\vsa}, \\
\vsa^{j+1} T_{j}:& \qquad  \mCS_{p+j+1,\vsa} \to \mCS_{p\vsa} ,\quad j\geq2,  &
\aver{T}^{-1} \mT^{-1} :& \qquad  \mCS_{p\vsa} \to \mCS_{p\vsa}.
\end{align}

From Lemma <ref> and the decomposition

\begin{align}
T_{01}^{-1} = Q \matrice{ 0 & 0 \\ 0 & - 2\xi c_t (\rho_t\zeta+\partial_\zeta)^{-1}(\rho_t\zeta-\partial_\zeta)^{-1} } Q ,\quad T_{02}^{-1}= T_0^{-1}-T_{01}^{-1},
\end{align}

and Lemma <ref>, we obtain the above first two bounds.

We define $T_1=T_{11}+T_{12}$ with $T_{11}$ the contribution that is linear in $D_\xi$ while $T_{12}$ accounts for the rest (no contribution in $T_1$ is quadratic in $D_\xi$). The above corresponding bounds then follow from Lemma <ref> for quadratic expressions in $D_\xi$, $D_\zeta$, $\xi$, and $\zeta$. The same lemma is used to bound $T_j$ for $j\geq2$. The final estimate is a repeat of Lemma <ref>.

It seems that there is no motivation for this lemma? I think we should conclude this section with a result of the following form: Which lemma? Section 5 really is on the transport equation and its solution.

Section 6 on the construction of the approximation and its error estimates. What you write below is in section 6.1 as constructing the various terms $a_j$ takes a bit of time.  It is better in section 6.

Given $f$, the hierarchy of equations

\begin{equation}
T_0 a_0 = 0, \quad T_1 a_0 + T_0 a_1 = 0, \quad T_2 a_0 + T_1 a_1 + T_0 a_2 = 0, \dots
\end{equation}

with $a_0(0,\xi,\zeta) = \fhat(\xi) \phi(\zeta)$ admits a solution $(a_0, a_1, a_2, \dots)$. Moreover, this solution satisfies the estimate

\begin{equation}
\|a_j\|_{\m\SSS_p} \leq C (T \lr{\nu}_T)^j \|\fhat\|_{\mC \SSS_p}.
\end{equation}

That seems very optimistic in the power of $\nu$ but I agree there should be a formula that works for all $j$ rather than only $J=1$. That formula will be quite ugly though.

The solution (<ref>) involves (i) a rescaling by $\mu_t$ in the Fourier variables or $\mu_t^{-1}$ in the physical variables that preserves the $L^2-$norm in $\xi_1$ as expected; (ii) a physically irrelevant phase shift by $\lambda_0$; and (iii) a multiplication by $e^{i\mu_2(t)\xi_1^2}$ corresponding to a convolution by $e^{i\mu_2^{-1}z_1^2/4}$ in the physical domain.

Do we agree with the above three effects?  I agree with them, but I think we should describe their manifestation arises in physical space instead. For instance there is (i) compression by the opposite factor; (ii) a phase shift; and (iii) dispersion at rate $\mu^{-1/2}$.

I'm not sure we need to verify all these effects in the numerics.

I heard about compressions in some numerical cases. This is the only term that I see contributing to this behavior. Since $\lambda_1$ may be proportional to time, the compression is potentially exponential in time. In the cases I have looked into so far it was always decaying like a power-law. Any suggestions for an exponential one? I believe that exponential compression is impossible (too bad!) because the compression factor is in fact $c_t = \cos(\phi_t)$, hence bounded above (by $1$) and below (by a positive number if $r_t$ does not go to $0$). I am guessing that the power law decay that you observe numerically has nothing to do with this compression effect. I think it is a dispersive effect coming from the term $\nu_t \xi_1^2$ in (<ref>) (hence the maximal amplitude decays like $\nu_t^{-1/2}$).

This refers to the profile along the edge, which becomes $\tau_t f(\tau_t x_1)$ and $\tau_t$ can be exponentially large in $t$ leading to a compression/dilation (and hence possible drop in amplitude) of the wavepacket along the direction of propagation. This is independent of the profile of $f(x)$ itself. I thought the profiles in the simulations were close to Gaussian.Thank you, I should have explained better. I rather meant for which magnetic potential? Well, that's a good question. $\mu_1$ is related to the integral of $j_t\gamma_t$ (with $j_t=\nu_{021}$ I believe). So any magnetic potential such that $\mu_1$ changes should lead to compression or dilation.  Right, last time, we tried a linear magnetic field $B(x) = B_0 x$ and that led to a power-law decay. Will look into this again, but I think anything that is not constant and not a power-law would be interesting, but I am not sure sth like that even exists $j_t=r_tc_t\partial_\tau\beta$ according to Alexis' formula. What is decaying like a power law? The max amplitude (this is what we are talking about) or the spatial profile of the wavepacket? If the former, then this decay can indeed be understood when $j_t$ is non-constant. Since we are looking at short intervals of time Fixed: so we can make them fixed-long and see decay?, I'm not sure we'll be able to distinguish exponential from power law decays/increases. Indeed, so far we looked at the amplitude... Okay, $r_t$ does not have anything to do with the magnetic field..., $c_t$ is just inversely proportional to $\vert B_t\vert$...so it is really about understanding $\partial_{\tau}\beta$ then... That's pretty simple:  $\p_\tau \beta = \p_\tau \big(B/|\nabla \kappa|\big)$: $\tau$ is tangent to $\Gamma$ and $\beta  = B/|\nabla \kappa|$ along $\Gamma$. ($\p_n \beta$ is the complicated one)

We have to do something better here... Sorry, but otherwise there is no point in me reading the paper. Right now Theorem 6.4 tells me I need to look at Theorem 6.2 to understand which objects I am looking at. Theorem 6.2 tells me (vaguely) that $a^J$ is constructed in section 6.1. I looked at it, I find $a^J=\sum_{j=0}^J\eps^{\frac j2}a_j$,  of order $\eps^{\frac{J+1}2}$ of $a$ solution of $Ta=0$. What is $T$?

Anyway, we cannot hope that anybody, including myself, can read this and not get demotivated..

Here is an attempt at improving things.

I still do not think this is good enough. Really, this section should be just a couple pages long... right now it is 7 (!) pages long, and quite chaotic: what is even the main result?

Below is an undetailed attempt at streamlining the process. OK. I let you finish. The structure is relatively straightforward. (i) construct $a_j$ (ii) get approximation on $\psi-\psi^J$ (iii) move back to $\tilde\Psi$ and $\Psi$, which are the only Theorems stated in that section. One could indeed start with the final result but quite frankly they are not hard to spot at all. The main simplification would be to merge props 6.2&6.3 treating all $J$'s at once rather than handling the effect of dispersion only when $J=1$. Besides this, the derivation is fairly straightforward. cI agree that it is straightforward; but right now this is not well conveyed.

Let $a_j$ be the functions constructed in Lemma <ref>, and define

\begin{equation}\label{eq-3o}
\tilde \Psi_J = S\fU_t \sum_{j=0}^J \epsi^{j/2} a_j, \qquad \Psi_J = e^{i\chi/\epsi} \tilde \Psi_J.
\end{equation}

I renormalized $S$ below by removing the factor $\epsi^{1/2}$. That would require a lot of changes throughout the paper (which I am not keen on doing) starting from the Taylor expansion of $Sh$. BTW that's if $S$ says $S$ and does not become $S^{-1}$. where $\fU_t = \bU_t \mV_t$ consists of the transformation introduced in <ref> and $S$ is the rescaling (<ref>). With $R$ the remainder emerging in (<ref>), we have

\begin{equation}
\left(\epsi D_t + \Di\right) \Psi_J = e^{i\chi/\epsi} \left(\epsi D_t + \tDi + R\right) \tilde \Psi_J.
\end{equation}

Moreover, using the definition (<ref>) of $\tilde \Psi_J$, we have:

\begin{align}
\left(\epsi D_t + \tDi + R\right) \tilde \Psi_J & = \left(\epsi D_t + \tDi \right) S \fU_t \sum_{j=0}^J \epsi^{j/2} a_j + R \Psi^J
\\