\begin{bmatrix} \frac{\partial  a}{\partial \zeta}(\zeta_0)(\delta \zeta) \\ \frac{\partial  c}{\partial \zeta}(\zeta_0)(\delta \zeta) \end{bmatrix}
\end{equation}

where the coefficient matrix has the following components:  $a_{11} = a_{11}^0+b_{11}$, where

a_{11}^0 =  \left\la \frac{\partial Q_{ a,  c}}{\partial  a} \Big|_{( a(\zeta_0),  c(\zeta_0))} , Q_{ a(\zeta_0),  c(\zeta_0)} \right\ra  = \frac{\partial}{\partial a}\Big|_{a(\zeta_0)} \|Q_{a,c(\zeta_0)} \|_{L_x^2}^2 =0

$$b_{11}=\left\la \frac{\partial Q_{ a,  c}}{\partial  a} \Big|_{( a(\zeta_0),  c(\zeta_0))} , \eta_0 \right\ra \quad \implies \quad |b_{11}| \leq \|\eta_0\|_{L_x^2}$$
Next, $a_{12}=a_{12}^0+b_{12}$, where

\begin{align*}
a_{12}^0 &=  \left\la \frac{\partial Q_{ a,  c}}{\partial  c} \Big|_{( a(\zeta_0),  c(\zeta_0))} , Q_{ a(\zeta_0),  c(\zeta_0)} \right\ra \\
&= \frac12 \frac{\partial}{\partial c}\Big|_{c(\zeta)} \| Q_{a,c} \|_{L_x^2}^2 = \frac12 \frac{\partial}{\partial c}\Big|_{c(\zeta)} \left( c \|Q\|_{L_x^2}^2 \right) =\frac12 \|Q\|_{L_x^2}^2
\end{align*}

$$b_{12} =  \left\la \frac{\partial Q_{ a,  c}}{\partial  c} \Big|_{( a(\zeta_0),  c(\zeta_0))} , \eta_0\right\ra
\quad \implies \quad |b_{12}| \leq \|\eta_0\|_{L_x^2}$$

Next, $a_{21} = a_{21}^0 + b_{21}$, where

\begin{align*}
a_{21}^0 &= \left\la (x- a(\zeta)) \frac{\partial Q_{ a,  c}}{\partial  a} \Big|_{( a(\zeta_0),  c(\zeta_0))} ,  Q_{ a(\zeta_0),  c(\zeta_0)} \right\ra \\
&= \frac12\int (x-a) \frac{\partial}{\partial a} [ Q_{a,c}(x)^2 ]\, dx  \Big|_{(a(\zeta_0),c(\zeta_0)}
= -\frac12\int (x-a) \frac{\partial}{\partial x} [ Q_{a,c}(x)^2 ]\, dx \Big|_{(a(\zeta_0),c(\zeta_0)} \\
&= \frac12 \| Q_{a(\zeta_0),c(\zeta_0)} \|_{L_x^2}^2
\end{align*}

by integration by parts,

$$b_{21} = \left\la \frac{\partial}{\partial a}\Big|_{(a(\zeta_0),c(\zeta_0))} [(x-a(\zeta)) Q_{a,c}], \eta_0 \right\ra \qquad \implies \quad |b_{12}|\lesssim \|\eta_0\|_{L_x^2}$$

Finally $a_{22} = a_{22}^0+b_{22}$, where

\begin{align*}
a_{22}^0 &=
\left\la (x- a(\zeta))\frac{\partial Q_{ a,  c}}{\partial  c} \Big|_{( a(\zeta_0),  c(\zeta_0))} ,  Q_{ a(\zeta_0),  c(\zeta_0)}\right\ra \\
& = \frac12 \frac{\partial}{\partial c} \Big|_{(a(\zeta_0),c(\zeta_0))} \la (x-a) Q_{a,c}, Q_{a,c} \ra = 0
\end{align*}

b_{22} =
\left\la (x- a(\zeta))\frac{\partial Q_{ a,  c}}{\partial  c} \Big|_{( a(\zeta_0),  c(\zeta_0))} , \eta_0 \right\ra \qquad \implies \quad |b_{22}| \lesssim \|\eta_0\|_{L_x^2}$$

\begin{bmatrix}
a_{11} & a_{12} \\ a_{21} & a_{22}
\end{bmatrix}
= \frac12 \|Q\|_{L_x^2}^2 \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} + O( \|\eta_0\|_{L_x^2})$$

from which it follows that

\begin{bmatrix}
a_{11} & a_{12} \\ a_{21} & a_{22}
\end{bmatrix}^{-1}
= 2\|Q\|_{L_x^2}^{-2} \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} + O( \|\eta_0\|_{L_x^2})$$

Solving (<ref>) by inverting this $2\times 2$ matrix

\begin{equation}
\label{E:est-62}
\begin{bmatrix} \frac{\partial  a}{\partial \zeta}(\zeta_0)(\delta \zeta) \\ \frac{\partial  c}{\partial \zeta}(\zeta_0) (\delta \zeta) \end{bmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22}
\end{bmatrix}^{-1}
\begin{bmatrix}
\la \delta \zeta , Q_{ a(\zeta_0),  c(\zeta_0)} \ra \\
\la \delta \zeta , (x- a(\zeta_0))Q_{ a(\zeta_0),  c(\zeta_0)} \ra
\end{bmatrix}
\end{equation}

which  gives the needed components of (<ref>).  Combining (<ref>), (<ref>), (<ref>), and (<ref>), we obtain

$$\eta(x) = \zeta(x) + \int_{s=0}^1 Q_{a(s\zeta),c(s\zeta)}'(x) \frac{\partial a}{\partial\zeta} \Big|_{s\zeta}(\zeta) \, ds - \int_{s=0}^1 [(x-a)Q_{a,c}]\Big|_{a(s\zeta),c(s\zeta)}'(x)\frac{\partial c}{\partial\zeta} \Big|_{s\zeta}(\zeta)  \, ds$$

For each $\am$, $\cm$, and corresponding $a$, $c$,

\begin{equation}
\label{E:est-55}
\| \eta\|_{H_x^{1/2}} \lesssim \| \zeta\|_{H_x^{1/2}}
\end{equation}

Taking $\am(t)$, $\cm(t)$ and correspondingly $a(t)$, $c(t)$, along the flow[For this, we need only assume that $a(t)\sim t$ and $\frac12<c(t)<2$],

\begin{equation}
\label{E:coeff-conv03}
\sup_n \| \eta \|_{L_I^2 L_{x\in (n,n+1)}^2} \lesssim  (\ln h^{-1})\sup_n \| \zeta \|_{L_I^2L_{x\in (n,n+1)}^2}  + h^{1/2}\| \zeta\|_{L_I^2L_x^2}
\end{equation}

(<ref>) follows directly from (<ref>) and the two equations after (<ref>).  To prove (<ref>), we will use that for $a,b >0$ and $\alpha,\beta \in \mathbb{R}$,

$$\sup_t \la t - \alpha \ra^{-a} \la t - \beta \ra^{-b}  \lesssim \la \alpha - \beta \ra^{-\min(a,b)}$$

We know that $a(t) \sim t$.  Starting with (<ref>) (see the approximation (<ref>) to help with conceptualization), apply the $L_I^2L_{x\in (n,n+1)}^2$ norm for fixed $n$, and estimate as:

\begin{align*}
\| \eta \|_{L_I^2 L_{x\in (n,n+1)}^2} &\lesssim \| \zeta \|_{L_I^2 L_{x\in (n,n+1)}^2} + \int_x \| \zeta(x,t) \la n - a(t) \ra^{-2} \la x- a(t) \ra^{-1} \|_{L_{t\in I}^2} \, dx \\
&\lesssim \| \zeta \|_{L_I^2 L_{x\in (n,n+1)}^2} + \int_x \| \zeta(x,t) \|_{L_{t\in I}^2} \sup_{t\in I} \la n - a(t) \ra^{-2} \la x- a(t) \ra^{-1}  \, dx \\
&\lesssim \| \zeta \|_{L_I^2 L_{x\in (n,n+1)}^2} + \int_x \| \zeta(x,t) \|_{L_{t\in I}^2} \la n -x \ra^{-1}  \, dx
\end{align*}

Split the $x$-integral into $|x-n|< h^{-1}$ and $|x-n|>h^{-1}$.  The region $|x-n|<h^{-1}$ is divided into unit-sized $x$-pieces producing the factor $\sum_{|m|<h^{-1}} \la m\ra^{-1} \lesssim \ln h^{-1}$.  In the region $|x-n|>h^{-1}$, we apply Cauchy-Schwarz and use $\| \la x\ra^{-1} \|_{L_{|x|>h^{-1}}^2} \leq h^{1/2}$.  Together, this yields

$$\| \eta \|_{L_I^2 L_{x\in (n,n+1)}^2} \lesssim \| \zeta \|_{L_I^2 L_{x\in (n,n+1)}^2} + (\ln h^{-1})\sup_m \| \zeta \|_{L_I^2L_{x\in (m,m+1)}^2}  + h^{1/2}\| \zeta\|_{L_I^2L_x^2}

From this, (<ref>) follows.

Define the remainder $\eta$ according to

\begin{equation}
\label{E:decomp2}
u = Q_{a, c} + \eta
\end{equation}

imposing orthogonality conditions

\begin{equation}
\label{E:orth2}
\la \eta, Q_{a,c} \ra =0 \,, \qquad \la \eta, (x-a) Q_{a,c} \ra = 0
\end{equation}

An implicit function theorem argument shows that there exists a unique choice of $(a,c)$ so that these orthogonality conditions hold.  This is the definition of the parameters $(a(t),c(t))$ and of the remainder $\eta$.

Starting with $\partial_t u = JE'(u)$, we substitute (<ref>) to obtain

$$\partial_t (Q_{a,c} + \eta) = J E'( Q_{a,c} + \eta)$$

Analogously to the derivation of (<ref>), we find

\begin{equation}
\label{E:eta2}
\partial_t \eta = - \dot a \partial_a Q_{a,c} - \dot c \partial_c Q_{a,c} + J E'(Q_{a,c}) + JE''(Q_{a,c}) \eta  - \frac12 \partial_x (\eta^2)
\end{equation}

We re-center the equation for $\eta$ by letting

$$w (y) = \eta (y+a) \quad \iff \quad \eta (x) = w(x-a)$$

The orthogonality conditions on $w$ read

\begin{equation}
\label{E:w-ortho}
\la w, Q_c \ra =0 \,, \qquad \la w, yQ_c \ra =0 .
\end{equation}

The equation for $w$ is

\begin{equation}
\label{E:w}
\begin{aligned}
\partial_t w = \;
&(\dot a-c+W(ha))\partial_yQ_c +(-\dot c c^{-1} + hW'(ha))\partial_y(yQ_c) + \partial_y(e_2Q_c) \\
&+\partial_y \mathcal{L}_c v + \partial_y(\dot a - c + W(hx)) w - \frac12 \partial_y w^2
\end{aligned}
\end{equation}

Here, (<ref>) is analogous to (<ref>).

For all $t$, if $\frac12 \leq  c  \leq 2$ and $\|w \|_{L_y^2} \ll 1$, then

\begin{equation}
\label{E:par2b}
\begin{aligned}
& |\dot a  -  c  + W(h a ) - \tfrac12 h^2W''(h a )  c ^{-1}  | \lesssim h^3 + \|w\la y \ra^{-1} \|_{L_y^2}^2 \\
&|\dot c   - hW'(h a ) c  - \tfrac12 h^3 W''(h a ) c ^{-1}| \lesssim h^4 + h^2(\ln h^{-1})\sup_{n\in \mathbb{Z}} \|w\|_{L_{n<y < n+1}^2} + \|w\la y \ra^{-1} \|_{L_y^2}^2
\end{aligned}
\end{equation}

Also, for a time interval $I$,

\begin{equation}
\label{E:par2bt}
\begin{aligned}
& \int_I |\dot a  -  c  + W(h a ) - \tfrac12 h^2W''(h a )  c ^{-1}  |\, dt \lesssim h^3|I| + \|w\la y \ra^{-1} \|_{L_I^2L_y^2}^2 \\
&\int_I |\dot c   - hW'(h a ) c  - \tfrac12 h^3 W''(h a ) c ^{-1}|\, dt \\
& \qquad \qquad \lesssim h^4|I| + h^2(\ln h^{-1})|I|^{1/2}\sup_{n\in \mathbb{Z}} \|w\|_{L_I^2L_{n<y < n+1}^2} + \|w\la y \ra^{-1} \|_{L_I^2L_y^2}^2
\end{aligned}
\end{equation}

Taking $\partial_t$ of the orthogonality condition $\la w, Q_ c  \ra =0 $, then exactly as in the proof of Lemma <ref>, we obtain

$$|\dot c   - hW'(h a ) c   - \frac12 h^3 W''(h a ) c ^{-1}| (4\pi -  c ^{-1}\la w, \partial_y(yQ_ c )\ra) $$

$$\lesssim h^4 +  h^2 (\ln h^{-1}) \sup_n \|w\|_{L_{n\leq y\leq n+1}^2}+ \|\la y \ra^{-1} w\|_{L_y^2}^2$$

from which the second inequality in (<ref>) follows.  The second inequality in (<ref>) also follows as in the proof of Lemma <ref>.

Now we prove the first inequality in (<ref>).
Taking $\partial_t$ of the orthogonality condition $0=\la w, y Q_ c  \ra$, we obtain

$$0 = \la \partial_t w, y Q_ c  \ra + \la w, y\partial_t Q_ c  \ra$$

For the first term, we substitute (<ref>), and for second term, we use that $\partial_tQ_ c  = \dot c   c ^{-1} \partial_y (y Q_ c) $, to obtain

\begin{align*}
0 &= ( \dot a  -  c  + W(h a )) \la \partial_yQ_ c , yQ_ c \ra && \leftarrow \text{I}\\
&\quad +(-\dot  c   c ^{-1} +hW'(h a )) \la \partial_y(yQ_ c ),yQ_ c  \ra  && \leftarrow \text{II}\\
&\quad + \la \partial_y(e_2Q_ c ), yQ_ c  \ra && \leftarrow \text{III}\\
&\quad + \la \partial_y\mathcal{L}_ c  w, yQ_ c  \ra && \leftarrow \text{IV}\\
&\quad + \la \partial_y(\dot  a  -  c  + W(h(y+ a )))w, y Q_ c  \ra && \leftarrow \text{V}\\
&\quad - \tfrac12 \la \partial_y w^2, y Q_ c  \ra && \leftarrow \text{VI}\\
&\quad + \dot  c   c ^{-1} \la w, y\partial_y (yQ_ c ) \ra && \leftarrow \text{VII}
\end{align*}

Given that $\la \partial_y Q_c, y Q_c \ra = - \frac12\|Q_c\|_{L_y^2}^2$ and
$\|Q_ c \|_{L^2}^2 = 8\pi  c$, we have

$$\text{I} = -4\pi  c( \dot a  -  c  + W(h a ))$$

Also, given that $\partial_y(yQ_ c )$ is even and $yQ_ c $ is odd, we have $\la \partial_y(yQ_ c ),yQ_ c  \ra =0$ and thus $\text{II}=0$.

To address term III, we carry out the Taylor expansion

\begin{align*}
e_2 &= W(h(y+ a )) - W(h a ) - hW'(h a ) y \\
&= \frac12 h^2W''(h a )y^2 + \frac16 h^3W'''(h a ) y^3 + \frac1{24}h^4W''''(h(y_*+ a ))y^4
\end{align*}

for some $y_*$ between $0$ and $y$, by the Lagrange form of the remainder.  Substituting

\begin{align*}
\text{III} &= -\la e_2, Q_ c  \partial_y (yQ_ c)  \ra \\
&= - \frac12 h^2 W''(h a ) \int y^2 Q_ c  \partial_y (yQ_ c)  \, dy - \frac16h^3 W'''(h a ) \int y^3 Q_ c  \partial_y (yQ_ c)  \, dy \\
& \qquad - \frac1{24}h^4W''''(h a ) \int y^4 Q_ c  \partial_y (yQ_ c)
\end{align*}

Since $\int y^2 Q_ c  \partial_y (yQ_ c)  \, dy= -\frac12 \int y^2Q_c^2 \, dy = -\frac12 c^{-1} \int z^2Q^2 \, dz = -4\pi c^{-1}$,

$$\text{III} =  2\pi h^2W''(h a ) c^{-1}  + O(h^4)$$

To address Term IV, we use (<ref>) and orthogonality condition $\la w, Q_c \ra =0$:

$$\text{IV} = - \la w, \mathcal{L}_ c  \partial_y(yQ_c)  \ra = - c \la w, \mathcal{L}_c \partial_c Q_c \ra =  c\la w, Q_c \ra = 0$$

In Term V, we expand

$$W(h(y+ a ))= W(h a ) + hW'(h a )y + e_2(y, a )$$

to yield

$$\text{V} = (\dot  a  -  c  + W(h a ))\la \partial_y w, yQ_ c \ra - hW'(h a ) \la w, y \partial_y(y Q_ c)  \ra - \la e_2w, \partial_y(y Q_ c)  \ra$$

Note that the middle term combines with Term VII

$$\text{V}+\text{VII} = (\dot  a  -  c  + W(h a ))\la \partial_y w,  yQ_ c  \ra +(\dot c   c ^{-1} - hW'(h a )) \la w, y \partial_y(y Q_ c)  \ra - \la e_2w, \partial_y(y Q_ c) \ra$$

By Taylor's theorem with the Lagrange form of the remainder

$$e_2(y,a) = \frac12 h^2 W''(h(y_*+a))y^2$$

Let $R>0$ such that $\supp W \subset [-R,R]$.  Then $-a - Rh^{-1} \leq y \leq -a + Rh^{-1}$.  This gives

$$\la e_2 w, \partial_y(yQ_c) \ra = \frac12 h^2 \int_{-a-Rh^{-1}}^{-a+Rh^{-1}} W''(h(y_*+a) w y^2\partial_y(yQ_c) \, dy$$

Since $\| y^2\partial_y(yQ_c)\|_{L^\infty_y} \lesssim 1$ and $\|W''\|_{L^\infty_y}\lesssim 1$ Cauchy-Schwarz gives

$$|\la e_2 w, \partial_y(yQ_c) \ra| \lesssim h^2 \|w\|_{L^2} (Rh^{-1})^{1/2} \lesssim h^{3/2} \|w\|_{L^2_y}$$

$$\text{VI} = \frac12 \la w^2, \partial_y(yQ_c) \ra \lesssim \|w\la y \ra^{-1} \|_{L_y^2}^2$$

Combining the estimates for Terms I-VII, we obtain

\begin{align*}
\indentalign |(\dot a -c + W(ha)(1+ \frac{1}{4\pi c} \la w, \partial_y(yQ_c)\ra) - \frac12 h^2 W''(ha) c^{-2}| \\
&\lesssim h^3 + \|w\|_{L_y^2}^2 + |\dot c c^{-1} - hW'(ha)| \|w\|_{L_y^2}
\end{align*}

By the second inequality in (<ref>), $|\dot c c^{-1} - hW'(ha)|\lesssim h^2$, and therefore, the corresponding term in the inequality above can be bounded by the other terms.  From this it follows that

$$ |(\dot a -c + W(ha) - \frac12 h^2 W''(ha) c^{-2}| \lesssim h^3 + \|w\|_{L_y^2}^2$$

which is the first inequality in (<ref>).  The first inequality in (<ref>) follows from the first inequality in (<ref>) after integrating in $t$.

There exists $\kappa \geq 1$, $\mu>0$, and $0<h_0\ll 1$ such that the following holds.
Let $0< h \leq h_0$ and suppose the initial data $u_0\in H_x^1$ satisfies

$$\| u_0(x) - Q_{0,1}(x) \|_{H_x^{1/2}} \leq h^{3/2}$$

Suppose that $u$ satisfying (pBO) with initial condition $u(x,0)=u_0(x)$ is decomposed as (<ref>) with remainder $\eta$ satisfying orthogonality conditions (<ref>).    For every $T>0$ such that $\frac12 \leq c(t) \leq 2$ for all $0\leq t \leq T$,
we have that the recentered remainder $w(y,t) = \eta(y+a(t),t)$ satisfies

\begin{equation}
\label{E:symp-estimates-eq1}
\begin{aligned}
&\|w\|_{L_{[0,T]}^\infty H_y^{1/2}}  \leq \kappa h^{3/2} e^{\mu hT} \\
&\sup_n \| w\|_{L_{[0,T]}^2L_{y\in(n,n+1)}^2} \leq \kappa h^{3/2} (\ln h^{-1}) e^{\mu hT}
\end{aligned}
\end{equation}

and the parameters $a(t)$, $c(t)$ satisfy the following bounds (<ref>).

From (<ref>) and (<ref>) in Corollary <ref>, combined with (<ref>) in Proposition <ref>, we immediately obtain (<ref>).  The ODE bounds (<ref>) hold by Lemma <ref>.

Theorem <ref> can now be proved as a consequence of Proposition <ref>.

By Proposition <ref>, we have the estimate (<ref>) for $w$.
The parameters $(a(t),c(t))$ in Proposition <ref> satisfy the bounds in Lemma <ref>.  Define $(A(s),C(s))$ by $a(t)=h^{-1}A(ht)$ and $c(t) = C(ht)$.  Then by (<ref>) and (<ref>), $(A(s),C(s))$ satisfy

\begin{equation}
\label{E:ODE11}
\begin{aligned}
&\int_0^s |\dot { A} -  C+ W( A) +\frac12 C^{-2} h^2 W''(A)| \, d\tau \leq \kappa^2 h^3 (\log h^{-1}) e^{2\mu s}\\
&\int_0^s |\dot { C}  - CW'( A) - \frac12 C^{-2} h^2 W'''(A) |\, d\tau \leq \kappa^2 h^3 (\log h^{-1}) e^{2\mu s}
\end{aligned}
\end{equation}

on $0\leq s \leq  \min(\frac14\mu^{-1} \ln h^{-1},S_0)$.     Now apply Lemma <ref> on ODE perturbation to compare the $(A,C)$ parameter dynamics with the so-called exact trajectory $(\hat A, \hat C)$ defined in (<ref>).   Specifically, we obtain that
$| A - \hat A| \lesssim h^3 e^{2\mu s}$ and $|C - \hat C| \lesssim h^3e^{2\mu s}$, and thus $|a-\hat a| \lesssim  h^2e^{2\mu ht}$ and $|c-\hat c| \lesssim h^3 e^{2\mu ht}$.   These bounds imply

$$\| Q_{\hat a, \hat c} - Q_{a, c} \|_{H_x^{1/2}} \lesssim h^2 e^{2\mu ht}$$

$$\|u - Q_{\hat a, \hat c} \|_{H_x^{1/2}} \leq \| u - Q_{a,c} \|_{H_x^{1/2}} + \|Q_{a,c} - Q_{\hat a, \hat c} \|_{H_x^{1/2}} =  \|w \|_{H_x^{1/2}} + \|Q_{a,c} - Q_{\hat a, \hat c} \|_{H_x^{1/2}}$$

By (<ref>),

$$\|u - Q_{\hat a, \hat c} \|_{H_x^{1/2}} \lesssim h^{3/2} e^{\mu ht}$$

Thus Theorem <ref> follows.

Suppose $X, \bar X :\mathbb{R} \to \mathbb{R}^d$ solve

\begin{align*}
&\dot X(s) = f(X(s)) + h^2 g(X,s) \\
&\dot{\bar X}(s) = f(\bar X(s))
\end{align*}

with the same initial condition $X(0)=\bar X(0)$,
where $f: \mathbb{R}^d \to \mathbb{R}^d$, $ g: \mathbb{R}^{d+1} \to \mathbb{R}^d$.  Suppose that the $d\times d$ matrix $f'(X)$ is uniformly bounded:  for all $X\in \mathbb{R}^d$,

$$\| f'(X)\|_{\ell^2} \leq \kappa$$

where $\ell^2$ is the square sum norm on the $d^2$ entries of the matrix.  Then

$$|X(s)-\bar X(s)|^2 \leq  h^4 \int_0^s e^{-(2\kappa+1)(s-s')}|g(X(s'),s')|^2 \, ds'$$

Let $V(s) = X(s)-\bar X(s)$.  Then ($| \bullet |$ is the usual square sum norm on $\mathbb{R}^d$)

\begin{equation}
\label{E:ODE1}
\frac{d}{ds} |V|^2 = 2V\dot V = 2V\cdot (f(X)-f(\bar X)) + 2h^2 V\cdot g(X,s)
\end{equation}

We have

$$f(X) - f(\bar X) = \int_{\sigma=0}^1 \frac{d}{d\sigma}[ f(\bar X+ \sigma V)] \, d\sigma = \left( \int_{\sigma=0}^1 f'(\bar X+ \sigma V) \, d\sigma \right)V$$

Then by Cauchy-Schwarz,

$$| f(X) - f(\bar X)| \leq \kappa |V|$$

Substituting this into (<ref>), and using that $2h^2 V \cdot g(X,s) \leq |V|^2 + h^4|g(X,s)|^2$, we obtain

$$\frac{d}{ds} |V|^2 \leq (2\kappa+1) |V|^2 + h^4|g(X,s)|^2$$

The standard integrating factor method completes the proof.

In our application,

$$X = \begin{bmatrix} A \\ C \end{bmatrix} \,, \qquad f(X) = \begin{bmatrix}
C-W(A) \\ CW'(A) \end{bmatrix}$$

$$f'(X) = \begin{bmatrix} -W'(A)  & 1 \\ CW''(A) & W'(A) \end{bmatrix}$$

Since $\frac12\leq C \leq 2$, this is uniformly bounded.

§ LINEAR LIOUVILLE THEOREM FOR (BO) ASYMPTOTIC STABILITY

In this section, we will prove Theorem <ref>. By Theorem <ref>,

\begin{equation}
\label{E:D1Dmod}
\sup_{y_0\in \mathbb{R}} \|  \la D_y \ra^{1/2} ((g_{\gamma,y_0}')^{1/2} v ) \|_{L_{[0,T]}^2L_y^2}^2 \lesssim_\gamma  \|v\|_{L_{[0,T]}^\infty L_y^2}^2
\end{equation}

uniformly in $T>0$.  This implies the conveniently stated estimate

\begin{equation}
\label{E:linlio1}
\sup_{|I|=1} \| v \|^2_{L^2_{t\in \mathbb{R}} L^2_{y \in I }} \lesssim  \|v \|^2_{L^{\infty}_{t\in \mathbb{R}} L^2_y }
\end{equation}

where the supremum is taken over all unit-length intervals $I\subset \mathbb{R}$.  From (<ref>), we will obtain

We have

\begin{equation}
\label{E:temperedbound}
\int_{t\in \mathbb{R}} \frac{1}{\langle t-t_0 \rangle^{4/5}} \int_{y\in\mathbb{R}} v^2 (t,y) dy dt < \infty
\end{equation}

uniformly in $t_0\in \mathbb{R}$.

By translation in time, it suffices to assume that $t_0=0$.
Split the integral into

\int_t \frac{1}{\langle t \rangle^{4/5}} \int_{|y| < \langle t \rangle^{3/5}} v^2 (t,y) dy dt + \int_t \frac{1}{\langle t \rangle^{4/5}} \int_{|y| > \langle t \rangle^{3/5}} v^2 (t,y) dy dt := \text{I}+\text{II}

and compute

\begin{align*}
\text{I} & = \int_t \frac{1}{ \langle t \rangle^{4/5}} \sum_n \int_{y \in [n, n+1], |y| \leq \langle t \rangle^{3/5}} v^2 (t,y) dy dt   \\
& =  \sum_n \int_t \frac{1}{ \langle t \rangle^{4/5}} \int_{y \in [n, n+1], |y| \leq \langle t \rangle^{3/5}} v^2 (t,y) dy dt
\end{align*}

The condition on the inner integral implies that $\la n \ra \lesssim \la t \ra^{3/5}$, from which it follows that $\la t \ra^{-4/5} \leq \la n \ra^{-4/3}$.  Therefore, we can continue the estimate as

\text{I} \lesssim \sum_n \frac{1}{\la n\ra^{4/3}} \int_t  \int_{y \in [n, n+1], |y| \leq \langle t \rangle^{3/5}} v^2 (t,y) dy dt \leq \sum_n \frac{1}{\la n\ra^{4/3}}  \sup_I \|v \|^2_{L^2_{[0,T]} L^2_{y\in I} }  < \infty

by (<ref>).  Moreover, by the uniform spatial decay hypothesis (<ref>), we have

\int_{|y| > \langle t \rangle^{3/5}} v^2 (t,y) dy \lesssim \frac{1}{\langle t \rangle^{3/5}}

from which we obtain

\text{II} \lesssim \int_t \frac{1}{\langle t \rangle^{7/5}} dt < \infty

Since $\text{I}<\infty$ and $\text{II}<\infty$, (<ref>) holds.

By Proposition 2 on p. 920 of Kenig & Martel [20], there exists $A\gg 1$ such that with

\begin{equation}
\label{E:monphi}
\phi (y)= \frac{\pi}{2} + \arctan (\frac{y}{A})
\end{equation}

the following holds:  for any $\lambda \in (0,1)$, $t\leq t_0$, and $y_0 >1$, we have the monotonicity estimate

\begin{equation}
\label{E:mon}
\begin{aligned}
\indentalign \int v^2 (y,t_0) ( \phi (y- y_0) -  \phi (-y_0)  ) dy \\
& \leq \int v^2 (y,t) ( \phi (y- y_0 -\lambda (t_0 - t)) -  \phi (-y_0 - \lambda (t_0 - t))  ) dy \\
& \qquad + C \int^{t_0}_t \frac{\| v(t') \|^2_{L^2_y}  }{ (y_0 + \lambda (t_0 - t'))^2 } dt'  = p_1(t)+p_2(t)+p_3(t)
\end{aligned}
\end{equation}

We have decomposed the right-side as

$$p_1(t) = \int_{y> \frac12(y_0+\lambda(t_0-t))} v^2 (y,t) ( \phi (y- y_0 -\lambda (t_0 - t)) -  \phi (-y_0 - \lambda (t_0 - t))  ) dy$$

$$p_2(t) = \int_{y< \frac12(y_0+\lambda(t_0-t))} v^2 (y,t) ( \phi (y- y_0 -\lambda (t_0 - t)) -  \phi (-y_0 - \lambda (t_0 - t))  ) dy$$

$$p_3(t) = C\int^{t_0}_t \frac{\| v(t') \|^2_{L^2_y}  }{ (y_0 + \lambda (t_0 - t'))^2 } dt'$$

Note that

$$p_3(t) \lesssim \int_t^{t_0} \la t'-t_0\ra^{-4/5} \|v(t')\|_{L_y^2}^2 \, dt' \; \sup_{t'\leq t_0} \left[ \frac{\la t'-t_0\ra^{4/5}}{(y_0+\lambda(t_0-t'))^2} \right]$$

Thus by (<ref>), $p_3(t) \lesssim y_0^{-6/5}$ uniformly in $t<t_0$.  Next, we will show that $\lim_{t\to -\infty} p_1(t) = 0$ and  $\lim_{t\to -\infty} p_2(t) = 0$.  Indeed, by (<ref>),

$$|p_1(t)| \lesssim \frac{1}{y_0 + \lambda(t_0-t)} \quad \implies \quad \lim_{t\to -\infty} p_1(t) = 0$$

$$|p_2(t)| \lesssim \int_{y} v^2 (y,t)  dy \; \sup_{y<\frac12(y_0+\lambda(t_0-t))} ( \phi (y- y_0 -\lambda (t_0 - t)) -  \phi (-y_0 - \lambda (t_0 - t)))$$

and from the formula (<ref>) for $\phi(y)$,

$$\sup_{y<\frac12(y_0+\lambda(t_0-t))} ( \phi (y- y_0 -\lambda (t_0 - t)) -  \phi (-y_0 - \lambda (t_0 - t))) \leq 2 \phi(-\tfrac12(y_0+\lambda(t_0-t)))$$

from which it follows that $\lim_{t\to-\infty} p_2(t) =0$.

From these estimates on $p_1(t)$, $p_2(t)$ and $p_3(t)$, we see that by taking $t\to -\infty$ in (<ref>), we obtain that for all $t_0\in \mathbb{R}$

\begin{equation}
\label{E:mon1}
\int v^2 ( y,t_0) ( \phi (y- y_0) -  \phi (-y_0)  ) dy \lesssim y_0^{-6/5}
\end{equation}

The whole argument leading to (<ref>) applies with $v(y,t)$ replaced by $v(-y,-t)$, so that we can also assert that (<ref>) holds with $v(y,t)$ replaced by $v(-y,-t)$.  Thus, for all $t_1\in \mathbb{R}$

\begin{equation}
\label{E:mon2}
\int v^2 ( -y,-t_1) ( \phi (y- y_0) -  \phi (-y_0)  ) dy \lesssim y_0^{-6/5}
\end{equation}

Changing variable $-y \mapsto y$, and using that $\phi(-y-y_0)-\phi(-y_0) = \phi(y_0) - \phi(y+y_0)$ (which follows from the formula (<ref>) for $\phi$), we have

\begin{equation}
\label{E:mon3}
\int v^2 ( y, -t_1) ( \phi (y_0) -  \phi (y+y_0)  ) dy \lesssim y_0^{-6/5}
\end{equation}

Taking $t_1=-t_0$, and adding (<ref>) and (<ref>), we obtain

\begin{equation}
\label{E:mon4}
\int v^2 ( y, t_0) \rho(y,y_0) dy \lesssim y_0^{-6/5}
\end{equation}

$$\rho(y,y_0) = \phi(y-y_0)- \phi(-y_0) - \phi(y+y_0) + \phi(y_0)$$

From the formula (<ref>) for $\phi$, it follows that $\rho$ is even in $y$ (that is $\rho(-y,y_0) = \rho(y,y_0)$) and $\partial_y \rho(y,y_0) \geq 0$ for $y>0$.  Since $\rho(0,y_0)=0$ and $\rho(y_0,y_0) \geq \frac{\pi}{6}$ whenever $y_0\geq \sqrt{3}A$, it follows that $\rho(y,y_0) \geq 0$ for all $y\in \mathbb{R}$ and $\rho(y,y_0) \geq \frac{\pi}{6}$ when $|y| \geq y_0$ (provided $y_0 \geq \sqrt{3}A$).  Thus from (<ref>)

$$\forall \; y_0 > \sqrt{3}A, \quad \int_{|y| > y_0} v^2 (t_0, y) dy \lesssim y_0^{-6/5}

from which we can integrate in $y_0$ and find that uniformly for all $t\in \mathbb{R}$

\begin{equation}
\label{E:finitevar}
\int_{y\in \mathbb{R}} |y| v^2 (t, y) dy \lesssim 1
\end{equation}

The (nonlocalized) virial identity obtained by computing $\partial_t \int y v(y,t)^2 \, dy$, substituting the equation (<ref>) for $v$, applying integration by parts in $y$, and integrating over $t_1\leq t \leq t_2$, is

$$ \int y v^2 (t_2) \, dy - \int y v^2 (t_1) \,dy  =  - \|v\|_{L_y^2}^2 - 2\|D_y^{1/2} v\|_{L_y^2}^2 + \int (Q  - yQ') v^2 \, dy$$

From this, it follows that

$$\| v \|_{L_{[t_1,t_2]}^2 H_y^{1/2}}^2 \lesssim \| |y|^{1/2} v \|_{L_{[t_1,t_2]}^\infty L_y^2}^2 + \| \la y \ra^{-1} v \|_{L_{[t_1,t_2]}^2L_y^2}^2$$

By (<ref>) and (<ref>), the right side is bounded uniformly for all $t_1<t_2$, so taking $t_1=0$ and $t_2\to +\infty$ implies $\| v \|_{L_{t>0}^2 H_y^{1/2}}<\infty$.  Hence, there exists a time sequence $t_n \rightarrow \infty$ along which  $ \|v(t_n) \|_{H_y^{1/2}} \rightarrow 0$ as $n \rightarrow +\infty$.
Now, from the fact that $\mathcal{L} Q' =0$ we can deduce that the quantity $ \langle \mathcal{L} v(t), v(t) \rangle$ is conserved in time.   Hence for any $t$,

\langle \mathcal{L} v(t), v(t)\rangle = \lim_{t_n \rightarrow +\infty} \langle \mathcal{L} v(t_n), v(t_n) \rangle =0

But since $v$ satisfies the orthogonality conditions (<ref>), it follows that $\langle \mathcal{L} v(t), v(t)\rangle  \gtrsim \|v(t) \|^2_{H^{1/2}} $.  Therefore, $v(t)\equiv 0$ for all $t$, as claimed.

[1] C.J.  Amick and J.F. Toland, Uniqueness and related analytic properties for the Benjamin-Ono equation–a nonlinear Neumann problem in the plane,  Acta Math. 167 (1991), no. 1-2, 107–126.

[2] T. B. Benjamin, Internal waves of permanent form in fluids of great depth, J. Fluid Mech. 29 (1967), pp. 559–562.

[3] D.P. Bennett, R.W. Brown, S.E. Stansfield, J.D. Stroughair and J.L. Bona, The stability of internal solitary waves, Mathematical Proceedings of the Cambridge Philosophical Society, 94 (1983), pp. 351–379.

[4] T.L. Bock and M.D. Kruskal, A two-parameter Miura transformation of the Benjamin-Ono equation, Volume 74A, Number 3-4, Physics Letters, 12 November 1979, pp. 173–176.

[5]  A.-P. Calderón, Commutators of singular integral operators, Proc. Nat. Acad. Sci. U.S.A. 53 (1965), pp. 1092–1099.

[6] K. M. Case, The $N$-soliton solution of the Benjamin-Ono equation, Proc. Nati. Acad. Sci. USA, Vol. 75, No. 8,  August 1978, pp. 3562–3563.

[7] Yu Deng, Invariance of the Gibbs measure for the Benjamin-Ono equation, J. Eur. Math. Soc. (JEMS) 17 (2015), no. 5, pp. 1107–1198.

[8] Yu Deng, Nikolay Tzvetkov, Nicola Visciglia, Invariant measures and long time behaviour for the Benjamin-Ono equation III, Comm. Math. Phys. 339 (2015), no. 3, pp. 815–857.

[9] Luiz Gustavo Farah, Justin Holmer, Svetlana Roudenko, Kai Yang, Asymptotic stability of solitary waves of the 3D quadratic Zakharov-Kuznetsov equation, .

[10] Luiz Gustavo Farah, Justin Holmer, Svetlana Roudenko, Kai Yang, Blow-up in finite or infinite time of the 2D cubic Zakharov-Kuznetsov equation, .

[11] R. Frank and E. Lenzmann, Uniqueness of non-linear ground states for fractional Laplacians in $\mathbb{R}$, Acta Math. 210 (2013), no. 2, pp. 261–318.

[12] Patrick Gérard, Thomas Kappeler, Petar Topalov, On the spectrum of the Lax operator of the Benjamin-Ono equation on the torus,  J. Funct. Anal. 279 (2020), no. 12, 108762, 75 pp.

[13] J. Ginibre and G. Velo, Proprietes de lissage et existence de solutions pour l'equation de Benjamin-Ono generalisee, C.R. Acad. Sci. Paris 308 (1989), pp. 309–314. Zbl 0694.35156

[14] Stephen Gustafson, Hideo Takaoka,  Tai-Peng Tsai, Stability in $H^{1/2}$ of the sum of $K$ solitons for the Benjamin-Ono equation, (English summary)
J. Math. Phys. 50 (2009), no. 1, 013101, 14 pp.

[15] Nakao Hayashi, Elena I. Kaikina, Neumann initial-boundary value problem for Benjamin-Ono equation, J. Differential Equations 252 (2012), no. 5, pp. 3520–-3540.

[16] Mihaela Ifrim and Daniel Tataru, Well-posedness and dispersive decay of small data solutions for the Benjamin-Ono equation, Ann. Sci. Éc. Norm. Supér. (4) 52 (2019), no. 2, pp. 297–335.

[17] A.D. Ionescu and C.E.  Kenig, Global well-posedness of the Benjamin-Ono equation in low-regularity spaces, J. Amer. Math. Soc. 20 (2007), no. 3, pp. 753–798.

[18] R.J. Iorio Jr., On the Cauchy problem for the Benjamin-Ono equation, Comm. Partial Diff. Eq. 11 (1986), pp. 1031–1181.

[19]  H. Koch and N. Tzvetkov, On the local well-posedness of the Benjamin-Ono equation in $H^s(\mathbb{R})$, Int. Math. Res. Not. 2003 (2003), no. 26, pp. 1449–1464.

[20] Carlos E. Kenig and Yvan Martel, Asymptotic stability of solitons for the Benjamin-Ono equation, Rev. Mat. Iberoam. 25 (2009), no. 3, pp. 909–970.

[21] C.E. Kenig and K. D. Koenig, On the local well-posedness of the Benjamin-Ono and modified Benjamin-Ono equations, Math. Res. Lett. 10 (2003), no. 5-6, pp. 879–895.

[22] Carlos Kenig, Gustavo Ponce, and Luis Vega, Well-posedness and scattering results for the generalized Korteweg-de Vries equation via the contraction principle, Comm. Pure Appl. Math. 46 (1993), no. 4, pp. 527–620.

[23] C. E. Kenig, G. Ponce, and L. Vega, Uniqueness properties of solutions to the Benjamin-Ono equation and related models, J. Funct. Anal. 278 (2020), no. 5, 108396, 14 pp.

[24] Kihyun Kim and Soonsik Kwon, Scattering for defocusing generalized Benjamin-Ono equation in the energy space $H^{1/2}(\mathbb{R})$, Trans. Amer. Math. Soc. 372 (2019), no. 7, pp. 5011–5067.

[25] Camille Laurent, Felipe Linares, and Lionel Rosier, Control and stabilization of the Benjamin-Ono equation in $L^2(\mathbb{T})$, Arch. Ration. Mech. Anal. 218 (2015), no. 3, pp. 1531–1575.

[26] F. Linares, A. Mendez, G. Ponce, Asymptotic Behavior of Solutions of the Dispersion Generalized Benjamin-Ono Equation, J. Dynam. Differential Equations 33 (2021), no. 2, pp. 971–984.

[27] Y. Martel, Linear problems related to asymptotic stability of solitons of the generalized KdV equations SIAM J. Math. Anal. 38 (2006), pp. 759–781.

[28] Y. Martel and F. Merle, Liouville theorem for the critical generalized Korteweg-de Vries equation, J. Math. Pures Appl. (9) 79 (2000), no. 4, pp. 339–425.

[29]  Yvan Martel and Didier Pilod Construction of a minimal mass blow up solution of the modified Benjamin-Ono equation, Math. Ann. 369 (2017), no. 1-2, pp. 153–245.

[30] Y. Matsuno, Higher-order nonlinear evolution equation for interfacial waves in a two-layer fluid system, Physical Review E Vol. 49, pp. 2091 (1994).

[31] Y. Matsuno, Multisoliton perturbation theory for the Benjamin-Ono equation and its application to real physical systems, Physical Review E Vol. 51, No. 2, pp. 1471–1483 (1995)

[32] Peter D. Miller and Alfredo N. Wetzel, Direct scattering for the Benjamin-Ono equation with rational initial data, Stud. Appl. Math. 137 (2016), no. 1, pp. 53–69.

[33] Peter D. Miller, Alfredo N. Wetzel, The scattering transform for the Benjamin-Ono equation in the small-dispersion limit, Phys. D 333 (2016), pp. 185–199.

[34] L. Molinet and D. Pilod, The Cauchy problem for the Benjamin-Ono equation in $L^2$ revisited, Anal. PDE 5 (2012), no. 2, pp. 365–395.

[35] Claudio Muñoz and Gustavo Ponce, On the asymptotic behavior of solutions to the Benjamin-Ono equation, Proc. Amer. Math. Soc. 147 (2019), no. 12, pp. 5303–5312.

[36] H. Ono, Algebraic solitary waves in stratified fluids, J. Phys.Soc. Japan 39 (1975), pp. 1082–1091.

[37] G. Ponce, On the global well-posedness of the Benjamin-Ono equation, Diff. Int. Eq. 4 (1991), pp. 527–542.

[38] Oscar Riaño, Svetlana Roudenko, Kai Yang, Higher dimensional generalization of the Benjamin-Ono equation: 2D case, 	.

[39] Svetlana Roudenko, Zhongming Wang, Kai Yang, Dynamics of solutions in the generalized Benjamin-Ono equation: a numerical study,

[40] J.C. Saut, Sur quelques generalizations de l'equation de Korteweg-de Vries, J. Math. Pures Appl. 58 (1979), pp. 21–61.

[41] Jean-Claude Saut, Benjamin-Ono and intermediate long wave equations: modeling, IST and PDE,  Nonlinear dispersive partial differential equations and inverse scattering, 95–160, Fields Inst. Commun., 83, Springer, New York, 2019.

[42]  Elias M. Stein and Rami Shakarchi, Real analysis. Measure theory, integration, and Hilbert spaces, Princeton Lectures in Analysis, 3. Princeton University Press, Princeton, NJ, 2005. xx+402 pp. ISBN: 0-691-11386-6.

[43] Mouhamadou Sy, Invariant measure and long time behavior of regular solutions of the Benjamin-Ono equation, Anal. PDE 11 (2018), no. 8, pp. 1841–1879.

[44] Terence Tao, Global well-posedness of the Benjamin-Ono equation in $H^1(\mathbb{R})$, J. Hyperbolic Differ. Equ. 1 (2004), no. 1, pp. 27–49.

[45] Nikolay Tzvetkov and Nicola Visciglia, Gaussian measures associated to the higher order conservation laws of the Benjamin-Ono equation, Ann. Sci. Éc. Norm. Supér. (4) 46 (2013), no. 2, pp. 249–299 (2013).

[46] Nikolay Tzvetkov and Nicola Visciglia, Invariant measures and long-time behavior for the Benjamin-Ono equation, Int. Math. Res. Not. IMRN 2014, no. 17, pp. 4679–4714.

[47] Nikolay Tzvetkov and Nicola Visciglia, Invariant measures and long time behaviour for the Benjamin-Ono equation II, J. Math. Pures Appl. (9) 103 (2015), no. 1, pp. 102–141.

[48] Yilun Wu, Jost solutions and the direct scattering problem of the Benjamin-Ono equation, SIAM J. Math. Anal. 49 (2017), no. 6, pp. 5158–5206.

[49] Yilun Wu, Simplicity and finiteness of discrete spectrum of the Benjamin-Ono scattering operator, SIAM J. Math. Anal. 48 (2016), no. 2, pp. 1348–1367.

[50] Katherine Zhiyuan Zhang, Benjamin-Ono soliton dynamics in a slowly varying potential,  Nonlinearity 33 (2020), no. 3, pp. 1064–1093.


# Atmospheric Response for MeV $\mathrm{\gamma}$ Rays Observed with Balloon-
Borne Detectors

Christopher M. Karwin NASA Postdoctoral Program Fellow
<EMAIL_ADDRESS>NASA Goddard Space Flight Center, Greenbelt, MD,
20771, USA  Carolyn Kierans NASA Goddard Space Flight Center, Greenbelt, MD,
20771, USA  Albert Y. Shih NASA Goddard Space Flight Center, Greenbelt, MD,
20771, USA  Israel Martinez Castellanos NASA Goddard Space Flight Center,
Greenbelt, MD, 20771, USA  Department of Astronomy, University of Maryland,
College Park, Maryland 20742, USA  Alex Lowell  Space Sciences Laboratory, UC
Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA  Thomas Siegert Julius-
Maximilians-Universität Würzburg, Fakultät für Physik und Astronomie, Institut
für Theoretische Physik und Astrophysik, Lehrstuhl für Astronomie, Emil-
Fischer-Str. 31, D-97074 Würzburg, Germany  Jarred Roberts  Department of
Astronomy & Astrophysics, UC San Diego, 9500 Gilman Drive, La Jolla CA 92093,
USA  Savitri Gallego Institut für Physik & Exzellenzcluster PRISMA+, Johannes
Gutenberg-Universität Mainz, 55099 Mainz, Germany  Adrien Laviron Laboratoire
Leprince-Ringuet, CNRS/IN2P3, École polytechnique, Institut Polytechnique de
Paris, 91120 Palaiseau, France  Andreas Zoglauer Space Sciences Laboratory, UC
Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA  John A. Tomsick Space Sciences
Laboratory, UC Berkeley, 7 Gauss Way, Berkeley, CA 94720, USA  Steven E. Boggs
Department of Astronomy & Astrophysics, UC San Diego, 9500 Gilman Drive, La
Jolla CA 92093, USA

###### Abstract

The atmospheric response for MeV $\gamma$ rays ($\sim$ 0.1 $-$ 10 MeV) can be
characterized in terms of two observed components. The first component is due
to photons that reach the detector without scattering. The second component is
due to photons that reach the detector after scattering one or more times.
While the former can be determined in a straightforward manner, the latter is
much more complex to quantify, as it requires tracking the transport of all
source photons that are incident on Earth’s atmosphere. The scattered
component can cause a significant energy-dependent distortion in the measured
spectrum, which is important to account for when making balloon-borne
observations. In this work we simulate the full response for $\gamma$-ray
transport in the atmosphere. We find that the scattered component becomes
increasingly more significant towards lower energies, and at 0.1 MeV it may
increase the measured flux by as much as a factor of $\sim 2-4$, depending on
the photon index and off-axis angle of the source. This is particularly
important for diffuse sources, whereas the effect from scattering can be
significantly reduced for point sources observed with an imaging telescope.

##  1 Introduction

As $\gamma$ rays with energies between $\sim$ 0.1 $-$ 10 MeV travel through
Earth’s atmosphere, they may undergo Compton scattering. This causes
attenuation and distortion of the original signal. In order to overcome these
atmospheric effects, observations in this energy band are typically made with
space-based telescopes. However, balloon-borne observations within Earth’s
atmosphere are still essential for the development of new telescope
technologies. Such development is quickly progressing, as exemplified by the
long-duration balloon flights of the Compton Spectrometer and Imager (COSI)
(Kierans et al., 2016; Tomsick et al., 2019; Beechert et al., 2022; Tomsick et
al., 2023) and the Sub-MeV/MeV gamma-ray Imaging Loaded-on-balloon Experiments
(SMILE) (Takada et al., 2022), as well as the recent short-duration balloon
flight of the Compton Pair telescope (ComPair) (Valverde et al., 2023), and
the upcoming balloon flight of the Antarctic Demonstrator for the Advanced
Particle-astrophysics Telescope (ADAPT) (Chen & Buckley, 2023). Moreover, a
number of other mission concepts have recently been proposed, and may very
well undergo balloon tests in the coming years, e.g., the All-sky Medium
Energy Gamma-ray Observatory eXplorer (AMEGO-X) (Caputo et al., 2022) and the
Galactic Explorer with a Coded aperture mask Compton telescope (GECCO)
(Orlando et al., 2022). In order to make an accurate assessment of a
telescope’s performance when operating at balloon altitudes, it is imperative
to have a complete understanding of how the atmosphere affects $\gamma$-ray
transport. Specifically, in this work we focus on the spectral response.

The scattering of MeV photons in the atmosphere can be characterized in terms
of two components. First, only a fraction of photons from a source will travel
through the atmosphere and reach the detector without scattering, which we
refer to as the transmitted photons. This component can be accounted for in a
straightforward manner by calculating the corresponding transmission
probability (TP). Generally, the TP depends on the initial energy of the
photon, the off-axis angle of the source, and the altitude of the
observations. Second, some fraction of photons from a source will reach the
detector after one or more scatters (even if not initially directed towards
the detector), which we refer to as the scattered photons. Accounting for this
component is much more challenging, as it requires tracking the $\gamma$-ray
transport in the atmosphere for all incident photons. Moreover, $\gamma$ rays
from astrophysical sources pass the Earth as plane waves, and in principle
this implies that the scattered photons can come from a surface area
effectively as large as the cross-section of Earth’s upper atmosphere, which
far exceeds the area of a detector. This can lead to difficulties when it
comes to acquiring appropriate statistics from Monte Carlo simulations of a
detector’s response.

The scattered component is important to account for when analyzing balloon-
borne observations, as it can produce a significant energy-dependent
distortion in the measured spectrum. In general, this tends to lead to more
photons towards lower energies. For diffuse continuum sources, such as the
Galactic diffuse continuum emission and extragalactic $\gamma$-ray background,
photons enter the detector from all directions, and therefore the spectral
distortion from the scattered photons becomes very significant. The effect is
not as crucial for point sources observed with imaging telescopes, such as
COSI, because they have the ability to decipher the photon’s incident
direction. However, the degree to which photons can be selected coming from
the source direction also has a strong dependence on the instrument’s angular
resolution.

Figure 1: Left: Schematic showing the simulation setup (not to scale). We use
a spherically symmetric mass model of the atmosphere, defined relative to
Earth’s surface. The atmosphere is characterized using NRLMSIS, corresponding
to the plot on the right. See text for more details. Right: Atmospheric
profile defining the mass model. The left axis shows the number density for
the primary species of the atmosphere, nitrogen (N), oxygen (O), argon (Ar),
helium (He), and hydrogen (H). The right axis shows the total mass density.

In this work we determine the full response for $\gamma$-ray transport in the
atmosphere via Monte Carlo simulations, which includes both the transmitted
and scattered components. As part of this, we publicly release the COSI
atmosphere simulation and analysis pipeline, cosi-atmosphere111The cosi-
atmosphere package is available at https://cosi-
atmosphere.readthedocs.io/en/latest/. Additionally, we provide atmospheric
response matrices calculated for altitudes between 25 $-$ 45 km (in 0.5 km
steps)222Instructions for accessing the atmospheric response matrices can be
found in the cosi-atmosphere documentation.. The goal of the cosi-atmosphere
package is to build a user-friendly Python-based library that can be employed
for atmospheric physics associated with MeV $\gamma$-ray astronomy. Note that
these tools are independent of any specific detector, and thus they can be
easily adapted for different observations. The package currently includes the
spectral response for balloon-borne observations, and we plan to extend the
tools to other relevant topics in the near future, including the determination
of atmospheric $\gamma$-ray backgrounds, which dominant the emission at
balloon-altitudes, as well as determination of the $\gamma$-ray albedo and
reflection components, for space-based observations.

The paper proceeds as follows. In Section 2 we detail the simulation setup,
which is based on a spherical geometry. In Section 3 we describe the
atmospheric response, including some specific applications, and validation of
the simulations. The summary and conclusions are given in Section 4. In
Appendix A we give more details about the simulations. In Appendix B we give
details for calculating the TP analytically. Appendix C provides an additional
example complementing Section 3. In Appendix D we present response
calculations using a simplified model of the atmosphere, consisting of a
rectangular geometry, which we show to be consistent with the spherical
geometry simulations. In order to use a concrete example, our calculations in
this work are based on the 2016 COSI balloon flight (Kierans et al., 2016),
which primarily motivates our choice of inputs for the atmospheric model.
Specifically, we use a representative date and geographical location of
2016-06-13 and $(\mathrm{lat,lon})=(-5.66^{\circ},-107.38^{\circ})$,
respectively.

##  2 Simulations

###  2.1 Atmospheric Mass Model

The atmospheric response is simulated with the COSI atmosphere pipeline, which
employs the Medium-Energy Gamma-ray Astronomy library (MEGAlib) software
package333MEGAlib is available at https://megalibtoolkit.com/home.html
(Zoglauer et al., 2006), based on Geant4 (Agostinelli et al., 2003). We create
a spherically symmetric mass model of Earth’s atmosphere, as depicted in the
left panel of Figure 1. The model is comprised of spherical atmospheric shells
having a thickness of 100 m, and extending from Earth’s surface444We use
Earth’s equatorial radius, which is slightly larger than the polar radius of
6357 km. ($R_{\earth}=6378$ km) to an altitude of 200 km. The atmosphere is
characterized using the latest version (v2.1) of the Naval Research
Laboratory’s Mass Spectrometer Incoherent Scatter Radar Model
(NRLMSIS)555NRLMSIS is available at https://swx-trec.com/msis (Picone et al.,
2002; Emmert et al., 2021, 2022), implemented in the COSI atmosphere pipeline
via the Python interface, pymsis666pymsis is available at
https://swxtrec.github.io/pymsis/. NRLMSIS is an empirical model of Earth’s
atmosphere that describes the average observed behavior of temperature and
density, from the ground to an altitude of roughly 1000 km. More specifically,
the model specifies the altitude profile of the number density for the primary
species of the atmosphere (i.e., nitrogen, oxygen, argon, and helium), as
shown in the right panel of Figure 1. In general, this is dependent on
location (latitude, longitude, altitude), time (year and day), and solar and
geomagnetic activity levels.

###  2.2 Photon Tracking

We simulate an isotropic source with a flat energy spectrum (i.e., constant
number of photons per energy bin) between 10 keV $-$ 10 MeV, using $10^{7}$
photons. The source is simulated using a surrounding sphere with a radius of
$R_{\earth}+200$ km. In order to mimic plane waves incident on the atmosphere,
photons are emitted perpendicularly from a disk located at the surrounding
sphere, and having the same radius (see left panel of Figure 1). With the
isotropic source, photons are directed towards Earth’s surface uniformly over
the sky at all incident angles. The electromagnetic processes that occur as
the $\gamma$ rays pass through the atmosphere are described with the Geant4
Livermore physics models, which includes Compton scattering, pair conversion,
photoelectric absorption, and Rayleigh scattering. Correspondingly, secondary
$\gamma$ rays may also be produced from pair production and subsequent
annihilation in the atmosphere. The simulations track the $\gamma$-ray
transport, including initial ($i$) and “measured” ($m$) values of the photon’s
energy ($E$), position ($\vec{r}$), and direction ($\hat{d}$), as indicated in
the left panel of Figure 1. The “measured” values are obtained by tracking the
photons whenever they cross a watched777In this context, “watched” just refers
to the volume being used to monitor the photon properties. volume, consisting
of a spherical shell at a radius of 33.5 km. With this method, our
calculations are completely independent of any specific detector. It should be
noted that photons in the simulations may cross the watched volume numerous
times, and for each crossing the tracked information is stored. In general
though, each successive crossing is increasingly unlikely, as most photons
that go to lower altitudes do not back scatter, and the ones that do scatter
back above the balloon altitude are highly unlikely to scatter back again.
There are certainly cases where analyzing multiple crossings for a single
photon becomes relevant. However, in this work we only consider the first
crossing, which is a reasonable simplification considering that an Earth limb
cut of $\sim 90^{\circ}$ is typically employed in the real data analysis.

###  2.3 Determination of Incident Angles

From the tracked values, the incident angle ($\theta$) relative to the surface
normal ($\hat{n}$) can be determined, for both initial and measured photons.
The normal vector at any point on a sphere can be calculated from the position
vector at that point, i.e., $\hat{n}=\vec{r}/\lVert\vec{r}\rVert$. The
incident angle is then given by the standard equation:

$\theta=\mathrm{cos}^{-1}(-\hat{n}\cdot\hat{d}),$ (1)

where both $\hat{n}$ and $\hat{d}$ are unit vectors. The normal vector for
each event is determined from the position of the measured photon, and this is
also used for determining the initial incident angle. In this scheme,
$\theta_{i}$ is undefined for photons that aren’t measured. This occurs in two
situations: 1) Photons are not emitted towards the watched volume, and do not
get scattered into it. 2) Photons are initially emitted towards the watched
volume, but get scattered away from it. While photons in the first case can be
disregarded, photons in the second case must be accounted for in order to get
the normalization of the response correct. The reason for this is that for any
given initial incident angle, we must count the photons that get scattered and
never reach the surface.

In order to determine the initial incident angle for unmeasured photons, we
can find where the trajectory defined by the initial photon parameters
intersects the watched volume. This will give us the normal vector that can be
used in Eq. 1. The standard vector equation of a sphere centered at the origin
with radius $r$ is given by

$||\vec{x}_{s}||^{2}=r^{2},$ (2)

where $\vec{x}_{s}$ are points on the sphere. The standard vector equation of
a line is given by

$\vec{x}_{l}=\vec{o}+D\hat{u},$ (3)

where $\vec{o}$ is the origin of the line (which in our case corresponds to
$\vec{r}_{i}$), $D$ is the distance from the origin, $\hat{u}$ is a unit
vector giving the direction of the line (corresponding to $\hat{d}_{i}$), and
$\vec{x}_{l}$ gives points on the line. Plugging Eq. 3 into Eq. 2 allows us to
solve for $D$:

$D=\frac{(-\hat{u}\cdot\vec{o})\pm\sqrt{(\hat{u}\cdot\vec{o})^{2}-||\hat{u}||^{2}(||\vec{o}||^{2}-r^{2})}}{||\hat{u}||^{2}},$
(4)

where we take the shortest distance because we want the first intersection.
Once the intersection is known, we can determine the initial incident angle.
From our $10^{7}$ simulated photons, $1.5\times 10^{6}$ (15%) are not
measured, and $4.7\times 10^{5}$ (4.7%) have no (real) solution to Eq. 4
(meaning that they were not initially directed towards the watched volume).
The total number of initially thrown photons used in the response
normalization is determined by the latter, which gives $9.53\times 10^{6}$.

###  2.4 Geometric Correction Factor

The goal of our simulations is to characterize the atmospheric response
independent of any assumption on the detector geometry. Therefore, we need to
apply a geometric correction factor,
$\zeta_{\mathrm{g}}(\theta_{i},\theta_{m})$, to the measured photons, in order
to account for the 2-dimensional detecting area used in the simulations. Such
a factor has also been employed in Pendleton et al. (1995), and here we follow
a similar approach. Because the curvature of Earth is very large, we can
approximate it as being locally flat over the general area where we expect a
majority of the scattered photons to travel from. Under this assumption, for a
flux, $F$, from a direction, $\theta$, the counts (N) through a horizontal
patch of area ($a$) of the watched surface has a cos($\theta$) factor due to
the projection:

$N=F(\theta)\times t\times a\times\mathrm{cos}(\theta),$ (5)

where $t$ is the exposure time. Accordingly, the ratio of measured flux
($F_{m}$) to initial flux ($F_{i}$) for respective directions $\theta_{m}$ and
$\theta_{i}$ is:

$\begin{split}\frac{F_{m}(\theta_{m})}{F_{i}(\theta_{i})}&=\frac{N_{m}}{a\mathrm{cos}(\theta_{m})t}\times\frac{a\mathrm{cos}(\theta_{i})t}{N_{i}}\\\
&=\frac{N_{m}}{N_{i}}\times\frac{\mathrm{cos}(\theta_{i})}{\mathrm{cos}(\theta_{m})}.\end{split}$
(6)

In order to properly normalize the response, we must account for the
difference in the projected area. Thus, this implies that the geometric
correction factor is given by:

$\zeta_{\mathrm{g}}(\theta_{i},\theta_{m})=\frac{\mathrm{cos}(\theta_{i})}{\mathrm{cos}(\theta_{m})}.$
(7)

When constructing the atmospheric response matrices, the number of observed
counts, $N_{m}(\theta_{i},\theta_{m})$, is weighted by this factor. We have
verified that the distributions of $\vec{r}$, $\hat{d}$, $E$, and $\theta$
obtained from the simulations are consistent with the expectations from the
simulation setup, as discussed in Appendix A.

##  3 Atmospheric Response

###  3.1 Characterizing the Response

For a given altitude ($A$), the atmospheric response
($\epsilon_{\mathrm{atm}}$) can be quantified in terms of initial and measured
values of energy, incident angle, and azimuth angle, corresponding to the six
parameters $E_{i},\ E_{m},\ \theta_{i},\ \theta_{m},\ \phi_{i},\ \mathrm{and}\
\phi_{m}$. In this work we consider two specific representations of the
response. For the first representation, we define an azimuth angle
($\phi_{1}$) with respect to the zenith of the detector, relative to the
photon’s initial position. Specifically, we use the difference between the
azimuth angles for the vectors $\hat{d_{m}}$ and $\hat{d_{i}}$, which in
spherical coordinates can be expressed as

$\phi_{1}=\phi_{dm}-\phi_{di},$ (8)

where $\phi_{dm}$ is the azimuth of $\hat{d_{m}}$, and $\phi_{di}$ is the
azimuth of $\hat{d_{i}}$. Thus, for this first case, the response is
represented in terms of $E_{i},\ E_{m},\ \theta_{i},\ \theta_{m},\
\mathrm{and}\ \phi_{1}$. This representation is particularly useful when
dealing with diffuse sources, such as the Galactic diffuse or extragalactic
gamma-ray background. In these cases, photons enter the detector from all
directions, and typically we can simplify the response by summing over all
azimuth angles. Indeed, one of the primary motivations for this work is to
quantify the atmospheric response for the scattered component, which is most
important when analyzing diffuse sources. This representation is also
applicable when dealing with non-imaging telescopes.

For the second representation, we characterize the response with respect to
the position of the source on the sky. To accomplish this, we use the
difference of incident directions:

$\Delta\theta_{\mathrm{s}}=\mathrm{cos}^{-1}(\hat{d_{i}}\cdot\hat{d_{m}}).$
(9)

Additionally, we define an azimuth angle ($\phi_{2}$) in the perpendicular
plane to the source position. The normal to the plane is given by
$\hat{d}_{i}$, and the initial point in the plane is obtained by the
orthogonal component ($\hat{z}_{\perp}$) of the projection of $\hat{z}$ onto
$\hat{d_{i}}$, where $\hat{z}$ is the zenith of the instrument. For the
initial direction ($\hat{\phi}_{2}$) we use the cross product:

$\hat{\phi}_{2}=\hat{d_{i}}\times\hat{z}_{\perp}.$ (10)

Thus, for this second case, the response is represented in terms of $E_{i},\
E_{m},\ \theta_{i},\ \Delta\theta_{\mathrm{s}},\ \mathrm{and}\ \phi_{2}$. This
representation is useful for analyzing point sources, observed with an imaging
telescope, such as COSI. In this case, photons can be selected coming from the
direction of the source, which is primarily limited by the angular resolution
of the instrument. This is discussed further in Section 3.5.

###  3.2 Convolving the Response

Figure 2: Energy dispersion matrix for all photons (left), transmitted photons
(middle), and scattered photons (right), for a 50∘ off-axis source. The y-axis
is the measured photon energy, and the x-axis is the initial photon energy.
The matrices are normalized by the total photons simulated in each bin of
$E_{i}$ and $\theta_{i}$. The total energy dispersion matrix is the sum of the
transmitted and scattered components (note that the components are shown with
different colorbar ranges).

In order to completely account for the atmospheric response, a full
convolution with the detector response ($R$) must be made. Considering our
first representation of the atmospheric response at a given altitude, the
predicted counts in energy bin $i$, for an off-axis angle $\theta_{i}$888In
this work we consider an instrument with a zenith pointing, in which case
“off-axis” and “zenith” are synonymous. and an exposure time $\Delta t$, can
be calculated as

$\begin{split}P_{i}=\Delta
t\iiint_{\phi_{1},\theta_{m},E_{m}}\epsilon_{\mathrm{atm}}(A,\theta_{i},\theta_{m},E_{i},E_{m},\phi_{1})\\\
\circledast R(Z,\phi,E_{i},E_{m})\circledast F(E_{i},\vec{\alpha})\ dE_{m}\
d\theta_{m}\ d\phi_{1},\end{split}$ (11)

where $Z$ and $\phi$ are the zenith and azimuth of the photon source,
respectively, $F$ is the input spectral model, and $\vec{\alpha}$ gives the
model parameters of the spectrum. Note that the convolution must map
$\theta_{m}$ to $Z$ and $\phi_{1}$ to $\phi$, in the working reference frame.
This convolution fully accounts for the fact that the scattered photons are
detected at different off-axis angles compared to the source, and thus
detected with a different detector response. A similar convolution would be
made for our second representation of the atmospheric response, except that we
would need to replace $\phi_{1}\rightarrow\phi_{2}$ and
$\theta_{m}\rightarrow\Delta\theta$.

###  3.3 Energy Dispersion Matrices

In order to examine some basic properties of the atmospheric response for both
the transmitted and scattered components, in this section we calculate energy
dispersion matrices. As our representative case we consider a source with an
off-axis angle of $50^{\circ}$, and we project the response onto the axes
$E_{i}$ and $E_{m}$. We bin the simulations using $4^{\circ}$ angular bins and
16 log-spaced energy bins ($\sim$ 5.33 bins per decade). This particular
binning is primarily motivated by the need to obtain sufficient statistics in
each bin. Additional photons can always be simulated if finer bins are
required, such as for the detailed analysis of spectral lines. Note that the
angular binning is comparable to the resolution of the COSI balloon
instrument: $5.9^{\circ}$ at 0.511 MeV; $3.9^{\circ}$ at 1.809 MeV (Beechert
et al., 2022).

Figure 2 shows the energy dispersion matrices for both the scattered and
transmitted photons, as well as the sum of the two components. We identify a
photon as having undergone at least one scatter if its measured incident angle
varies from its initial incident angle. Of course with this approach we are
limited by the angular resolution of our incident angle bins, i.e., we cannot
resolve a photon that scatters within $4^{\circ}$. The matrices are normalized
by the total number of photons simulated in each bin of $E_{i}$ and
$\theta_{i}$, as discussed in Section 2.3. The color scale gives the ratio of
photons that are detected with a measured energy, $E_{m}$, given an initial
energy, $E_{i}$. For the transmitted component, the ratio is a proper
probability. However, for the scattered component, the measured incident angle
is different than the initial incident angle, and thus the ratio is not a true
probability (as it is not normalized to unity).

The transmitted photons are shown in the middle panel of Figure 2. As can be
seen, this component resides completely along the main diagonal, with a
probability that steadily increases towards the upper energy bound. The
photons are along the main diagonal because they are the ones that reach the
detector without scattering, and therefore there is no energy loss.

The scattered photons are shown in the right panel of Figure 2. The
distribution of the scattered component is substantially different than that
of the transmitted component. Most notably, there is a significant number of
off-diagonal photons, and the ratio is highest towards lower energies. A
diagonal component is still present, but this is mostly an artifact of the
coarse energy binning. Based on standard Compton dynamics, the photon’s energy
after a scattering ($E^{\prime}$) is a function of the scattering angle
$\theta$:

$E^{\prime}=\frac{E_{i}}{1+\frac{E_{i}}{m_{e}c^{2}}(1-\cos\theta)}~{},$ (12)

where $m_{e}c^{2}\simeq 0.511\,\mathrm{MeV}$ is the rest-mass energy of the
electron. For $E_{i}/{m_{e}c^{2}}\ll 1$, $E^{\prime}\sim E_{i}$ for all values
of $\theta$, leading to all scattered low-energy photons piling up close to
the diagonal. The scattered component also has a prominent feature in the
energy bin containing the 0.511 MeV line. This is due to pair production and
subsequent annihilation in the atmosphere. The electron-positron pairs are
produced from $\gamma$-ray photons interacting with the Coulomb field of
electrons or nuclei in the atmosphere.

The detection fraction can be obtained by projecting the energy dispersion
matrix onto the initial energy axis, as shown in Figure 3. Note that for the
transmitted component, the detection fraction is equivalent to the TP. These
results clearly show that the effect from the scattered component is most
dominant for initial energies between $\sim 0.2-0.4$ MeV. Moreover, the
detection fraction of the scattered component exceeds that of the transmitted
component for energies below $\sim$ 0.6 MeV. For comparison, we also show the
TP for the transmitted component calculated analytically. As can be seen, the
analytical calculation is in very good agreement with the simulations, which
is another important validation of our simulation pipeline. We note that the
results from the simulation appear to be slightly higher than the analytical
result, but this is due to the coarse binning in angle and energy (note that
data points are plotted at the geometric mean). Details of the analytical
calculation are provided in Appendix B.

Figure 3: Detection fraction for the three components, as specified in the
legend, for a 50∘ off-axis source. The detection fraction is the projection of
the energy dispersion matrix (see Figure 2) onto the initial energy axis, and
for the transmitted component it is equivalent to the transmission
probability.

###  3.4 Approximating the Atmospheric Response

Figure 4: Left: Correction factor (as defined in Eq. 14) for a range of power
law spectral models ranging from $1-3$, for a 50∘ off-axis source. The solid
line is for the transmitted component only (note that the models mostly
overlap), and the dashed line is for the total component, which includes
transmitted and scattered photons. The black dash-dot line is the TP, obtained
by projecting the energy dispersion matrix onto the initial energy axis.
Right: Correction factor ratio (as defined in Eq. 15) for the same range of
models shown in the left plot.

As discussed in Section 3.2, the proper way to correct for atmospheric effects
is to start with a full convolution of the atmospheric response with the
detector response. However, this requires detailed aspect information of the
instrument during the observations, and in general the calculation can be
complex. It would therefore be helpful to have a simplified approach that
could be used to estimate atmospheric effects. To this end, here we calculate
the ratio of the predicted counts to the model counts, for a given spectral
input, which we refer to as the correction factor. The predicted counts are
obtained by forward-folding the atmospheric response with the input spectrum,
where we use the energy dispersion matrices discussed in the previous section.
Thus, in this case we are assuming that the source is at an average off-axis
angle, and we are integrating over all azimuth angles. Moreover, we are
neglecting the fact that the scattered photons will be measured at different
incident angles, and thus with a different detector response. Despite these
simplifications, this approach provides a reasonable approximation of the
atmospheric effects for diffuse sources, such as the Galactic diffuse
continuum emission (Karwin et al., 2023).

As a toy example, we can consider a source with three generic energy bins. For
a given altitude and off-axis angle, the forward-folding is obtained as
follows:

$\epsilon_{\mathrm{atm}}\vec{F}=\vec{P},$ (13)

where $\epsilon_{\mathrm{atm}}$ is the energy dispersion matrix (as shown in
Figure 2), $\vec{F}$ is the model flux (i.e., the source flux incident on the
atmosphere), and $\vec{P}$ is the predicted flux (i.e., the resulting flux
incident on the instrument after accounting for atmospheric effects), both in
units of $\mathrm{ph\ cm^{-2}\ s^{-1}}$. Symbolically, the matrix
multiplication can be written out as follows:

$\begin{bmatrix}0&0&\epsilon_{33}\\\ 0&\epsilon_{22}&\epsilon_{23}\\\
\epsilon_{11}&\epsilon_{12}&\epsilon_{13}\end{bmatrix}\begin{pmatrix}F_{1}\\\
F_{2}\\\ F_{3}\end{pmatrix}=\begin{pmatrix}P_{3}\\\ P_{2}\\\
P_{1}\end{pmatrix},$

where the indices specify the energy channel. For the energy dispersion
matrix, the first index gives the row (corresponding to the measured energy
axis), and the second index gives the column (corresponding to the initial
energy axis). Note that the upper triangle of the energy dispersion matrix is
always zero, as the off-diagonal entries are found in the lower triangle, due
to photons which lose energy after scattering. From this we obtain

$\displaystyle P_{1}=$ $\displaystyle\
\epsilon_{11}F_{1}+\epsilon_{12}F_{2}+\epsilon_{13}F_{3}$ $\displaystyle
P_{2}=$ $\displaystyle\ \epsilon_{22}F_{2}+\epsilon_{23}F_{3}\ $
$\displaystyle P_{3}=$ $\displaystyle\ \epsilon_{33}F_{3}$

Writing out the equations in this way is helpful for gaining insight into the
dynamics of how the scattered component can produce an energy-dependent
distortion in the predicted spectrum. For each energy channel the predicted
counts is the sum of the flux from all channels at and above the given energy,
weighted by the respective ratios. In the case of the transmitted component,
the off-diagonal entries are all zero, and the predicted flux in a given
energy bin is a fraction of the model flux in that same bin (corresponding to
the TP). However, for the scattered component, the off-diagonal entries are
not all zero, and thus the predicted counts in a given energy bin will have
contributions from higher energy bins.

The energy-dependent correction factor, $c$, for a given energy bin, $E_{i}$,
can now be defined as

$c(E_{i})=P_{i}/F_{i}.$ (14)

The left panel of Figure 4 shows the correction factor for a range of power
law spectral models with different assumptions on the spectral index. The
correction factor is shown for both the transmitted component, as well as the
total component. For the former, the correction factor is equivalent to the
TP, and thus it is the same for all spectral indices. On the other hand, the
total correction factor has a dependence on the spectral index, due to the
scattered photons. Harder sources will generally have more photons at higher
energies that get scattered down to lower energies. Indeed, this is the exact
trend that we find.

As shown by the correction factor, the scattered component can have a
significant impact on the atmospheric response, compared to including only the
transmitted component. In order to quantify this further, we define the
correction factor ratio, $R(E_{i})$, as

$R(E_{i})=\frac{c_{\mathrm{tot}}(E_{i})}{c_{\mathrm{tran}}(E_{i})},$ (15)

where $c_{\mathrm{tot}}$ and $c_{\mathrm{tran}}$ are the correction factors
for the total and transmitted components, respectively. The right panel of
Figure 4 shows the correction factor ratio for the same range of power law
spectral models shown in the left panel. As can be seen, the scattered
component is most important towards lower energies. Note that photoelectric
absorption becomes dominant for energies below $\sim$ 100 keV, and the
atmosphere rapidly becomes opaque to $\gamma$-ray photons. In this regime, the
simulations have very low statistics, and the correction factor ratio becomes
very large. In this work we are mainly interested in energies $>0.1$ MeV, and
lower energies will be the focus of a future study.

For balloon-borne observations of a diffuse source at a given altitude, Eq. 14
can be applied to correct for the atmospheric response. This would be
equivalent to scaling the effective area. Alternatively, if the transmitted
component has already been accounted for (i.e., in the determination of the
effective area), Eq. 15 can be applied to correct for the scattered component.
However, we again stress that these corrections only provide approximations,
and in order to completely correct for the atmospheric response, a full
convolution with the detector response must be made, as described in Section
3.2.

So far we have only considered an off-axis angle of $50^{\circ}$. More
generally, a similar trend is found for all off-axis angles. For higher off-
axis angles, the detection fraction for both the transmitted and scattered
components decreases, although the correction factor ratio increases,
essentially due to the fact that there is more atmosphere for the photons to
traverse. As another example, in Appendix C we show the results for an on-axis
source.

###  3.5 Point Sources and Imaging Telescopes

Figure 5: Top: Distribution of $\Delta\theta_{\mathrm{s}}$ for an on-axis
source. Middle: 2-dimensional distribution of $\Delta\theta_{\mathrm{s}}$ and
initial energies, for measured energies between $0.750-1.155$ MeV. The image
is smoothed with a 2D Gaussian kernel, with standard deviations of 1 and 0.1
(pixels) in the x and y directions, respectively. The dashed green contour is
at the level of the max array value divided by 65. Bottom: Correction factor
ratio calculated for different max values of $\Delta\theta_{\mathrm{s}}$, as
specified in the legend.

For point sources observed with an imaging telescope, such as COSI, Eqs. 14
and 15 can also be applied as an approximation of the atmospheric response.
However, in this case, photons can be selected coming from the direction of
the source, mainly limited by the angular resolution of the instrument. The
approximation should therefore be made using the alternative representation of
the atmospheric response (in terms of $E_{i},\ E_{m},\ \theta_{i},\
\Delta\theta_{\mathrm{s}},\ \mathrm{and}\ \phi_{2}$). This allows us to
calculate the energy dispersion matrices using a slice of
$\Delta\theta_{\mathrm{s}}$, where we can choose a max value corresponding to
the instrument’s angular resolution.

As a specific example, here we consider an on-axis source
($\theta_{i}=0^{\circ}$), and we use the angular resolution of the COSI
balloon instrument ($5.9^{\circ}$ at 0.511 MeV). The top panel of Figure 5
shows the distribution of $\Delta\theta_{s}$. As can be seen, a majority of
the photons are in the first bin, corresponding to the transmitted component.
Additionally, the distribution shows a sharp cutoff at $90^{\circ}$, due to
the fact that we are only considering the first time a photon crosses the
watched volume. The middle panel of Figure 5 shows the 2-dimensional
distribution of $\Delta\theta_{\mathrm{s}}$ and initial energies, for measured
energies between $0.750-1.155$ MeV. Quite interestingly, this plot reveals a
characteristic property of the scattering; that is, as we observe at higher
angular distances from the source location, a majority of the measured photons
in a given energy band originate at increasingly higher energies. Finally, the
bottom panel of Figure 5 shows the correction factor ratio calculated for
different max values of $\Delta\theta_{\mathrm{s}}$, where we consider
1$\times$, 3$\times$, and 5$\times$ the angular resolution of the COSI balloon
instrument at 0.511 MeV ($5.9^{\circ}$, $17.7^{\circ}$, and $29.5^{\circ}$,
respectively). The calculations are made using a power law spectral model with
a photon index of 2.0. As can be seen, the ability to select photons within a
limited angular distance of the source position substantially reduces the
effects of scattering. In terms of the correction factor ratio, the effect is
$\lesssim 30\%$ when including up to 5$\times$ the angular resolution.

As a check on our calculations, we have verified that for the second
representation of the atmospheric response used here, we obtain consistent
results for the detection fraction (i.e., Figure 12) in the limiting case of
including all angles, for both transmitted and scattered components. Indeed,
the correction factor ratio for all angles in Figure 5 is consistent with the
results obtained with the other representation of the atmospheric response.

###  3.6 Validation of the Simulations

As an additional check on our simulations, we compare our results to estimates
of the atmospheric response from Takada et al. (2011), which was used to
successfully fit the growth curve of the SMILE balloon flight. To do this we
calculate the scattering ratio, defined as

$\lambda(a,Z)=\frac{F_{s}(a,Z)}{F_{t}(a,Z)+F_{s}(a,Z)},$ (16)

where $F_{t}(a,Z)$ and $F_{s}(a,Z)$ are the flux of the transmitted and
scattered components, respectively, for an atmospheric depth, $a$, and zenith
angle, $Z$. We use an atmospheric depth of $7.8\ \mathrm{g\ cm^{-2}}$,
corresponding to an altitude of 33.5 km and temperature of 234.35 K. We assume
a power law spectral model with an index of 2, and consider zenith angles
between $0^{\circ}-20^{\circ}$. These are essentially the same selections that
were used in Takada et al. (2011), with the exception that their atmospheric
depth was actually $8.0\ \mathrm{g\ cm^{-2}}$. Figure 6 shows the comparison
of the scattering fraction. The figure also includes a number of other
estimates from the literature (see Takada et al., 2011, and references
therein). Note that all literature results are from simulations or
calculations. Overall, our results are in very good agreement with the results
from Takada et al. (2011). We note that near 10 MeV our results are slightly
lower, but this is also the upper energy bound of our analysis, and thus the
comparison here is not very reliable.

Figure 6: Scattering fraction as a function of measured photon energy. Results
from this work are shown in black, which are compared to a number of other
calculations/simulations from the literature, as specified in the legend.

As another sanity check on our simulations we have analyzed the atmospheric
response using a rectangular mass model with a narrow beam as the $\gamma$-ray
source. In fact, although unknown to us when originally developing the
simulations, a similar approach was actually taken in Pendleton et al. (1995)
for balloon-borne measurements of the supernova 1987A in the hard X-ray
continuum. One of the main benefits of the rectangular mass model is that it
can provide a somewhat simpler and intuitive approach to studying the
atmospheric response, and it is less computationally intensive. Overall, we
find that the results from the rectangular mass model are consistent with
those from the spherical mass model. More details for this are provided in
Appendix D.

We have also successfully applied our atmospheric corrections to data from the
2016 COSI balloon flight. Detection of the Galactic diffuse continuum emission
during the flight was recently reported in Karwin et al. (2023). In that case,
application of the correction factor ratio was found to bring the measured
flux in closer agreement with previous measurements, as would be expected.

In summary, the simulations presented in this work have been validated in the
following ways:

  1. 1.

We have verified that the photon distributions from the simulations are in
accordance with expectations from the simulation setup.

  2. 2.

The TP calculated from simulations is in excellent agreement with the
analytical calculation.

  3. 3.

The scattering ratio is in very good agreement with other calculations from
the literature.

  4. 4.

The results from the spherical mass model were shown to be in good agreement
with results from the simplified rectangular mass model.

  5. 5.

The atmospheric corrections were successfully applied to real observational
data in Karwin et al. (2023).

##  4 Summary and Conclusion

In this work we have simulated the full response for $\gamma$-ray transport in
the atmosphere. The simulations are run with the COSI atmosphere pipeline,
which employs MEGAlib, based on Geant4. The atmosphere is characterized using
the latest version of NRLMSIS. We characterized the transport of $\gamma$ rays
in the atmosphere in terms of two components. The first component is due to
photons that scatter and never reach the detector, thereby causing an
attenuation of the original signal. The second component is from photons that
reach the detector after scattering one or more times.

Using the energy dispersion matrices from our simulations, we have calculated
the detection fractions for both photon components. Additionally, we have
defined the correction factor and correction factor ratios, which approximate
the effects of the scattered component on the measured spectrum, i.e., the
extent to which the scattered photons cause an energy-dependent distortion.
For the transmitted component, the detection fraction is equivalent to the TP,
and it gives the probability that a photon will reach the detector. Since the
photons that reach the detector do not undergo any scattering, they arrive at
the detector with the same energy and direction as they started with. The
detection fraction for the scattered component is analogous to that of the
transmitted component. The main difference is that the photons arrive at the
detector with different energies and directions compared to their starting
values.

Accounting for the scattered component is most important for diffuse sources
because photons enter the detector from all directions. We find that the
detection fraction for the scattered photons is highest for initial energies
between $\sim 0.2-0.4$ MeV. In general, the contribution from the scattered
component depends on the photon index of the source. Harder sources have more
photons at higher energies, which lose energy as they are scattered. Thus, the
scattered component is more dominant for harder sources. The end result is an
energy-dependent spectral distortion, which is highest towards lower energies.
At 0.1 MeV the scattered component may increase the flux (with respect to only
accounting for attenuation) by as much as a factor of $\sim 2-4$, depending on
the photon index and off-axis angle of the source.

For point sources observed with imaging telescopes, such as COSI, the effect
from scattering is not as important because they have the ability to decipher
the direction of the incident photons, mainly limited by the angular
resolution of the instrument. When including photons out to an angular
distance of $29.5^{\circ}$ (5$\times$ the angular resolution of the COSI
balloon instrument at 0.511 MeV), the effect from scattering is $\lesssim
30\%$, in terms of the correction factor ratio.

These results highlight the importance of accounting for photons that scatter
into the detector when making balloon-borne observations, especially for
diffuse sources. The simulation and analysis pipeline described in this work
is publicly available and is readily applicable to observations of MeV
$\gamma$ rays in the atmosphere.

## Acknowledgements

The COSI balloon program was supported through NASA APRA grants NNX14AC81G and
80NSSC19K1389. We also acknowledge support for this work under NASA APRA grant
80NSSC21K1815. This work is partially supported under NASA contract
80GSFC21C0059, and it is also supported in part by the Centre National
d’Etudes Spatiales (CNES). CMK’s research was supported by an appointment to
the NASA Postdoctoral Program at NASA Goddard Space Flight Center,
administered by Oak Ridge Associated Universities under contract with NASA.
CMK is pleased to acknowledge conversations with Alex Moiseev which helped to
improve the quality of the analysis.

##  Appendix A Photon Distributions

In Figure 7 we show the photon distributions of position and direction from
the simulation, where the first row is for the initial photons, and the second
row is for the measured photons. Column 1 shows a 3-dimensional scatter plot
of the photon position. The photons are distributed over a spherical region,
as expected. The second column shows the radial distribution. For the initial
photons we show the count rate, with respect to the radius of the surrounding
sphere disk. The rate is constant with a mean value of $0.074\mathrm{\ ph\
km^{-2}}$, and goes to zero exactly at the radius of the surrounding sphere
($r_{\mathrm{sphere}}=R_{\earth}+200$ km). This is exactly as expected. For
the measured photons the radius (with respect to Earth’s center) is exactly at
the location of the watched volume (6411.6 km). Columns 3 and 4 show all-sky
HEALPix maps of position and direction, respectively. We use an NSIDE of 16,
corresponding to an approximate angular resolution of $3.7^{\circ}$. These
distributions show that the position and direction are uniformly distributed
over the sky. Moreover, with this angular resolution we would expect $\sim
3300$ (initial) photons per pixel, consistent with what is shown.

Figure 7: The top row shows initial photons and the bottom row shows measured
photons. The columns are as follows: Column 1: 3-dimensional scatter plot of
photon position, where the color scale shows the latitude. Column 2 top:
Radial distribution for the initial photons, where we show the rate (photons
per area), with respect to the radius of the surrounding sphere disk. Column 2
bottom: Radial distribution for the measured photons, showing the counts,
where the radius is defined relative to Earth’s center (i.e., the origin of
the coordinate system). The dotted orange line corresponds to Earth’s radius,
and the dotted green line corresponds to the altitude of the surrounding
sphere. Column 3: HEALPix map showing the position on the sky. Column 4:
HEALPix map showing the direction on the sky (i.e., where the direction vector
is pointing).

The distribution of energies for both initial and measured photons is shown in
Figure 8. For initial photons, the distribution is flat across all energies,
in accordance with the simulated spectrum. The distribution of measured
photons shows more photons at lower energies, which is a consequence of the
energy loss from scattering. Note that in total there are $10^{7}$ initial
photons and $8.5\times 10^{6}$ measured photons.

Figure 8: Distributions of energy for both initial and measured photons.

Figure 9: Left: Distributions of incident angles for both measured (un-
weighted and weighted by the geometric correction factor) and initial photons,
integrated over all energies. The dashed grey curve shows the analytical
solution for the flux of a constant field through a hemisphere, with input
parameters corresponding to our simulation setup. Right: Distribution of
measured incident angles, for an initial angle of $50^{\circ}$ and energy of 8
MeV. For comparison, we show both the weighted and un-weighted histograms.

The left panel of Figure 9 shows the distributions of incident angles for both
measured (un-weighted and weighted by the geometric correction factor) and
initial photons. The distributions have a characteristic shape, peaking at
$45^{\circ}$ and falling off at lower and higher angles, including a sharp
cutoff at $90^{\circ}$. These characteristics can be well understood by
considering the flux of a constant field through a hemisphere. Specifically,
we consider the upper hemisphere from the setup shown in the left panel of
Figure 1, with the positive $\hat{z}$ direction pointing towards the top of
the page. The flux is calculated as

$\Phi=\iint_{S}\vec{F}\cdot\hat{n}\,dA,$ (A1)

where $\vec{F}$ is the count rate from the surrounding sphere disk given by
$\vec{F}=-F_{0}\hat{z}$, with $F_{0}=0.074\ \mathrm{ph\ km^{-2}}$. In
spherical coordinates we have $\hat{n}=\hat{r}$ and
$dA=R^{2}\,\mathrm{sin}\theta\,d\theta\,d\phi$, where $R$ is the radius of our
watched volume, and $\theta,\phi$ are co-latitude and longitude, respectively.
We can represent $\vec{F}$ in spherical coordinates as

$\vec{F}=-F_{0}\,(\mathrm{cos}\theta\,\hat{r}-\mathrm{sin}\theta\,\hat{\theta}).$
(A2)

Plugging everything into Eq. A1, and integrating over the azimuth angle, the
magnitude of the flux through the hemisphere can now be written as

$\Phi=2\pi R^{2}F_{0}\int\mathrm{sin}\theta\,\mathrm{cos}\theta\,d\theta,$
(A3)

from which we obtain

$\frac{d\Phi}{d\theta}(\theta)=2\pi
R^{2}F_{0}\,\mathrm{sin}\theta\,\mathrm{cos}\theta\,d\theta.$ (A4)

This equation is plotted with the dotted grey line in the left panel of Figure
9, and as can be seen, there is excellent agreement with the simulations.

In order to more clearly exemplify the effect of the weights, in the right
panel of Figure 9 we show the distribution of measured incident angles, for an
initial incident angle of $50^{\circ}$ and energy of 8 MeV, for both the
weighted and un-weighted histograms. As can be seen, the weights are important
for high off-axis angles, due to the effect of the projected area becoming
increasingly smaller.

##  Appendix B Analytical Calculation of Transmission Probability

For a narrow beam of mono-energetic photons, the change in $\gamma$-ray beam
intensity ($I$) at some distance ($x$) in a material can be expressed as:

$dI(x)=-I(x)\sigma ndx,$ (B1)

where $\sigma$ is the interaction cross section (units of $\mathrm{cm}^{2}$),
and $n$ is the number density of the material (units of
$\\#/\mathrm{cm}^{3}$). Integrating both sides along the line of sight gives:

$\int_{I_{0}}^{I_{f}}\frac{dI}{I}=\int_{x_{0}}^{x_{f}}-\sigma n\mathrm{d}x\
\implies\frac{I_{f}}{I_{0}}=\mathrm{exp}\bigg{(}\int_{x_{0}}^{x_{f}}-\sigma
n\mathrm{d}x\bigg{)}.$ (B2)

We can replace the quantity $\sigma n$ with $\eta(E)\rho(x)$, where $\eta(E)$
is the mass attenuation coefficient as a function of energy (units of
$\mathrm{cm^{2}/g}$), and $\rho(x)$ is the mass density (units of
$\mathrm{g/cm^{3}}$). Thus, the TP can be calculated as

$\mathrm{TP}\equiv\frac{I_{f}}{I_{0}}=\mathrm{exp}\bigg{(}-\eta(E)\int_{x_{0}}^{x_{f}}\rho(x)dx\bigg{)}.$
(B3)

This calculation is implemented in the cosi-atmosphere package. The mass
density is given by our atmospheric model. Note that the mass density is a
function of radius, but this can easily be mapped to any point, $x$, along the
line-of-sight, as given in Eq. B3 (e.g., using the law of cosines). Data for
the mass attenuation coefficients are taken from the National Institute of
Standards and Technology (NIST)
XCOM999https://physics.nist.gov/PhysRefData/Xcom/html/xcom1.html program. We
use a simplified atmospheric mixture consisting of 78% $\mathrm{N}_{2}$ and
22% $\mathrm{O}_{2}$. This serves as a close approximation, particularly for
altitudes $\lesssim$100 km. Figure 10 shows the mass attenuation coefficients
as a function of energy, for different relevant interactions. As can be seen,
incoherent scattering (i.e., Compton scattering) is dominant in our energy
range of interest. Pair production dominates above $\sim 30$ MeV, and photo
absorption dominates below $\sim 50$ keV. Coherent scattering (i.e., Rayleigh
scattering) also becomes important at low energies.

Figure 10: Mass attenuation coefficients as a function of energy, based on an
atmospheric composition consisting of 78% $\mathrm{N}_{2}$ and 22%
$\mathrm{O}_{2}$. The different curves show different interactions, as
specified in the legend. Pair production interactions are shown separately for
interactions with a nuclear field (NF) and an electron field (EF). The total
component is shown including coherent scattering (wCS).

##  Appendix C Spherical Mass Model: On-Axis Source

Figure 11 shows the energy dispersion matrices for an on-axis source.
Qualitatively, they are consistent with the results for the $50^{\circ}$ off-
axis source presented in Section 3.3. The corresponding detection fraction is
shown in the left panel of Figure 12. Again, they are qualitatively consistent
with the $50^{\circ}$ off-axis case. The main difference is that the detection
fraction is overall higher for the on-axis source, as expected. We note that
some statistical variation is evident in the correction factor curves. Indeed,
the on-axis case has the lowest statistics, as shown in Figure 9. Further
simulations should resolve this and produce smooth curves. In the middle and
right panels of Figure 12 we show the corresponding correction factor and
correction factor ratio, respectively.

Figure 11: Energy dispersion matrix for all photons (left), transmitted
photons (middle), and scattered photons (right), for an on-axis source. The
y-axis is the measured photon energy, and the x-axis is the initial photon
energy. The matrices are normalized by the total photons simulated in each bin
of $E_{i}$ and $\theta_{i}$. The total energy dispersion matrix is the sum of
the transmitted and scattered components.

Figure 12: Left: Transmission probability for the three components, as
specified in the legend. For reference, the TP calculated analytically is also
overlaid. Middle: Correction factor for a range of power law spectral models
ranging form 1 $-$ 3\. The solid line is for the transmitted component only,
and the dashed line is for the total component, which includes transmitted and
scattered photons. The black dash-dot line is the TP, obtained by projecting
the energy dispersion matrix onto the initial energy axis. Right: Correction
factor ratio for the same range of models shown in the left plot.

##  Appendix D Rectangular Mass Model

As an alternative approach to the spherical mass model of the atmosphere, a
rectangular geometry can be employed. In this setup, instead of using
spherical atmospheric shells, the atmosphere is described using planar slabs.
Additionally, a narrow beam is used instead of an isotropic source. Otherwise,
the atmosphere is still modeled using NRLMSIS, and the response is defined and
analyzed in a similar way.

Figure 13: Schematic depicting the setup of the rectangular mass model, as
described in the text.

Figure 14: Left: 2-dimensional scatter plot of measured photon positions, for
a $50^{\circ}$ off-axis source, using a rectangular mass model. The black
circles show the measured positions of the transmitted photons, and the blue
circles show the measured positions of all photons (transmitted + scattered).
Middle: Radial distribution of the measured photons. The left y-axis shows
counts/area (corresponding to the blue curves), and the right y-axes shows the
total counts (corresponding to the orange curve). Right: Angular distribution
of the measured photons. For reference, we show both weighted and un-weighted
histograms.

The logic behind the rectangular simulations can be understood by considering
the schematic in Figure 13. A typical (non-beamed) astrophysical source
radiates isotropically. The curvature of the wave increases with the square of
the distance from the source, and thus it can be approximated as a plane wave
once it eventually passes Earth. The plane wave can be described as the
superposition of many narrow beams. Three such beams are depicted in Figure
13. If we consider the middle beam (narrow beam 2), some photons from the beam
will reach the detector without scattering, and other photons will scatter and
never reach the detecting area. This corresponds to the transmitted component,
as we have defined in this work. If we consider the two adjacent beams (narrow
beams 1 and 3), the same thing also occurs. Thus, the same distribution of
photons that scatter out of the detector area also scatter back into it, from
the superposition of all nearby beams. This is simply a consequence of the
circular geometric symmetry of the problem. The scattered photons that enter
the detector correspond to the scattered component, as we have also defined in
this work. By watching the entire plane at 33.5 km, we can characterize the
scattered photons using a single narrow beam source (radius of 1 cm) placed at
200 km.

We simulate a $50^{\circ}$ off-axis source using the rectangular mass model.
Photons are classified as having scattered if their measured incident angle
varies by more than $0.2^{\circ}$ from the initial incident angle. Note that
the angular resolution here is much better than that used for binning the
spherical mass model. The left panel of Figure 14 shows the positions of the
measured photons, where we show separately the transmitted and total
(transmitted + scattered) components. Correspondingly, the middle panel of
Figure 14 shows the radial distribution of the measured photons, and the right
panel shows the distribution of measured incident angles. The radial
distribution shows a bimodal feature. First, there is a flat distribution (in
counts/area) between $\sim 0-1$ cm, corresponding to the beam radius of the
simulated photons. This part of the distribution corresponds to the
transmitted component. The second part of the distribution corresponds to the
scattered component, and the intensity gradually falls off with increasing
radius. In terms of total counts, the scattered component peaks near $\sim 10$
km, before quickly falling off. This indicates that a majority of the photons
the contribute to the scattered component originate from within $\sim 10$ km
of the detector. Because the radius of Earth is so large, the curvature of its
surface within 10 km is minimal101010Moving along a flat surface compared to a
spherical surface within a 10 km radius leads to an altitude difference
($\Delta h$) of $\Delta h/R_{E}=10^{-6}$., and thus we can approximate the
surface as being locally flat in the region where a majority of the scattered
photons originate from. For this reason, we expect the rectangular mass model
to give comparable results as those obtained with the spherical mass model,
and indeed, this is what we find.

Figure 15: Energy dispersion matrix for all photons (left), transmitted
photons (middle), and scattered photons (left), for a 50∘ off-axis source,
using a rectangular mass model. The y-axis is the measured photon energy (Em),
and the x-axis is the initial photon energy. The matrices are normalized by
the total photons simulated in each energy bin. The total energy dispersion
matrix is the sum of the transmitted and scattered matrices.

Figure 15 shows the energy dispersion matrices. They are very similar to the
results with the spherical mass model. Note that with this simulation setup it
is easier to obtain high statistics because we simulate a single initial
incident angle at a time, which allows for finer resolution. The detection
fraction is shown in the left panel of Figure 16. The middle and right panels
of Figure 16 show the correction factor and correction factor ratio,
respectively. Overlaid in these plots are the corresponding results from the
spherical mass model. As can be seen, they are generally in very good
agreement. The total correction factors match well, although the transmitted
component is slightly higher at low energy for the spherical mass model. This
is most likely due to the coarser binning (in angle and energy) that is used
for the spherical mass model. Likewise, the correction factor ratio for the
rectangular mass model is higher towards low energy, which can be attributed
to the difference in the binning.

As a further example, Figure 17 shows the correction factors and correction
factor ratios for an on-axis source. Again, there is very good agreement
between the rectangular and spherical mass models. This overall agreement
serves as another validation of the spherical geometry simulations, and at the
same time, it also shows that the rectangular geometry can be used as an
alternative/complementary approach to approximating the atmospheric response.
For example, this would be particularly useful to study the behaviour of
resolved $\gamma$-ray lines, for which more statistics in much smaller energy
bins is required.

Figure 16: Left: Detection fraction for the three components, as specified in
the legend, for a 50∘ off-axis source, using a rectangular mass model. For
reference, the TP calculated analytically is also overlaid. Middle: Correction
factor for power law spectral models ranging form $1-3$, for a 50∘ off-axis
source, using a rectangular mass model. The solid line is for the transmitted
component only, and the dashed line is for the total component, which includes
transmitted and scattered photons. The dotted curves show the results from the
spherical mass model. Right: Correction factor ratio for the same range of
models shown in the left plot.

Figure 17: Left: Detection fraction for the three components, as specified in
the legend, for an on-axis source, using a rectangular mass model. For
reference, the TP calculated analytically is also overlaid. Middle: Correction
factor for power law spectral models ranging form $1-3$, for an on-axis
source, using a rectangular mass model. The solid line is for the transmitted
component only, and the dashed line is for the total component, which includes
transmitted and scattered photons. The dotted curves show the results from the
spherical mass model. Right: Correction factor ratio for the same range of
models shown in the left plot.

## References

  * Agostinelli et al. (2003) Agostinelli, S., Allison, J., Amako, K., et al. 2003, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 506, 250, doi: https://doi.org/10.1016/S0168-9002(03)01368-8
  * Beechert et al. (2022) Beechert, J., Lazar, H., Boggs, S. E., et al. 2022, Nuclear Instruments and Methods in Physics Research A, 1031, 166510, doi: 10.1016/j.nima.2022.166510
  * Caputo et al. (2022) Caputo, R., et al. 2022, J. Astron. Telesc. Instrum. Syst., 8, 044003, doi: 10.1117/1.JATIS.8.4.044003
  * Chen & Buckley (2023) Chen, W., & Buckley, J. H. 2023, PoS, ICRC2023, 841, doi: 10.22323/1.444.0841
  * Emmert et al. (2021) Emmert, J. T., Drob, D. P., Picone, J. M., et al. 2021, Earth and Space Science, 8, e01321, doi: 10.1029/2020EA001321
  * Emmert et al. (2022) Emmert, J. T., Jones, M., Siskind, D. E., et al. 2022, Journal of Geophysical Research (Space Physics), 127, e2022JA030896, doi: 10.1029/2022JA030896
  * Karwin et al. (2023) Karwin, C. M., et al. 2023, ApJ, 959, 90, doi: 10.3847/1538-4357/ad04df
  * Kierans et al. (2016) Kierans, C. A., et al. 2016, PoS, INTEGRAL2016, 075, doi: 10.22323/1.285.0075
  * Orlando et al. (2022) Orlando, E., Bottacini, E., Moiseev, A. A., et al. 2022, J. Cosmology Astropart. Phys, 2022, 036, doi: 10.1088/1475-7516/2022/07/036
  * Pendleton et al. (1995) Pendleton, G. N., Paciesas, W. S., Fishman, G. J., Meegan, C. A., & Wilson, R. B. 1995, ApJ, 439, 963, doi: 10.1086/175233
  * Picone et al. (2002) Picone, J. M., Hedin, A. E., Drob, D. P., & Aikin, A. C. 2002, Journal of Geophysical Research (Space Physics), 107, 1468, doi: 10.1029/2002JA009430
  * Takada et al. (2011) Takada, A., Kubo, H., Nishimura, H., et al. 2011, ApJ, 733, 13, doi: 10.1088/0004-637X/733/1/13
  * Takada et al. (2022) Takada, A., Takemura, T., Yoshikawa, K., et al. 2022, ApJ, 930, 6, doi: 10.3847/1538-4357/ac6103
  * Tomsick et al. (2019) Tomsick, J., Zoglauer, A., Sleator, C., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 98.  https://arxiv.org/abs/1908.04334
  * Tomsick et al. (2023) Tomsick, J. A., et al. 2023, PoS, ICRC2023, 745, doi: 10.22323/1.444.0745
  * Valverde et al. (2023) Valverde, J., Kirschner, N., Metzler, Z., et al. 2023, arXiv e-prints, arXiv:2308.12464, doi: 10.48550/arXiv.2308.12464
  * Zoglauer et al. (2006) Zoglauer, A., Andritschke, R., & Schopper, F. 2006, NAR, 50, 629, doi: 10.1016/j.newar.2006.06.049