\section{Introduction}
\label{sec:intro}

The Standard Model (SM) of particle physics has proven to be remarkably successful. However, there is an abundance of empirical evidence that it is incomplete. While there are a variety of channels we can explore to discover new physics, colliders provide a clean and controlled experimental environment to identify the particle content of beyond the SM phenomena.

As we look to advance the energy and intensity frontier of collider physics, a possibility with growing interest is the construction of a TeV-scale $\mu^+ \mu^-$ collider \cite{Antonelli:2015nla,Long:2020wfp,MICE:2019jkl,Delahaye:2019omf,Neuffer:1983jr,Delahaye:2013jla}. Such a muon collider ($\mu C$) is a particularly compelling option as it affords a complementary physics program to that of a high-energy hadron collider like the LHC.
For example, with a $\mu C$ we gain access to direct couplings of both electroweak-mediated and second-generation processes~\cite{Chen:2016wkt, Buttazzo:2018qqp, Chiesa:2020awd, Costantini:2020stv, Han:2020uid, Han:2020pif, Han:2020uak, Han:2021kes, AlAli:2021let, Asadi:2021gah, Ruiz:2021tdt}.
Additionally, with increased available center-of-mass energy, we can expand our discovery prospects for massive new physics.

Since the cost of a $\mu C$---or any future high-energy collider---is substantial, it is prudent to consider possible auxiliary experiments that extend the physics program of the collider facility. An economical extension with remarkable and complementary discovery potential is a beam dump. In a beam-dump experiment, the high energy muon beam is `dumped' into a dense material to greatly increase the total rate of interaction at the price of center-of-mass energy. This experimental setup can therefore test couplings too small to be probed at the main collider by several orders of magnitude in a slightly lower mass range \cite{Bjorken:2009mm,Essig:2010gu, Kanemura:2015cxa, Chen:2017awl, Sakaki:2020mqb,Asai:2021xtg,Sieber:2021fue}.

In this letter we propose the construction of a beam-dump experiment to be included in the design of a future $\mu C$. We consider benchmark models of moderate-mass, weakly coupled new vector particles that are inaccessible at any other terrestrial experiment. First, we consider the dark photon scenario~\cite{Holdom:1985ag, Pospelov:2007mp, Arkani-Hamed:2008hhe, Pospelov:2008zw}, for which similar past proposals have focused on electron beams; see, e.g.,~\cite{Reece:2009un, Bjorken:2009mm, Essig:2010xa}. Here, the main novelty of a proposed $\mu C$ is high energy, which can provide access to  a different range of masses and (due to the large boost) lifetimes than lower-energy electron beam dumps. We also consider a model for which a muon beam is uniquely well-suited, namely the gauged flavor symmetry $L_\mu-L_\tau$~\cite{Foot:1990mn, He:1990pn,  Ma:2001md,Huang:2021nkl}. In this case, the gauge boson is produced much more copiously  from a muon beam than an electron or proton beam.
We present the projected reach of such an experiment for several generic experimental configurations.

\section{Production from a High Energy Muon Beam}
\label{sec:prodMuBeam}
\subsection{Preliminaries}

In what follows we restrict our attention to the new physics scenario of a new vector particle, which we generically call $Z'$, coupling to a current including muons. The effective Lagrangian of interest for a new $U(1)'$ gauge boson is
\begin{multline}
\mathcal{L} \supset \mathcal{L}_\text{SM}  -\frac{1}{4}Z'_{\mu\nu}Z'^{\mu\nu} + \frac{1}{2} m_{Z'}^2 Z'^\mu Z'_\mu
- \\ \sum_{l\in e, \mu, \tau}\left( i g Q_l \bar{l} \gamma^\mu Z'_\mu l  + i g Q'_l \nu_l^\dagger \bar{\sigma}^\mu Z'_\mu \nu_l \right).
\label{eq:fullLag}
\end{multline}
We consider two $U(1)$ models: dark photons $(Q_l =1, Q_l' =0, g = \epsilon e)$ and the gauged flavor symmetry $L_\mu-L_\tau$ $(Q_{\mu, \tau} = Q_{\mu, \tau}' = \pm1)$. In the latter case, kinetic mixing will also be present (at least through loops of muons and taus; see, e.g.,~\cite{Bauer:2018onh}), but the effect is both small and dependent on details of the UV completion, so we will neglect it in our discussion.

The photon--dark photon interaction is often defined via kinetic mixing in a basis of quantized charges, but the above formulation can be derived through the field redefinition $A_\mu \rightarrow A_\mu + \epsilon Z'_\mu$ to generate a coupling to the EM current. Both the dark photon and gauged flavor symmetry models are of interest at a high-energy muon collider, although they are differently motivated. Dark photons can be produced directly at any charged particle collider, but the extended reach of our proposed experiment comes from the increased center-of-mass energy.
The $L_\mu-L_\tau$ model, on the other hand, would benefit uniquely from a muon collider, as this would be the first experiment capable of direct production at high energies.
\begin{figure}[h!]
\centering
\includegraphics[width=0.4\textwidth]{figures/WWdiag}
\includegraphics[width=0.4\textwidth]{figures/WW22}
\caption{Top: the dominant bremsstrahlung production for a vector particle $Z'$ at a muon beam dump. Bottom: the same production process in the Weizs\"acker-Williams approximation.}
\label{fig:WW}
\end{figure}
\subsection{Cross Section}
The dominant production mechanism is the $2\rightarrow3$ bremsstrahlung process shown in the top of \fig{fig:WW}, where the incoming high-energy muon exchanges a virtual photon with a nucleon in the target and radiates a $Z'$~\cite{Kim:1973he, Tsai:1973py, Tsai:1986tx, Bjorken:2009mm}.
To compute this cross section, we use the Weizs\"acker-Williams approximation~\cite{vonWeizsacker:1934nji, Williams:1934ad}. For relativistic incoming muons, the exchanged photon is nearly on-shell. We can therefore approximate the full scattering process $\left(\mu (p) + N(P_i) \rightarrow \mu (p') + N'(P_f)+ Z'(k)\right)$ with the $2 \rightarrow2$ process $\left(\mu (p) + \gamma(q) \rightarrow \mu (p') + Z'(k)\right)$ shown in the bottom of \fig{fig:WW}, evaluated at minimum virtuality $t_\text{min} \equiv -q^2_\text{min}$, weighted by the effective photon flux. The cross section in the lab frame is
\begin{multline}
\frac{d\sigma(p + P_i \rightarrow p' + k + P_f)}{dE_{Z'}d \cos \theta_{Z'}} = \left( \frac{\alpha_\textsc{EM} \chi}{\pi}\right) \left( \frac{x E_0 \beta_{Z'}}{(1-x)} \right)\\ \times \frac{d \sigma (p+q \rightarrow p' + k)}{d(p \cdot k)}\bigg \rvert_{t=t_\text{min}},
\label{eq:fullWW}
\end{multline}
where $E_0$ is the energy of the incoming muon beam, $x\equiv E_{Z'}/E_0$ is the fraction of energy of the $Z'$, $\theta_{Z'}$ is the angle of emission, and $\beta_{Z'} \equiv \sqrt{1-m_{Z'}^2/E_{Z'}^2}$. The effective photon flux is parameterized by $\chi$, defined as
\begin{equation}
\chi \equiv \int_{t_\text{min}}^{t_\text{max}} dt\, \frac{t-t_\text{min}}{t} G_2(t),
\end{equation}
where $G_2(t)$ is the electric form factor of the target atom, including both atomic and nuclear and elastic and inelastic effects, following the approximation in~\cite{Bjorken:2009mm}.

This approximation scheme is valid in the regime of highly relativistic beam particles and emitted vector particles:
\begin{equation}
\frac{m_\mu}{E_0}, \ \ \frac{m_{Z'}}{xE_0}, \ \ \theta_{Z'} \ll 1.
\end{equation}
After integrating out the angular dependence, the differential cross section in $x$ is
\begin{equation}
\frac{d\sigma(2\rightarrow3)}{dx} = \frac{8\alpha_\textsc{EM}^2 \alpha_g Q_\mu^2 \chi  \beta_{Z'}}{m_{Z'}^2 \frac{1-x}{x} + m_\mu^2 x} \left( 1 - x + \frac{x^2}{3}\right),
\label{eq:cs2to3}
\end{equation}
where $\alpha_g \equiv \frac{g^2}{4\pi}$. The normalized distribution of \eq{eq:cs2to3} is plotted in \fig{fig:sigmaDist} for various vector masses $m_{Z'}$. For all values of mass we find that the probability of emission has support primarily in the highly relativistic regime.
\begin{figure}[h!]
\includegraphics[width=0.4\textwidth]{figures/sigmaDistLog}
\caption{Distributions of the fraction of the beam energy $x$ carried away in the emission of the $Z'$ for various masses. The vertical lines indicate where the boost factor $\gamma=100$ for the corresponding mass, thus for all masses the majority of emission occurs to the right of these lines in the highly relativistic regime.}
\label{fig:sigmaDist}
\end{figure}

Note that in the context of an electron beam dump experiment, the term proportional to the beam particle mass in the denominator of \eq{eq:cs2to3} can be neglected.
Here, on the other hand, we are interested in $m_{Z'}$ near the muon mass, so we cannot drop this term.
\subsection{Signature}
The signal of interest for this experimental setup is a dilepton final state ($e^+e^-$ or $\mu^+\mu^-$). From the Lagrangian defined in \eq{eq:fullLag}, the decay rate to massive leptons is
\begin{equation}
\Gamma\left( Z' \rightarrow l^+ l^-\right) = \frac{g^2 Q_l^2}{12\pi} m_{Z'}\left( 1 + \frac{2m_l^2}{m_{Z'}^2} \right)\sqrt{1 - \frac{4 m_l^2}{m_{Z'}^2}}.
\label{eq:zPrimeWid}
\end{equation}
We restrict our attention to this signature for ease of detection prospects.

Before we proceed with the details, we can check that a multi-TeV muon collider provides the right environment to extend the search boundaries of these models. Taking the dark photon as a benchmark, the approximate decay length (estimated from \eq{eq:zPrimeWid})
\begin{multline}
l_{Z'} \equiv c \gamma \tau_{Z'} \approx x \left( \frac{E_\text{0}}{\text{TeV}} \right)  \left( \frac{\text{GeV}}{m_{Z'}} \right)^2 \left( \frac{10^{-7}}{\epsilon} \right)^2 \times 10 \text{ m},
\label{eq:cgammatau}
\end{multline}
where $x \sim 1$.
This suggests that with a modest-size experiment, a TeV beam dump can dramatically expand the reach of these vector models.

\subsection{Number of Signal Events}
The differential number of signal $Z'$ events per energy fraction $x$ and position $z$ along the beamline is given by the equation
\begin{multline}
\frac{dN}{dx\, dz} = N_\mu \frac{X_0}{m_T} \int_{E_{Z'}}^{E_0} dE_1 \int_0^T dt\, I(E_1; E_0, t)
\\ \times \left( \frac{E_0}{E_1} \frac{d\sigma}{dx'}\right)_{x' = \frac{E_{Z'}}{E_1}} \frac{dP(z-\frac{X_0}{\rho}t)}{dz},
\label{eq:mastEq}
\end{multline}
where $N_\mu$ is the number of incident muons on target, $m_T$ is the mass of a target atom, and $X_0$ and $\rho$ are the unit radiation length and density of the target respectively.\footnote{In the literature, $1/m_T$ is often written as $N_0 / A$ where $N_0$ is Avogadro's number and $A$ is the atomic number of the target. This assumes that $X_0$ is measured in particular units, g/cm$^2$.}
We parameterize the position along the length of the target in terms of the dimensionless parameter $t$ from $0$ to $T$, such that the full length of the target is $L_\text{tar} = \frac{X_0}{\rho}T$.

The energy of the beam particle after radiative losses in the material is modeled by the function $I(E_1; E_0, t)$. However, the effective radiation length of a muon in reasonable target materials (e.g., water or lead) is 50~m to 1~km~\cite{Zyla:2020zbs}. For the proposed experimental setup, we can safely assume zero radiative losses and replace
\begin{equation}
I(E_1; E_0, t) = \delta(E_1 - E_0).
\end{equation}

The differential probability of $Z'$ decay is given by
\begin{equation}
\frac{dP(z)}{dz} = \frac{1}{l_{Z'}} e^{-z/l_{Z'}}.
\end{equation}
As indicated in \eq{eq:cgammatau}, the decay length $l_{Z'}$ is a function of the $Z'$ mass $m_{Z'}$ and energy $x E_0$.
After integrating over the muon beam energy and target thickness, \eq{eq:mastEq} simplifies to
\begin{multline}
\frac{dN}{dx} = N_\mu \frac{\rho\, l_{Z'}}{m_T} \frac{d\sigma}{dx} \times \\
\left( e^{L_\text{tar}/l_{Z'}} - 1\right)e^{-(L_\text{tar}+L_\text{sh})/l_{Z'}} \left(1 - e^{-L_\text{dec}/l_{Z'}} \right).
\label{eq:masterEq}
\end{multline}
where the various length scales are illustrated in \fig{fig:expSetup}.

\begin{figure}[h!]
\includegraphics[width=0.4\textwidth]{figures/beamdumpExp_corporate}
\caption{Schematic of a muon beam dump experiment. The lengths are not drawn to scale. }
\label{fig:expSetup}
\end{figure}
\section{Experimental Setup}
\label{sec:expSetup}
The experimental setup of a beam dump experiment is shown in \fig{fig:expSetup}. The high energy muon beam is dumped into a solid material target of length $L_\text{tar}$. Immediately after the target is a shield of length $L_\textrm{sh}$ dedicated to removing the residual beam and any background. This includes both a region with a strong magnetic field of length $l_m$ and a veto mechanism to identify or remove any remaining particles. Beyond the shielding is the fiducial decay region of length $L_\text{dec}$ before the detector. To be detected, signal events must be produced in the target and decay in fiducial region.

From \eq{eq:masterEq} it is clear that the length scales of the various components of the experiment have dramatic impact on the number of signal events observed. New particles that are produced too early get shielded or vetoed and are missed, but particles that live too long escape the experimental hall altogether.

For concreteness, we choose experimental parameters such as length scales at reasonable orders of magnitude. The purpose of this letter is to demonstrate the feasibility of beam dump experiments at future muon colliders. Since we cannot know the available technology or experimental design constraints at this point, we do not attempt to optimize the experimental design.
\subsection{Target}
In this work we consider a target made of water at standard temperature and pressure, with a length of $L_\text{tar} = 10$ m. Other economical and common choices for a beam dump target include pressurized water, lead or tungsten. Higher density materials will increase the cross section as the effective photon flux increases, but may be too expensive or infeasible to install at the beam dump site.
\subsection{Shielding}
Since the $Z'$ and its subsequent decay products are highly boosted, it is imperative to remove the central muon beam for a near zero-background signal. To do so, the shielding region must remove the muon beam from the geometric acceptance of the detector. Since high-energy muons are extremely difficult to absorb, we instead propose to deflect the beam using a strong magnetic field.

The size of the detector is set by the opening angle of the dilepton decay products of the boosted $Z^{\prime}$. The decay products are produced extremely forward with an emission angle in the lab frame of $\theta_\textrm{max} \lesssim m_{Z'}/xE_0$. Note that, in the parameter space of  interest where $m_{Z'} > m_\mu$, the maximum angle of $Z'$ emission $\theta^{Z'}_\textrm{max} \sim \frac{\sqrt{m_{Z'} m_\mu}}{E_0}$ is parametrically smaller and therefore negligible.

The shielding magnet must be sufficiently strong to divert a roughly TeV-scale muon beam at a greater angle $\theta_\text{mag} > \theta_\textrm{max}$. If we assume the magnetic field is constant over a length $l_m$ just after the target, the field strength $B$ must be approximately
\begin{equation}
\frac{B}{\mathrm{Tesla}} \sim \theta_\text{max} \times \frac{E_0}{\text{GeV}} \times \frac{\text{meter}}{l_m}.
\end{equation}
For the parameter space of interest, this corresponds to sub-Tesla magnetic fields over a few meters, which is comparable to or smaller than similar vetos proposed in other future beam-on-target experiments \cite{SHiP:2015vad}.

Since bending a high-energy beam of charged particles will emit some radiation, there must also be a passive veto beyond the magnetic field. This veto can be made of a combination of dense material to absorb photons, such as iron or lead, as well as a charged particle tracker to catch any charged remnants.

We assume that the total length of this shielding region will be contained in $L_\text{sh} = 10$ m, but we discuss the consequences of a larger $L_\text{sh}$ later with the results.
\subsection{Detector}
Once a $Z'$ has been produced in the target and propagates beyond the shielding region, it must decay in the fiducial region before the detector in order to be observed. As previously mentioned, our signal of interest is a dilepton resonance ($e^+e^-$ or $\mu^+\mu^-$). This signal can be observed with a relatively simple detector setup: a charged particle and muon tracker. Again, as the state-of-the-art detection technology is not yet known, we do not provide a detailed description of the system. The size of the detector $d$ must be roughly
\begin{equation}
d \sim \theta_\text{max} L_\text{dec}.
\label{eq:detectorSize}
\end{equation}

From \eq{eq:masterEq}, it is clear that a longer experimental hall increases the number of $Z'$ particles that decay in acceptance. As a reasonable benchmark we set $L_\text{dec} = 100$ m, which corresponds to a detector size of $d\lesssim$ 2 m for the relevant range of masses. Note that for this exploratory study, we consider a minimally instrumented scenario. However, one could envision more sophisticated detection scenarios that involve instrumenting along the fiducial region and accounting for missing energy signals as well.

\section{Reach for New Gauge Forces}
\label{sec:reach}
We present the reach of both a dark photon and $L_\mu - L_\tau$ $Z'$ model. For concreteness we consider the reach with a 1.5 TeV beam (corresponding to a 3 TeV collider), a standard benchmark in $\mu C$ literature \cite{Delahaye:2019omf}.

For the dark photon scenario, we show the existing constraints from $e^+ e^-$ or $\mu^+\mu^-$ resonance searches at BaBar~\cite{BaBar:2014zli}, NA48~\cite{Na482:2015wmo}, the A1 Experiment at the Mainz Microtron~\cite{Merkel:2014avp}, KLOE~\cite{KLOE-2:2011hhj, KLOE-2:2012lii, KLOE-2:2014qxg, KLOE-2:2016ydq},  and LHCb~\cite{LHCb:2019vmc}; previous beam dump experiments, such as E141~\cite{Riordan:1987aw} and E137~\cite{Bjorken:1988as, Batell:2014mga, Marsicano:2018krp} at SLAC, E774 at Fermilab~\cite{Bross:1989mp}, CHARM~\cite{CHARM:1985anb, Gninenko:2012eq} and NuCal~\cite{Blumlein:1990ay, Blumlein:2011mv, Blumlein:2013cua}; as well as constraints from Supernova 1987A~\cite{Chang:2016ntp} in gray.
Additionally we plot the projected reach from other future experiments including Belle-II~\cite{Belle-II:2018jsg}, LHCb~\cite{Ilten:2015hya, Ilten:2016tkc}, SHiP~\cite{Alekhin:2015byh}, and AWAKE~\cite{Caldwell:2018atq}.

The projected sensitivity of the $\mu C$ beam dump to a dark photon $Z'$ is shown in \fig{fig:dpReach}. Since the number of muons delivered to target cannot be known at this stage, we provide three reach curves reflecting conservative to optimistic projections of $N_\mu$. A discussion of these choices can be found in the appendix.
\begin{figure}[h!]
\centering
\includegraphics[width=0.45\textwidth, left]{figures/dpReach.pdf}
\caption{Contour plots indicating 5 signal events detected with $N_\mu =10^{18}, \ 10^{20}, 10^{22}$ with a beam energy $E_0 = 1.5$ TeV. The dips in the contours near $m_{Z'} = 1$ GeV occur when there is resonant production in the $Z' \rightarrow$ hadron decay channel, thus reducing the dilepton branching ratio.}
\label{fig:dpReach}
\end{figure}
The experiment would expand the reach in parameter space not only beyond existing constraints, but also in complementary regions to other future experiments. The gain in coverage occurs mainly in the directions of larger coupling and mass. This is due to the high energy of the beam, and therefore production of highly boosted $Z'$ particles.

The boundary of the discovery region at large couplings occurs when the $Z'$ decays too early and is vetoed. However, a relativistic $Z'$ will live longer in the lab frame and therefore can decay in the fiducial region. At this unprecedented beam energy, $Z'$s at higher masses than before are sufficiently boosted to survive past the veto. If the shielding length $L_\text{sh}$ is extended, then sensitivity degrades in the large coupling regime while leaving the bottom edge unmoved. However, if the beam energy $E_0$ is increased, both the upper and lower boundaries are shifted upwards to higher couplings.
\begin{figure}[h!]
\centering
\includegraphics[width=0.45\textwidth, left]{figures/lmultaureachConst}
\caption{Same as \fig{fig:dpReach} but for the $L_\mu-L_\tau$ model. The sensitivity is bounded at lower $m_{Z'}$ by the dimuon production threshold, since the electron channel is not open in this model.}
\label{fig:lmultauReach}
\end{figure}

The parameter space of the $L_{\mu} - L_{\tau}$ model is notably unconstrained in the region $g \lesssim 10^{-3}$ and $m_{Z'} \gtrsim10$ MeV.
Existing constraints come from measurements of the primordial abundances of light nuclei~\cite{Escudero:2019gzq}, from observations of SN1987A \cite{Croon:2020lrf},
from measurements of the anomalous magnetic moment of the muon,\footnote{Given the current discrepancy between the theoretical prediction~\cite{Aoyama:2020ynm} and experimental measurement of $(g-2)_{\mu}$~\cite{Muong-2:2021ojo}, we take the $5\sigma$ upper limit as a constraint, and show the $2\sigma$ preferred region in green in \fig{fig:lmultauReach}.} limits from neutrino trident production~\cite{Altmannshofer:2014pba, CCFR:1991lpl} and searches for $e^+e^- \to \mu^+\mu^- Z'(\mu^+\mu^-)$ at BaBar~\cite{BaBar:2016sci}.
The current bounds are shown in grey in \fig{fig:lmultauReach}. We also show the projected limits from other muon beam experiments, M$^3$~\cite{Kahn:2018cqs} and NA64$\mu$~\cite{Sieber:2021fue}.
Note that other proposed experiments such as Ref.~\cite{Krnjaic:2019rsv} might have comparable reach to NA64$\mu$ and M$^3$.
The reach plot, drawn with the same values of $N_\mu$ and beam energy $E_0$, is shown in \fig{fig:lmultauReach}. The coverage is completely separated from other constraints on this model due to the novelty of both the beam of muons and the energy of the beam.

\section{Conclusions and Future Work}
\label{sec:conc}
A future multi-TeV muon beam dump experiment would provide a window into previously unexplored parameter space for a variety of motivated new physics models. The sensitivity improvements beyond other similar proposed experiments, such as Refs.~\cite{Bjorken:2009mm, Essig:2010gu, Kanemura:2015cxa, Chen:2017awl, Sakaki:2020mqb,Asai:2021xtg,Sieber:2021fue}, stem from two features of such an experimental setup: the increased beam energy and direct coupling to muons.

 As discussed, the unprecedented beam energy translates to boosted new particles with extended lifetimes. Therefore couplings that would otherwise be too large to be detected at previous beam dump experiments would be accessible. Additionally, a muon collider is uniquely well suited to study models with couplings to muons. In this letter we computed the reach of the gauged $L_{\mu} - L_{\tau}$ symmetry as a motivated example, but this broadly applies to more general dark sectors with non-universal fermion interactions.

In the proposed detection strategy, we've taken a minimalist approach to instrumentation. However, in the event of an observed resonance, additional detectors to identify the rate into taus could be used to determine the underlying theory.
If this experiment were to confirm the existence of a gauged flavor symmetry, this would be significant for several areas of particle physics. A gauged $L_\mu - L_\tau$ symmetry could explain the near-maximal mixing between muon neutrinos and tau neutrinos~\cite{Ma:2001md}. It has recently been observed that an $SU(3)$ extension of this group can give rise to a complete model of lepton masses~\cite{Alonso-Alvarez:2021ktn}.
Finally, we\footnote{One of us, anyway.} cannot resist noting that every measured gauge coupling to date is an $O(1)$ number, while \fig{fig:lmultauReach} shows that a discovery of $L_\mu - L_\tau$ at the beam dump would necessarily imply a tiny gauge coupling $\lesssim 10^{-5}$, which would become a powerful constraint on UV physics~\cite{Arkani-Hamed:2006emk}.

While we have focused on models with new vectors, the improved reach would be similarly impressive for other dark sectors with non-vector mediators. Additional searches should include models with new muon-philic scalars or pseudoscalars such as axions. These particles could be produced either directly from the muon or from photon fusion by effective interactions. A detailed study of the reach for these models will be presented in a future publication.

In the current era of particle physics, no stone can be left unturned when searching for new phenomena. A future muon collider and a corresponding beam dump experiment would greatly enhance our sensitivity in novel and complementary regimes to our current and past experimental program.

\section*{Acknowledgments}
We are grateful to Cliff Cheung, Matheus Hostert, Simon Knapen, Johannes K. L. Michel, Cristina Mondino, Clara Murgui, Simone Pagan Griso, Matthew Strassler, and Jesse Thaler for useful discussions.
CC, SH and MR are supported by the DOE Grant DESC0013607. CC is also supported by an NSF Graduate Research Fellowship Grant DGE1745303.
SH and MR are also supported in part by the Alfred P. Sloan Foundation Grant No.~G-2019-12504.
RKM is supported by the National Science Foundation under Grant No.~NSF PHY-1748958 and NSF PHY-1915071.
MR is also supported by the NASA Grant 80NSSC20K0506.

\section{Notes on the Expected Number of Muons on Target}
\label{app:luminosity}

The number of muons available at a beam dump experiment associated with a high-energy collider depends significantly on the beam design, for which at this time there are only rough proposals.
Parametrically, the number of muons in a collider beam, per unit time, is given by
\begin{equation}
{\dot N}_{\mu} = f_0 \times n_b \times N_{\mu / \textrm{bunch}}
\end{equation}
where $f_0$ is the muon source repetition rate (i.e., the rate at which the muons are generated by the source, which could be either a proton or a positron beam), $n_b$ is the number of colliding bunches in each beam, and $N_{\mu /\textrm{bunch}}$ is the number of muons per bunch.
Assuming the MAP design parameters~\cite{Delahaye:2013jla, Neuffer:2018yof} ($f_0 = 5\,\textrm{Hz}$, $N_{\mu / \textrm{bunch}} = 2 \times 10^{12}$, $n_b = 1$), this translates to ${\dot N}_{\mu} \sim 10^{20} / \textrm{year}$ generated for collisions.
The actual number of muons available at the beam dump target of course will be smaller, and will depend on how frequently the beam is dumped.
To our knowledge, this frequency is not discussed in any existing design proposals for a muon collider, and will have to be optimized to maximize the number of high-energy collisions while mitigating neutrino radiation and maintaining a high quality beam. Assuming, however, that the muons are dumped after approximately a muon lifetime, $N_{\mu} \sim 10^{20}$ in a multi-year run may be achievable.

While the MAP proposal envisions a proton driver as the source of the muon beam, an alternative possibility is to use a positron source, as in the LEMC scenario. A positron source would have a much smaller number of muons per bunch, but compensates with a much higher repetition rate, $f_0$. The current proposals obtain much smaller luminosities at a high-energy collider, and correspondingly lower number of muons on target for an associated beam dump ($N_{\mu} \lesssim 10^{18}$ assuming the LEMC parameters in ref.~\cite{Neuffer:2018yof}), but future technological advances may significantly increase the muon yield.

\bibliographystyle{utphys}


\section{Appendix}
\label{sec:appendix}

\subsection{Sequence Length Distribution} \label{app:token-seq-len}
As described in~\Secref{sec:token-sub}, we try four different tokenizers to substitute for our \texttt{RoBERTa}~\cite{liu2019roberta} model that uses the Byte-Pair Encoding (BPE) \citep{sennrich2015neural} tokenizer. Specifically, we substitue with the WordPiece tokenizer~\cite{wu2016google} used by \texttt{BERT}~\cite{devlin2018bert} (i.e., $\texttt{BERT}$ Tokenizer in \Tabref{tab:badglish-examples}) and the SentencePiece tokenizer~\cite{kudo2018sentencepiece} used by \texttt{Albert}~\cite{lan2019albert} (i.e., $\textbf{\texttt{Albert}}$ \textbf{Tokenizer} in \Tabref{tab:badglish-examples}). Additionally, we substitute with two new non-English tokenizers including the French \texttt{FlauBERT}~\cite{le2020flaubert} ($\textbf{\texttt{FlauBERT}}$ \textbf{Tokenizer} in \Tabref{tab:badglish-examples}) and the Dutch \texttt{DutchBERT}~\cite{devries2019bertje} ($\texttt{DutchBERT}$ Tokenizer in \Tabref{tab:badglish-examples}). As shown in \Figref{fig:tokenizer-seq-len}, we plot the distributions of sequence lengths as a measure of the heterogeneity introduced by new tokenizers to ensure variences across tokenized sequence lengths. Specifically, we see there are inferior tokenizers such as $\texttt{FlauBERT}$ Tokenizer with a 22.15\% increase in sequence length. Our results are consistent with previous findings~\cite{rust2020good} where sequence length distributions are closer

\subsection{Training Set-up} \label{app:training-setup}

\paragraph{Downstream Task.} We use the GLUE benchmark~\cite{wang2018glue} to evaluate model performance, which covers nine different NLP tasks. We report scores on the development sets for each task by fine-tuning our pre-trained or mid-tuned models. We fine-tune for 5 epochs for the smaller datasets (WNLI and MRPC) and 3 epochs for the others. For the performance metrics, we use Matthew's Correlation for CoLA, Pearson correlation for STS-B, and accuracy for all the other datasets.

\paragraph{Hyperparameter and Infrastructure.} For each of the mid-tuning and fine-tuning experiments, we collect averaged results from 3 runs with distinct random seeds. We tune our models with two learning rates $\{2e^{-5}, 4e^{-5}\}$, and report the best results from these two learning rates. Fine-tuning with 9 GLUE tasks takes about 8 hours on 4 NVIDIA Titan 12G GPUs. Mid-tuning with our subset of WikiText-103M corpus takes about 18 hours with the same infrastructure.

\input{figures/tokenizer-seq-len}

\subsection{GLUE Task Performance} \label{app:per-task-perf}
\Tabref{tab:per-task-perf} shows performance break-down for individual GLUE task under different transformations as described in \Secref{sec:transformation-type}.

\input{tables/per-task-results}

\section{Introduction}
What makes it hard for neural networks to learn new languages? Despite their strengths, large-scale pretrained language models (LLMs) require very large datasets for pretraining, making it challenging to train LLMs from scratch for low-resource languages~\cite{devlin2018bert, liu2019roberta, lacoste2019quantifying, clark2020electra}. For such languages, an appealing approach is to \textit{transfer} knowledge from an LLM trained for a high-resource language, since such models may transfer knowledge even under extreme conditions~\cite{papadimitriou2020learning, tamkin2020}.
A range of training strategies and model architectures have been explored to enable such crosslingual transfer of English LLMs, leveraging different techniques including adaptive pretraining~\cite{reimers2020making}, and embedding re-training~\cite{artetxe2020cross, tran2020english}. However, progress is challenging as it is hard to distinguish factors that effect such transfer.

\input{figures/main-figure-remake}

In this paper, we present a set of \textit{controlled transfer studies} for determining what makes cross-lingual transfer hard. We focus on three factors salient to crosslingual transfer: the \textbf{embedding layer}, the \textbf{tokenizer}, and \textbf{syntactic shifts}. We construct a set of systematically transformed versions of GLUE (t-GLUEs) targeting each of these factors, and observe how the finetuning performance of a pre-trained English \texttt{RoBERTa}~\cite{liu2019roberta} model degrades as a result of each of these transformations.
Crucially, our method allows us to disentangle the effects of correlated factors: while all factors would change at once if we transfer between natural languages, our transformations allow us to pinpoint what causes difficulties for transfer.

First, we focus on vocabulary effects: what is it that makes the embeddings important in cross-lingual transfer? We separately test the effects of having a consistent word embeddings layer, the
morphological consistency of the tokenizer, and the all-around quality of tokenization. Surprisingly, we find that the effect of word embedding alignment overshadows any aspect of tokenizer quality that we test. Second, we test how all of these vocabulary-level factors compare to the effects of structural syntactic perturbations, and find that structural perturbations have a significantly smaller effect on transfer ability than vocabulary effects.

\section{Related Work}

Cross-lingual transfer studies on multi-lingual models such as Multilingual BERT~\cite{devlin2018bert} demonstrate the utility of multilingual training for producing models with parallel representations that ease zero-shot transfer~\cite{pires2019multilingual}.
However, it is unlikely that a unified model can cover the more than 6,900 languages in the world \cite{feng2020language, wang2020negative}.
Enabling cross-lingual transfer from one language to another remains an important tool for expanding large-scale models to more languages, and
we aim to examine such cross-lingual transfer with well-defined experiments covering distinct linguistic factors.

Our experiments build off previous efforts that try to enable crosslingual transfer from pretrained monolingual LLMs to new languages~\cite{artetxe2018generalizing, artetxe2020cross, tran2020english, reimers2020making, gogoulou2021cross}. For example, BERT shows effective but limited transfer when attached to multi-lingual tokenizers~\cite{reimers2020making}. Recent work further improves this approach through continued multi-stage pretraining on foreign languages, which helps the English LLMs acquire new word embeddings~\cite{tran2020english}. In this work, we aim to provide quantitative insights into how these factors, such as static embeddings similarity, affect crosslingual transfer through controlled studies.

\input{tables/badglish-examples}

\section{Transformed English (t-Englishes)}
\label{sec:t-englishes}

\label{sec:transformation-type}
\Figref{fig:overview} shows our main experimental paradigm. We study cross-lingual knowledge transfer from English to t-Englishes: versions of English that are systematically altered through different types of transformations. Here we explain and motivate each t-English variant that we study, organized by the main component they consider in crosslingual transfer: the embedding layer, the tokenizer, and syntactic shifts. A full list with examples can be found in \Tabref{tab:badglish-examples}.

\subsection{Embedding Layer}

Previous works have consistently found that good embeddings are crucial for enabling effective cross-lingual transfer~\cite{tran2020english, artetxe2020cross}. However, these gains could be due to several factors, including better initialization statistics \citep{raghu2019transfusion}, or to a learned alignment between the learned embeddings and the pretraiend transformer layers \citep{wu2021identifying}. We consider a baseline setting, where we reinitalize the embedding layer such that the optimization during fine-tuning has to fully reconstruct the semantic space, and compare it against two other perturbations which preserve the statistics of the embedding weights:
1) a \textbf{token swap} scenario, where we scramble the identities of tokens so that the meaning of each token is not represented by its row in the word embeddings matrix, and 2)  a more challenging \textbf{word swap} scenario where the identities of whole words are scrambled before tokenization. In the word swap scenario, the morphological and sub-word consistency of the tokenizer is lost: if the word ``cat'' gets swapped to ``audible'', and the word ``bravery'' gets swapped to ``visible'', then after tokenization the ``-ible'' subword does not represent a consistent morpheme.

\subsection{Tokenizer}

\label{sec:token-sub}
Tokenizers have been shown to play an important role in crosslingual transfer for multi-lingual models~\citep{rust2020good}. In this section, we investigate how changes in tokenization between the source and target languages impacts downstream performance, holding all other factors constant. To do so, we substitute the Byte-Pair Encoding (BPE) \citep{sennrich2015neural} of the \texttt{RoBERTa} model with the WordPiece tokenizer~\cite{wu2016google} used by \texttt{BERT}~\cite{devlin2018bert} and the SentencePiece tokenizer~\cite{kudo2018sentencepiece} used by \texttt{Albert}~\cite{lan2019albert}. We also experiment with two non-English tokenizers, which produce lower-quality tokenizations for English:
(see~\Appref{app:token-seq-len}): the French \texttt{FlauBERT}~\cite{le2020flaubert} and the Dutch \texttt{DutchBERT}~\cite{devries2019bertje}. Examples of each tokenization scheme are presented in \Tabref{tab:badglish-examples}. For new tokenizers, word embedding weights are reinitialized accordingly.

\subsection{Syntactic Shifts}

\label{sec:syntactic-shifts}
While syntax is an crucial aspect of language, studies have also shown syntactic typology to be surprisingly non-predictive of transfer quality \citep{pham-etal-2021-order}, and other studies have shown LLMs to be largely word-order invariant \citep{sinha2021masked}.
We investigate a set of syntactic transformations that isolate syntactic word-order shifts from the other factors that can vary between languages such as tokenization, static embeddings, and morphological representation.
We use the open-source package \texttt{Galactic-Dependencies}~\cite{wang2016galactic} to transform the word order of our English training corpora to match the noun-phrase order and the verb-phrase order of French and Japanese ($\{\textbf{N}_{\text{fr}}, \textbf{V}_{\text{fr}}\}$ and $\{\textbf{N}_{\text{ja}}, \textbf{V}_{\text{ja}}\}$ in \Tabref{tab:badglish-examples}) and also perform a mixed transformation with French noun order and Japanese verb order ($\{\textbf{N}_{\text{fr}}, \textbf{V}_{\text{ja}}\}$ in \Tabref{tab:badglish-examples}).
We use the open-source package \texttt{Stanza}~\cite{qi2020stanza} for sentence segmentation and parsing
before scrambling word orders. We also test the the stronger transformation of completely reversing the word order, as well as randomly shuffling the word order to produce a lower bound on the amount of useful syntactic information.

\section{When Does Continued Pretraining Help?}

Instead of finetuning LLMs directly on target language tasks, continued pretraining on a target language corpus may help with model adaptation~\cite{artetxe2020cross, tran2020english}. Here, we examine how continued pretraining fares as a potential solution to the different shifts we consider in \Cref{sec:t-englishes}; this corresponds to the setting where there might exist a moderate amount of unlabeled data available for a low-resource language. Formally, we continued pretraining \texttt{RoBERTa} using the masked language modeling objective on a t-English corpus before evaluating its finetuning performance with downstream tasks. We use a subset of WikiText-103M corpus~\cite{merity2016pointer} containing approximately 15\% of examples for mid-tuning.\footnote{For comparison, the pretraining data for $\text{BERT}_{\text{BASE}}$ contains 3.3B tokens~\cite{devlin2018bert}. Here, we have about 15M tokens which is about 0.45\% of its pretraining data.} For each t-English variant, we ensure the number of steps for continued pretraining are kept constant.

\section{Results}
We use the GLUE benchmark~\cite{wang2018glue} to evaluate model performance, which consists of nine different NLP tasks. We produce t-GLUEs by transforming both the training and validation sets, and report scores on the transformed validation sets after finetuning our pre-trained or continued pretrained models. Each experiment is run three times with different random seeds. We include details in \Appref{app:training-setup}.

\subsection{Good Embeddings Are Most of What You Need}

As expected, we find a very large drop in GLUE performance (-40\%) from reinitializing the embeddings layer (\Cref{fig:perf-identity-swap}), indicating the importance of good embeddings to crosslingual transfer. However, the token swapping results reveal that only half of that drop (-21\%) is attributable to the token-specific information stored in each embedding, with the remainder due to merely having embeddings that lie within a plausible ``language space'' for the model. Surprisingly, scrambling the whole words rather than tokens and breaking the morphological consistency of the resulting tokenization does not have a larger effect than the one-to-one token swap, suggesting minimal additional contribution of subword information in this regime.

Interestingly, we find that continued pretraining closes the gap between the reinitialized and swapped conditions. While this eliminates over half of the performance drop of for the reinitialized t-English (-40\% to -18\%) the gain is much smaller for the swapped varaints (-21\% to -18\%). This suggests that other solutions to the word alignment problem besides continued pretraining will be needed in resource-scarce settings.

\input{figures/identity-swap}

\input{figures/tokenizer-substitution}

\input{figures/syntactic-shifts}

\subsection{Bad Embeddings Can Ruin A Good Tokenizer}\label{sec:token-results}
Surprisingly, we see virtually no effect of using different different tokenizers (\Cref{fig:tokenizer}), despite some of the tokenizers being very different from the original BPE tokenizer.
Similarly, continued pretraining also fails to reveal differences when different tokenizers are used.
These findings suggest that performance drops seen during crosslingual transfer are mostly due to the information loss of re-initializing the embedding weights, rather than the quality of tokenization.
Tokenization may thus be a ``lower-order bit'' for crosslingual transfer, which has little impact until good word embeddings are learned.

\subsection{Syntax Matters, But Not Too Much}  \label{sec:syntax-results}
As shown in \Figref{fig:perf-syn-shifts}, a \texttt{RoBERTa} model finetuned on randomly ordered GLUE sentences experiences a performance drop of 13\%. This is a markedly smaller drop than the word-shuffling experiments (-19\%), indicating that a bag-of-words classifier with good embeddings performs better than a model which experiences no syntactic shifts but that lacks good word embeddings.

This performance gap closes appreciably as we perform more structured syntactic shifts such as reversing the sentence (a drop of 10\%), or systematically permuting word orders using the dependency tree (a drop of between 7\% and 9\%). Rather than being invariant to word orders across natural language understanding tasks~\cite{sinha2021masked, pham-etal-2021-order}, we instead find that BERT-based models are in fact sensitive to word order, at least for the tasks in the GLUE benchmark. In addition, we find that continued pretraining can close the performance gap to all but a few percentage points for tree-based structural shifts. These results suggest that syntactic shifts have real but limited impact on crosslingual transfer, compared to embedding layer effects.

\section{Conclusions}
In this paper, we propose a novel paradigm to study cross-lingual transfer through transformations which simulate the linguistic changes across languages. In contrast to what prior work implies, we find significant effects from syntactic shifts and no effect from tokenizer shifts on cross-lingual transfer. However, our results suggest that solving the embedding alignment problem is the "high-order bit" for crosslingual transfer: it has the largest impact on finetuning performance and is the least improved by continued pretraining. Thus, future progress on solving this problem in large-scale transformers may have outsized impact.

\section{Acknowledgements}

We would like to thank Nelson Liu, Mirac Suzgun and Tol\'{u}l\d{o}p\d{\'{e}} \`{O}g\'{u}nr\d{\`{e}}m\'{i} for useful discussions and comments on drafts.


\section{Introduction}

One of the key outstanding questions in the modern S-matrix program is understanding the renormalization group (RG) from an on-shell point of view. The Wilsonian picture is fundamentally Euclidean, and referring to ``coarse graining" in position space, or ``integrating out high momentum modes" in momentum space, crucially relies on concepts far removed from the on-shell philosophy. Alternatively, much of the technical content of the RG was earlier realized using the seemingly more complicated ideas of  Bogoliubov--Parasiuk--Hepp--Zimmermann \cite{10.1007/BF02392399,Hepp1966,Zimmermann1969} on ultraviolet (UV) divergences, which have been given a deeper Hopf-algebraic understanding \cite{Connes:1998qv,Connes:1999yr,Connes:2000fe} exhibiting a closer connection to on-shell ideas.  More recently,  analogous RG-like structures were found for infrared (IR) divergences in the context of QCD \cite{Mueller:1979ih,Collins:1980ih,Sen:1981sd,Korchemsky:1988si,Sterman:2002qn,Gardi:2009qi,Agarwal:2021ais} and later soft-collinear effective theory, see, e.g., \cite{Collins:1989gx,Becher:2009qa,StewartLectures,Becher:2014oda,Ma:2019hjq,Beekveldt:2020kzk}.

It is then natural to ask whether these perspectives on UV and IR divergences can be unified into a single picture, one which can connect the emerging geometric underpinnings of scattering amplitudes to the physics of UV divergences and the RG. In this note, we take a step in this direction, by describing a beautiful class of polytopes associated with any Feynman graph, whose facet structure concisely captures and calculates both the leading UV and IR divergences in a unified way.

Our starting point is the treatment of Feynman integrals in dimensional regularization from a worldline perspective, where instead of integrating over the loop momenta, the only variables are the Schwinger proper times $\alpha=e^\tau$ measuring lengths of each propagator. From this point of view, the limits $\tau \to \mp\infty$ correspond to shrinking and expanding edges, which are naturally associated to the UV and IR singularities respectively. They are thus put on the same footing and the unified analysis of all divergences boils down to examining the asymptotic behavior of the integrand at infinity. Mathematicians have studied this subject, known as \emph{tropical geometry}, for decades, see, e.g., \cite{maclagan2015introduction} for an introduction. We will explore the physical meaning and consequences of this geometry in our setting.

\begin{figure}
    \centering
    \includegraphics[]{figpdfs/fig1-figure0.pdf}
    \caption{\label{fig:scalings}Example scalings of Schwinger parameter responsible for: (a) collinear, (b) soft/collinear, (c) power-law IR, (d) nested soft/collinear, and (e) nested UV singularities in four dimensions. In each case the divergence comes from shrinking/expanding the relevant worldline edge lengths as the product of the variables $\rho,\sigma, \ldots$ to the indicated powers, as they all scale  $ \to \infty$.  Massless and massive particles are denoted with single and doubled lines respectively.}
\end{figure}

As we will see, all of the elaborate structure of UV and IR divergences of a Feynman integral associated with a graph $G$ is captured by the \emph{Feynman polytope} $\mathbf{S}_G$. It can be constructed out of the \emph{Symanzik polytopes} (previously explored in \cite{brown2011multiple,Brown:2015fyf,Schultka:2018nrs,Borinsky:2020rqs,Mizera:2021icv} for IR-safe kinematics) defined fully by the combinatorics of $G$ and the kinematic data of external momenta and particle masses.

Elements of polyhedral geometry of Feynman integrals previously appeared in sector decomposition \cite{Binoth:2000ps,Binoth:2003ak,Bogner:2007cr,Kaneko:2009qx,Kaneko:2010kj,Borowka:2015mxa,Schlenk:2016epj,Heinrich:2021dbf}, expansion by regions \cite{Pak:2010pt,Ananthanarayan:2018tog,Semenova:2018cwy,Ananthanarayan:2020fhl,Tellander:2021xdz}, and blowups of integration domains \cite{Bloch2006,Brown:2015fyf}.\footnote{Tropical geometry has also recently emerged as a useful tool in studying other, complementary, aspects of scattering amplitudes, see, e.g., \cite{Tourkine:2013rda,Panzer:2019yxl,Cachazo:2019ngv,Drummond:2020kqg,Arkani-Hamed:2020cig,Drummond:2019cxm,Arkani-Hamed:2019mrd,Arkani-Hamed:2019plo,Arkani-Hamed:2020tuz,Chicherin:2020umh,Cachazo:2020wgu,He:2021non}.}
Sector decomposition is a powerful general algorithm for studying singularities of generalized hypergeometric function (see, e.g., \cite{Mastrolia:2018uzb,delaCruz:2019skx,Klausen:2019hrg,Feng:2019bdx,Abreu:2019wzk,Tellander:2021xdz,Mizera:2021icv} for recent progress).  But as such, its application to the study of Feynman diagrams does not make use of any special properties that one might expect would arise for the integrals which arise from physics, meaning one has to laboriously compute $\mathbf{S}_G$ for every $G$ from scratch. One of our central results is that the structure of the integrals handed to us by physics, and the corresponding polyhedral geometry, is indeed very special, allowing us to give an analytic solution to this problem for any scalar Feynman diagram.

The polytopes begin their lives described as the convex hull of the exponent vectors of Symanzik polynomials. We show that they also have a beautifully simple \emph{facet description}, cut out by linear inequalities. It is this dual definition that makes it possible to understand and compute the leading UV and IR divergences.

An overarching theme in the emerging connection between combinatorial geometry and scattering amplitudes has been the realization that the singularities associated with ``factorization" of amplitudes are reflected in an analogous ``factorization" seen in the facet structure of the geometries. This is also seen in the Symanzik polytopes and ${\bf S}_G$, whose boundary structure transparently encodes the combinatorics of shrinking subgraphs of $G$, thereby precisely encoding the ``Russian doll" hierarchical separation of worldline edge-lengths in $G$ that have the potential to generate divergences.

The definition of ${\bf S}_G$ as cut out by inequalities further endows each face of $\mathbf{S}_G$ with a number: its UV/IR degree of divergence. The way these faces fit into the polytope dictates how the individual divergences are related to one another. In the UV, this prescription recovers the standard picture of nested and overlapping divergences, in agreement with \cite{Brown:2011pj,Brown:2015fyf}. In the IR, it reveals novel patterns of divergences, including the familiar soft/collinear divergences together with their intricate generalizations. Examples are shown in Fig.~\ref{fig:scalings}.

As an extension of \cite{Arkani-Hamed:2019mrd}, we explain how the leading divergent term of every Feynman integral can be computed in terms of volumes of dual cones. In this formalism, the divergence is viewed as an additional redundancy of the integrand in certain directions at infinity. We give a prescription for computing the coefficient of this divergence as a procedure of ``modding out'' by the enhanced symmetry. We apply this tool to various examples, including the rectangular fishnet diagrams, resulting in a formula for the leading IR divergence at any loop order.

Our goal in this note is to take a geodesic path to the central constructions, and illustrate their utility and power with a few examples. A more leisurely description of these ideas as well as further generalizations and examples will appear in \cite{AHHM}.

\section{Schwinger Parametrization}
A scalar Feynman integral in parameter space corresponding to a graph $G$ is given by
\begin{equation}\label{eq:I}
    \mathcal{I}_G = \Gamma(d) \int \frac{\mathrm{d}^\mathrm{E} \alpha}{\text{GL}(1)} \frac{1}{\mathcal{U}^{\mathrm{D}/2-d}\mathcal{F}^{d}},
\end{equation}
where $d = \mathrm{E}-\mathrm{L}\mathrm{D}/2$ is the superficial degree of divergence for a scalar diagram $G$ with $\mathrm{E}$ edges and $\mathrm{L}$ loops in dimension $\mathrm{D}$.  The functions in the integral are the Symanzik polynomials which admit a graphical expansion in terms of spanning trees $T$ and $2$-trees $T'$
\begin{equation}
    \mathcal{U} = \sum_T \prod\limits_{e \notin T}\alpha_e, \quad
    \mathcal{F} = \sum_{T'} p_{T'}^2 \!\!\prod\limits_{e \notin T'}\alpha_e - \mathcal{U} \sum_{e=1}^{\mathrm{E}} m_e^2\alpha_e,
\end{equation}
where $p_{T'}^\mu$ is the total momentum flowing across the $2$-tree, see, e.g., \cite{nakanishi1971graph}. Modding out by overall-scale $\mathrm{GL}(1)$ redundancy means we can set $\alpha_\mathrm{E}=1$. Including Feynman integrals with non-trivial numerators does not pose difficulty and will be considered in \cite{AHHM}.

Below, we will use $\gamma$ to mean a subgraph of $G$, and $G/\gamma$ to denote $G$ with each connected component of $\gamma$ contracted to a vertex. To avoid confusion, we add the subscript ${}_G$ to denote the quantity for a specific diagram $G$, e.g., $\mathrm{L}_\gamma$ denotes the number of loops in $\gamma$.

\section{\label{sec:general}Tropical Geometry}

Feynman integrals \eqref{eq:I} fall into a general class of integrals whose properties can be understood using tropical geometry \cite{gelfand2009discriminants,Nilsson2013,berkesch2013eulermellin,Arkani-Hamed:2019mrd}. In terms of the variables $\tau_e = \log\alpha_e$, the integrand of \eqref{eq:I} can be approximated at infinity, where the UV/IR divergences come from, by $e^{\mathrm{Trop}}$ with the piecewise-linear function
\begin{equation}\label{tropdef}
\mathrm{Trop} = \tau_G + (d_G -\tfrac{\mathrm{D}}{2})\max(\mathcal{U}_G) - d_G \max (\mathcal{F}_G),
\end{equation}
where $\tau_\gamma = \sum_{e \in \gamma}\tau_e$ and the maxima are taken over all the monomials in $\mathcal{U}_G$ and $\mathcal{F}_G$. Equivalently, one can capture the same information by the Minkowski sum/difference $\mathbf{S}_G = \mathbf{U}_G \oplus c \mathbf{F}_G$, of the Newton polytopes
\begin{equation}\label{eq:UF}
    \mathbf{U}_G = \mathrm{Newt}(\mathcal{U}_G), \qquad \mathbf{F}_G = \mathrm{Newt}(\mathcal{F}_G),
\end{equation}
which we refer to as \textit{Symanzik polytopes} and study closely in Sec.~\ref{sec:Symanzik}.
To obtain the most refined information about the regions where divergences can possibly come from, it is sufficient to consider $c>0$, see App.~\ref{app:general}.
Every \emph{ray} (boundary of domains of linearity) of $\mathrm{Trop}$ defines the normal to a facet of the Feynman polytope $\mathbf{S}_G$ and hence the two pictures can be used interchangeably.

\paragraph*{\bf Leading divergences.} Divergences come from cones (faces) generated by $\mathrm{Trop} \geqslant 0$ rays. We focus on logarithmic cases with $\mathrm{Trop} = 0$. Such $m$-dimensional cones (codimension-$m$ faces) give $\log^m$ or $1/\varepsilon^m$ UV/IR singularities. The leading divergence is controlled by cones $C_F$, or \emph{maximally-divergent faces} $F$ for which $m$ is the largest. As a refinement of the standard tropical rules \cite{Arkani-Hamed:2019mrd}, we find
\begin{equation}\label{tropvolume}
\mathcal{I}_G = \Gamma(d_G) \sum_{F} \mathrm{vol}(C_{F})\, \mathcal{J}_{F} + \ldots.
\end{equation}
It involves two ingredients: the volume of the cone $C_{F}$ (bounded by $\mathrm{Trop} \geqslant -1$) gives a numerical factor $\propto 1/\varepsilon^{m}$, while all the kinematic dependence is captured by the modded integral
\begin{equation}
\mathcal{J}_F =  \int \frac{\mathrm{d}^\mathrm{E} \alpha}{\text{GL}(1)^{m+1}}\, \frac{1}{\mathcal{U}_G^{\mathrm{D}/2-d}\mathcal{F}_G^{d}} \bigg|_{C_{F}^\perp}.
\end{equation}
Here we use the tropical approximation for the integrand in the direction $C_F^\perp$ orthogonal to the cone $C_F$, which by definition has an enhanced $\mathrm{GL}(1)^m$ symmetry. So extracting the coefficient of this divergence is precisely modding out by the extra redundancy!
Strictly speaking, the volume and the modded integral are not uniquely defined, but their product is, as illustrated on examples in App.~\ref{app:further-examples}.

\section{\label{sec:Symanzik}Symanzik Polytopes}

The polytopes $\mathbf{U}_G$ and $\mathbf{F}_G$ have been defined as a convex hull of vertices via the Newton polytopes \eqref{eq:UF}.  But the singularities of Feynman integrals are instead tied to their \emph{facet} description, cutting out the polytopes with linear inequalities.

Now, there is a systematic algorithm for finding the facets of an $n$-dimensional polytope given its vertices: checking whether a collection of vertices is a facet by asking whether the remaining vertices are all on the same side of the plane formed by this set. But clearly, the amount of work to do this grows without bound as the number of vertices and/or the dimension of the polytope become large.

Thus, what it means to ``understand" infinite families of polytopes, is to be able to analytically describe them in {\it both} the convex hull and facet descriptions. As we now show, the Symanzik polytopes can be ``understood" in this sense, and can be given an extremely simple facet presentation.

\paragraph*{\bf Inequalities.}
$\mathbf{U}_G$ is the Newton polytope of a homogeneous degree $\mathrm{L}_G$ polynomial and therefore is naturally defined on the hyperplane $a_G = \mathrm{L}_G$, where $\vec a$ is a point in the Newton polytope and $a_\gamma = \sum_{e \in \gamma}a_e$.
On this hyperplane, $\mathbf{U}_G$ is defined by
\begin{equation}
\label{udefn}
    \mathbf{U}_G\!:  \quad a_G = \mathrm{L}_G \quad\mathrm{and}\quad a_\gamma \geqslant \mathrm{L}_\gamma \quad\text{for all } \gamma\text{'s}.
\end{equation}
Analogously, $\mathbf{F}_G$ is cut out by
\begin{align}
\label{fdefn}
    \mathbf{F}_G\!:\quad  a_G = \mathrm{L}_G{+}1  \quad\text{and}\;\, &\begin{cases}
    a_\gamma \geqslant \mathrm{L}_\gamma & \text{if}\quad\mathcal{F}_{G/\gamma} \neq 0,
    \\
    a_\gamma \geqslant \mathrm{L}_\gamma{+}1 &  \text{if}\quad\mathcal{F}_{G/\gamma} = 0,
    \end{cases}\nonumber\\
    & \qquad\qquad\qquad\;\text{for all } \gamma\text{'s}.
\end{align}
Graphs $G/\gamma$ with $\mathcal{F}_{G/\gamma} = 0$ are called \emph{scaleless} and define inequalities crucial for the study of IR divergences. In contrast, for generic kinematics (say, all particles massive), the definition of $\mathbf{F}_G$ only differs from that of $\mathbf{U}_G$ by the choice of hyperplane and indeed is the Minkowski sum $\mathbf{F}_G = \mathbf{U}_G \oplus {\bf\Delta}_{\mathrm{E}-1}$, where ${\bf\Delta}_{\mathrm{E}-1}$ is the $(\mathrm{E}{-}1)$-dimensional simplex.

\paragraph*{\bf Boundary structure.}

Saturating the above inequalities defines faces of $\mathbf{U}_G$ and $\mathbf{F}_G$ with the following factorization properties:
\begin{equation}
\partial\mathbf{U}_G \supset \mathbf{U}_\gamma \times \mathbf{U}_{G/\gamma} \quad\text{on}\quad a_\gamma = \mathrm{L}_\gamma
\end{equation}
and
\begin{align}\label{eq:del-F}
    \partial\mathbf{F}_G \supset
    \begin{cases}
    \mathbf{U}_\gamma \times \mathbf{F}_{G/\gamma} & \text{on}\quad a_\gamma = \mathrm{L}_\gamma,\\
    \mathbf{F}_\gamma \times \mathbf{U}_{G/\gamma} &
    \text{on}\quad a_\gamma = \mathrm{L}_\gamma{+}1.
    \end{cases}
\end{align}
It follows that $\mathbf{S}_G$ factorizes in the same way as $\mathbf{F}_G$. Furthermore, we claim that \emph{facets} (codimension-$1$ faces) of $\mathbf{S}_G$ are those faces of \eqref{eq:del-F} for which every instance of $\mathbf{U}_{\gamma'}$ in the factorization above has $\gamma'$ 1VI (recall that a diagram is called 1VI if it cannot be disconnected by removal of a single vertex) and, in the $\mathcal{F}_{G/\gamma} = 0$ case, no subgraph $\tilde \gamma \subset \gamma$ already has $\mathcal{F}_{G/\tilde \gamma} = 0$ and $\mathrm{L}_{\tilde\gamma} = \mathrm{L}_\gamma$.

\paragraph*{\bf Compatibility of facets.} Feynman integrals motivate us to study the properties of $\mathbf{S}_G$.  In particular, it is natural to ask if two facets labelled by $\gamma_1$ and $\gamma_2$ are compatible, i.e., they meet.  The following identity
\begin{align}
    a_{\gamma_1}+a_{\gamma_2} = a_{\gamma_1 \cup \gamma_2}+a_{\gamma_1 \cap \gamma_2}
\end{align}
furnishes a necessary condition for compatibility:
\begin{multline}
\label{compatibility}
     (1+c)(\mathrm{L}_{\gamma_1}+\mathrm{L}_{\gamma_2}) +(c)_{\gamma_1}+(c)_{\gamma_2} \geqslant\\
     (1+c)(\mathrm{L}_{\gamma_1\cup\gamma_2}+\mathrm{L}_{\gamma_1\cap\gamma_2}) +(c)_{\gamma_1\cup\gamma_2}+(c)_{\gamma_1\cap\gamma_2}.
\end{multline}
The $c$'s in the parentheses are only included when the $\mathbf{F}_G$ facet is of $a_\gamma \geqslant \mathrm{L}_\gamma{+}1$ type.

Note that the above condition means Feynman polytopes are in general \emph{not} generalized permutohedra, in contrast with \cite{Brown:2015fyf,Schultka:2018nrs,Borinsky:2020rqs}.

\section{Tropicalization}

Recall that rays in the Trop space are dual to facets of the corresponding polytope.
Given the explicit form of the $\text{Trop}$ function \eqref{tropdef}, one can simply plug in the ray corresponding to a given inequality and get the associated value of $\text{Trop}$ along that ray.  However, equipped with the explicit \textit{constants} in the inequalities for $\mathbf{U}_G$ and $\mathbf{F}_G$, we already know the value of Trop along the facet rays: they are given by the constant in the inequality.

\paragraph*{\bf Degrees of divergence.} The two polytopes are defined by intersecting the same inequalities with different subspaces.  Because of projective invariance, we can uniformly shift the points on $\mathbf{U}_G$ and $\mathbf{F}_G$, which has no effect on the $\mathrm{Trop}$ function, with
\begin{align}
    \label{hyperplaneshifts}
    \vec a \to  \vec a- \tfrac{\mathrm{L}_G}{\mathrm{E}_G}\vec 1_G, \hspace{7mm}  \vec a \to  \vec a-\tfrac{\mathrm{L}_G+1}{\mathrm{E}_G}\vec 1_G,
\end{align}
respectively,
where $\vec 1_G$ is a vector with all entries equal to $1$. At the level of the integrand, this is equivalent to distributing the overall $\alpha_1\alpha_2\cdots \alpha_\mathrm{E}$ between the $\mathcal{U}_G$ and $\mathcal{F}_G$.  Note that no change of integration variables is implied here, but the inequalities in \eqref{udefn} and \eqref{fdefn} are shifted.
The resulting inequalities for $\mathbf{U}_G$ are:
\begin{align}
    \label{shiftedineqsU}
    a_\gamma  \geqslant \frac{\mathrm{L}_\gamma \mathrm{E}_G-\mathrm{L}_G \mathrm{E}_\gamma}{\mathrm{E}_G}
\end{align}
and for $\mathbf{F}_G$:
\begin{align}
\label{shiftedineqsF}
a_\gamma \geqslant
   \begin{cases}
   \frac{\mathrm{L}_\gamma \mathrm{E}_G-(\mathrm{L}_G+1)\mathrm{E}_\gamma}{\mathrm{E}_G}& \text{if}\quad\mathcal{F}_{G/\gamma} \neq 0,\\
    \frac{(\mathrm{L}_\gamma+1)\mathrm{E}_G-(\mathrm{L}_G+1)\mathrm{E}_\gamma}{\mathrm{E}_G}&\text{if}\quad\mathcal{F}_{G/\gamma} = 0,
   \end{cases}
\end{align}
for the two types of $\mathbf{F}_G$ facets.  The value of Trop along a ray is given by substituting the saturation of that ray's inequality in the maxima. Since our conventions identify the value of Trop with the constant in the outward-pointing normal inequality, we substitute the constants in \eqref{shiftedineqsU} and \eqref{shiftedineqsF} with a minus sign, see App.~\ref{app:general}. This allows us to compute the value of Trop along any extremal ray.  The two types of faces are distinguished by their inequality in $\mathbf{F}_G$ and give two cases:
\begin{align}\label{eq:Trop-gamma}
    \text{Trop} = \begin{cases}
    -d_\gamma & \text{if}\quad \mathcal{F}_{G/\gamma} \neq 0, \\
    d_{G/\gamma} & \text{if}\quad \mathcal{F}_{G/\gamma} = 0.
    \end{cases}
\end{align}
The first case is familiar and associated with superficially UV divergent subgraphs when $d_\gamma \leqslant 0$.  We will see that the second case is associated with IR divergences when $d_{G/\gamma} \geqslant 0$ and $G/\gamma$ is scaleless.

\section{Generalized Weinberg's theorem}
The above geometric picture makes Weinberg's theorem \cite{PhysRev.118.838} obvious: if all subdiagrams $\gamma \subset G$ have $d_{\gamma}>0$ (with appropriately defined $d_\gamma$ for non-scalar diagrams), the diagram $G$ is UV convergent.
With \eqref{eq:Trop-gamma}, we are now equipped to state its IR counterpart: if all subdiagrams $\gamma \subset G$ for which $G/\gamma$ is scaleless have $d_{G/\gamma} < 0$, the diagram $G$ is IR finite. Together, the two theorems give conditions for any diagram to be UV/IR finite, cf. \cite{Lowenstein:1975rg}.

\section{\label{sec:UV} Ultraviolet}
\begin{figure}[!t]
	\centering
	\includegraphics[]{figpdfs/fig71-figure0.pdf}
	\quad
    \raisebox{-.15\height}{\includegraphics[]{figpdfs/fig72-figure0.pdf}}
	\caption{\label{eyeofsauron}Two UV-divergent Feynman diagrams $G_L$ and $G_R$ considered in Sec.~\ref{sec:UV}.
	}
\end{figure}

We now provide some illustrative examples.  As a UV example, we consider the diagram $G_L$ in Fig.~\ref{eyeofsauron} (left). It contributes a leading logarithmic divergence in $\mathrm{D} = 4-2\varepsilon$ and therefore has a two-dimensional divergent Trop = 0 space (recall that one divergence comes from the overall $\Gamma(d_{G_L}) = \tfrac{1}{3\varepsilon} + \dots$).  In order to calculate the leading divergence, we first need to identify the divergent rays spanning the Trop = 0 space.  These correspond to the superficially divergent subgraphs, of which there are three: $\gamma_{1256}$, $\gamma_{3456}$, $\gamma_{56}$ where the subscripts denote the edges of the subgraph.  Next, we need to know the compatibility of these divergent rays in $\mathbf{S}_{G_L}$.  This is given by the compatibility criterion (\ref{compatibility}) and demonstrates that $\gamma_{56}$ is compatible with each of the other two subgraphs, but $\gamma_{1256}$ and $\gamma_{3456}$ are not compatible with each other.  This means that we have two two-dimensional cones and they share a ray.  We must compute the contribution to the leading divergence from each two-dimensional cone separately and add them.

We may now leverage eq.~(\ref{tropvolume}) to obtain
\begin{align}
    \label{eyedivergence}
    \mathcal{I}_{G_L} = \Gamma(d_{G_L}) \frac{1}{d_{\gamma_{56}}}\left(\frac{1}{d_{\gamma_{1256}}}+\frac{1}{d_{\gamma_{3456}}}\right)+\ldots,
\end{align}
where the finite integral is accounted for as it is merely a product of bubble integrals which equal one. Since $d_{\gamma_{1256}} = d_{\gamma_{3456}} = 2\varepsilon$ and $d_{\gamma_{56}} = \varepsilon$, the leading UV divergence contributes $\mathcal{I}_{G_L} = \frac{1}{3\varepsilon^3} + \ldots$.  The structure of the rays in this example captures the familiar RG intuition of sequentially shrinking loops.  One simplicial cone is generated by rays $r_{1256}$ and $r_{56}$ and captures the shrinking of the two corresponding subgraphs, one nested in the other.  Any ray in this cone corresponds to a direction of scaling integration variables where the bubble $\gamma_{56}$ shrinks at least as fast as the collective $\gamma_{1256}$, capturing the physics of shrinking successive UV subgraphs to effective vertices.

\begin{figure}[!t]
    \centering
    \includegraphics[]{figpdfs/fig81-figure0.pdf}
	\includegraphics[]{figpdfs/fig82-figure0.pdf}
    \caption{Divergent rays of the diagrams from Fig.~\ref{eyeofsauron}. Left: Two two-dimensional cones which share a single ray $r_{56}$, reflecting the fact that there is an overlapping divergence. Right: Triangulation of the three-dimensional cone into two simplicial cones separated by the red plane.  In this triangulation each simplex has the interpretation of nested shrinking subgraphs with ordering implied by the rays.}
    \label{eyeinbubblerays}
\end{figure}

\paragraph*{\bf Leading divergences and the RG.}
In general, the maximal cones Trop = 0 cones need not be simplicial. A given cone's contribution then need not be interpreted in the familiar nested shrinking picture consistent with our familiar RG intuitions.  Nonetheless, there is a triangulation of the cone that does this.  Consider the diagram $G_R$ from Fig.~\ref{eyeofsauron} (right), whose rays are depicted in Fig.~\ref{eyeinbubblerays} (right).  The triangulation divides the cone into two simplices, with rays that again imply an ordering of scaling rate.  That is, any ray in the simplex generated by $r_{345678}, r_{3478}, r_{78}$ has components $\gamma_{78}$ shrink the fastest, followed by $\gamma_{34}$ and then $\gamma_{56}$.  The nested shrinking is therefore one natural way of organizing the calculation of the leading the divergence, though not uniquely so.  Using \eqref{tropvolume}, we can verify that the two triangulations produce the same leading contribution:
\begin{align}
    &\frac{1}{d_{\gamma_{345678}}d_{\gamma_{78}}}\left(\frac{1}{d_{\gamma_{3478}}}+\frac{1}{d_{\gamma_{5678}}}\right)\nonumber\\ =&\,\frac{1}{d_{\gamma_{3478}}d_{\gamma_{5678}}}\left(\frac{1}{d_{\gamma_{78}}}+\frac{1}{d_{\gamma_{345678}}}\right) = \frac{1}{3\varepsilon^3},
\end{align}
and hence the leading contribution in this case is $\mathcal{I}_{G_R} = \Gamma(d_{G_R}) \frac{1}{3\varepsilon^3} + \ldots = \frac{1}{12\varepsilon^4}+\dots$.
Similarly, the phenomenon of overlapping divergences for the diagram $G_L$ is illustrated in Fig.~\ref{eyeinbubblerays} (left), where the ray $r_{56}$ is shared between two two-dimensional cones.

The fact that leading UV divergences are organized in a ``Russian doll" structure of shrinking nested subgraphs, reflects the Wilsonian way of organizing the generation of logarithmic divergences scale-by-scale. This is a simple but important first step in understanding the origin of the renormalization group. Of course the full structure of the RG goes beyond this basic fact about leading divergences of individual diagrams, and organizes the UV divergences into leading, subleading etc. logarithms captured by running coupling constants. We will discuss a geometric perspective on this understanding in \cite{AHHM}.

\section{\label{sec:IR} Infrared}

\begin{figure}[!t]
    \centering
    \includegraphics[scale=.7]{figpdfs/fig2-figure0.pdf}
    \caption{Left: Collinear divergence illustrated for the box diagram $G_B$ with one massless corner (single line) and the remaining ones massive (doubled line). A subgraph $\gamma_{14}$ corresponding to the graph's unique divergent ray in $\mathrm{D} = 4-2\varepsilon$.  Right: The resulting $G_B/\gamma_{14}$ is scaleless.}
    \label{IRcorner}
\end{figure}

Combinatorics of leading divergences is much richer in the IR.
The rays for which $\text{Trop} = -d_{G/\gamma}$ correspond to scaleless $G/\gamma$'s.  We are particularly interested in $G/\gamma$'s for which  $d_{G/\gamma} = 0$, as they correspond to logarithmic IR divergences.
Such $G/\gamma$'s can then be blown up into a graph $G$ by inserting a $\gamma$ into it. By studying the possible structures of $G/\gamma$'s, we can understand the ways in which IR divergences can appear.

All scaleless graphs must have massless internal edges.  There are then two ways to make scaleless graphs: momentum conservation or masslessness.  Either the only spanning $2$-trees isolate a single massless external momentum or they isolate all external momenta and vanish by momentum conservation.  In order to have a divergence, the sign of Trop tells us that  we want $d_{G/\gamma}$ to be as convergent as possible; the more superficially convergent $d_{G/\gamma}$, the more divergent the contribution to the Feynman integral $\mathcal{I}_G$.  But we also need $G/\gamma$ to be scaleless. There is a basic tension between being scaleless and being superficially convergent.  We can make a graph arbitrarily superficially convergent by adding external lines, i.e., bivalent vertices in the amputated graph, but we cannot have more than two such vertices with external kinematics and be scaleless.

\paragraph*{\bf In which dimensions do IR divergences exist?}
This discussion allows us to bound the dimensions in which IR divergences appear.  In particular, we can ask: what is the highest dimension in which a marginal IR divergence may exist?  For $n > 3$, we will maximize $\mathrm{E}-\mathrm{L} \mathrm{D}/2$ when we have a bubble filled with trivalent vertices.  In this case, with $2v$ trivalent vertices, we have
\begin{align}
    d_{G/\gamma} = \frac{(4-\mathrm{D}) + v(6-\mathrm{D})}{2}.
\end{align}
We see that as $v\to \infty$ then $\mathrm{D} \to 6$ becomes marginal, therefore one will not find IR divergences in six or higher dimensions.
Since adding numerators can only soften IR divergences, these bounds apply to a general theory.
Fig.~\ref{IRdimensions} catalogues the leading examples of marginal scaleless graphs in the relevant dimensions.

It is worth noting that the existence of marginal IR divergences in a given dimension implies the existence of power IR divergences in lower dimensions.
In particular, scalar theories such as $\phi^3$ exhibit power IR divergences in $\mathrm{D}=4$ \cite{AHHM}, see Fig.~\ref{fig:scalings}c, giving an additional raison d'être for numerators in gauge theories: softening IR divergences to logarithmic.

\begin{figure}[!t]
    \centering
    \includegraphics[]{figpdfs/fig9-figure0.pdf}
    \caption{Leading $G/\gamma$ graphs which are scaleless and marginal in the relevant dimensions.  We write $\mathrm{D}\to 6$ as marginality is only attained in the limit of an infinite number of rungs.}
    \label{IRdimensions}
\end{figure}

\paragraph*{\bf Soft/collinear divergences.}

The standard collinear divergence comes from a ray $r_{14}$ corresponding to shrinking $\gamma_{14}$ (or expanding $\gamma_{23}$), as illustrated in Fig.~\ref{fig:scalings}a. It is singular because $G_B/\gamma_{14}$ is scaleless, see Fig.~\ref{IRcorner}.
The cone is therefore a single ray $r_{14}$, which according to \eqref{tropvolume} gives, in the conventions of Fig.~\ref{IRcorner}
\begin{align}\!
\vcenter{\hbox{\includegraphics[width=7mm]{figpdfs/inline1-figure0.pdf}}}\!\!&=\!  -\frac{\Gamma(d_{G_B}\!)}{d_{G_{\!B}\!/\gamma_{14}}}\!\! \int\!\!\! \frac{\mathrm{d}^4 \alpha/\text{GL(1)}^2}{(s\alpha_1\alpha_3{+}t\alpha_2\alpha_4{+}p_1^2\alpha_1\alpha_2{+}p_3^2\alpha_3\alpha_4)^2} \!+\!\dots\nonumber\\
&= -\frac{1}{\varepsilon}  \int_{0}^{\infty}\!\! \frac{\mathrm{d} \alpha_1 \mathrm{d} \alpha_2 }{(s \alpha_1+t\alpha_2+p_1^2\alpha_1\alpha_2+p_3^2)^2} + \dots,\nonumber\\
&= -\frac{1}{\varepsilon}  \frac{\log\left(\frac{st}{p_1^2 p_3^2}  \right)}{s t-p_1^2p_3^2 } +\dots,
\label{eq:one-corner-box}
\end{align}
where we used the two $\mathrm{GL}(1)$'s to fix $\alpha_3 = \alpha_4 = 1$ and $d_{G_B} = 2{+}\varepsilon$ with $d_{G_B/\gamma_{14}} = \varepsilon$,
see App.~\ref{app:further-examples} for more details.

The soft/collinear divergence is illustrated in Fig.~\ref{fig:scalings}b (see also \cite{YelleshpurSrikant:2019khx}) and comes from a two-dimensional cone generated by two IR rays, say $r_{14}$ and $r_{34}$, sharing an edge. Its leading divergence is
\begin{align}
\vcenter{\hbox{\includegraphics[width=7mm]{figpdfs/inline2-figure0.pdf}}} &= \frac{\Gamma(d_{G_B})}{d_{G_B\!/\gamma_{14}}\, d_{G_B\!/\gamma_{34}}} \int  \frac{\mathrm{d}^4 \alpha/\text{GL(1)}^3}{(s\alpha_1\alpha_3+t\alpha_2\alpha_4)^2} + \dots \nonumber\\
&= \frac{1}{\varepsilon^2 st} +\dots,
\end{align}
where we used the enhanced $\mathrm{GL}(1)^3$ to fix $\alpha_2=\alpha_3=\alpha_4=1$. Generalizations to nested IR divergences can be quite intricate, see, e.g., Fig.~\ref{fig:scalings}d.

\paragraph*{\bf Fishnet diagrams.}

For scalar integrals in $\mathrm{D}=4$, it is natural to study the generalization to \emph{fishnet} diagrams, i.e., planar arrays of $M \times N$ squares \cite{Gurdogan:2015csr,Chicherin:2017frs,Loebbert:2020tje}, see Fig.~\ref{fishnetexample} for $M=N=4$ with two massless corners.
As discussed in Sec.~\ref{sec:general},
calculation of the divergent contribution from a graph is performed by identifying the rays generating divergent cones and applying \eqref{tropvolume}.  In the case of IR-divergent fishnets, the combinatorics of these rays is rich, but admits a simple description.  It suffices to study the case with two diagonally-opposed massless corners, as the combinatorics of IR rays associated with neighboring corners is trivial, i.e., the geometry is a product.  This means that for any assignment of the external masses, we can treat the diagonals independently and multiply the result.  A single leading-dimensional cone of the Trop = 0 space is in correspondence with a staircase walk (dashed lines) separating the fishnet into two Young tableau-shaped sets of edges (blue and red) as in Fig.~\ref{fishnetexample}.
Then the compatible rays $\gamma$ are given by the complements of
all the Young \emph{subdiagrams} based at the top-left (blue) or bottom-right corner (red), including the dashed bordering edges.  For example, the blue shading with the dashed lines is the $\gamma$ corresponding to the tableaux given by the red shading.

\begin{figure}[!t]
    \centering
    \scalebox{0.5}{
    \includegraphics[]{figpdfs/fig31-figure0.pdf}}
    \hspace{1.25cm}\scalebox{1.15}{
   \includegraphics[]{figpdfs/fig32-figure0.pdf}}
    \caption{Left: Half of a top-dimensional contribution to a four-point fishnet diagram with opposite massless corners.  The IR subgraphs (red and blue plus respective sub-tableaux) are complementary to those appearing in the inequalities, i.e., these are the edges becoming long/soft.  We additionally label for the red portion one of the standard tableaux corresponding to a cell in the triangulation. Right: Scaling parameters $\rho_i \to \infty$ of the lower-right edges for a square labeled by integer $i$ in the standard tableaux.}
    \label{fishnetexample}
\end{figure}

Since such cones are in general not simplicial, in order to leverage \eqref{tropvolume} one needs a triangulation.  For concreteness, consider the red subdiagram in Fig.~\ref{fishnetexample}.  One triangulation is given by cells corresponding to the \emph{standard} tableaux, see Fig.~\ref{fishnetexample}. Scaling of the bottom and right edges of every red box labeled with integer $i$ is given by $\rho_i \rho_{i+1} \cdots \rho_6$, where each $\rho_i$ parameterizes scaling along a ray generating the simplicial cell.  Each cell clearly contributes the volume $\frac{1}{6!}$, so the non-trivial combinatorics is in the counting of the standard tableaux, given by the hook-length formula.  We claim that summing over the volume contributions from all the cells yields
\begin{equation}
    \label{catalan}
    \frac{2^{MN}}{(MN)!} C_{M,N} = \!\!\sum\limits_{\lambda_{i} > \lambda_{i+1}}\frac{f^{(\lambda_1, \dots, \lambda_M)}f^{(\bar \lambda_M, \dots, \bar \lambda_1)}}{(\lambda_1{+}\dots{+} \lambda_M)!(\bar\lambda_1{+}\dots{+} \bar\lambda_M)!}  ,
\end{equation}
where $\bar \lambda_i = N-\lambda_i$ and $f^{(\lambda_1, \dots, \lambda_M)}$ counts the number of the standard Tableaux with shape $\lambda = (\lambda_1, \dots, \lambda_M)$, and $C_{M, N}$ is the multi-dimensional Catalan number \cite{MultiCatalan} with standard definition:
\begin{equation}
    C_{M,N} = (M N)!\frac{\text{sf}(M-1)\text{sf}(N-1)}{\text{sf}(M+N-1)}
\end{equation}
with the superfactorial $\text{sf}(N) = 0!1!2!\cdots N!$.

For the fishnet with one massless corner, there is no sum over partitions and we merely have to count the number of standard tableaux for the $M \times N$ rectangle, which is given by $C_{M, N}$:
\begin{equation}
\label{3massfish}
\raisebox{-2em}{\includegraphics[]{figpdfs/fig3massfish-figure0.pdf}}
\;=\; \frac{F_{M,N}(s, t, p_i^2)}{(MN)!} \frac{C_{M,N}}{\varepsilon^{MN}} + \ldots,
\end{equation}
where $F_{M,N}(s, t, p_i^2)$ is a non-trivial function of kinematics computed by the $\text{GL}(1)^{MN}$-modded integral in the way described above, see \eqref{eq:one-corner-box} for a special case. And in the case of two diagonally opposed massless edges, we employ \eqref{catalan} yielding:
\begin{equation}
\label{2massfish}
\raisebox{-2em}{\includegraphics[]{figpdfs/fig2massfish-figure0.pdf}}
\;=\; \frac{2^{MN}}{(MN)!} \frac{C_{M, N}}{\varepsilon^{MN} s^N t^M} + \ldots,
\end{equation}
The leading contribution to the fully massless fishnet is obtained by squaring the kinematic-independent part of \eqref{2massfish} and keeping the same kinematic dependence.\footnote{Similarly, the contribution with two neighboring massless corners is obtained by squaring the kinematic-independent part of \eqref{3massfish} and doing the appropriate modded integral giving $1/(s^N t^M)$. The contribution with the three massless corners is obtained by multiplying kinematic-independent parts of \eqref{3massfish} and \eqref{2massfish} and again obtaining $1/(s^N t^M)$ from the modded integral.} See \cite{AHHM} for details and \cite{Basso:2021omx,Olivucci:2021pss} for alternative perspectives.

\section{Outlook}

In this work we introduced polytopes capturing the structure of UV/IR divergences of Feynman integrals. Recently, Symanzik polytopes also played a role in the study of kinematic singularities, known as anomalous thresholds \cite{Mizera:2021icv}, and it would be interesting to relate the two approaches. Other future direction include a systematic approach to extracting subleading divergences, inclusion of kinematic numerators, as well as summing over multiple Feynman diagrams under the same integral, further pursued in \cite{AHHM}.

\begin{acknowledgments}
We thank Francis Brown for introducing us to Symanzik polytopes and for sharing his unpublished notes \cite{BrownIHES}. We also thank Nick Early, Hofie Hannesdottir, Song He, Thomas Lam, Erik Panzer,  Giulio Salvatori, and Simon Telen for stimulating discussions. This work is supported by the grant DE-SC0009988 from the U.S. Department of Energy.
S.M. gratefully acknowledges the funding provided by Frank and Peggy Taplin.
\end{acknowledgments}


\section{Introduction}

Emission spectroscopy of exoplanets allows us to probe into the deeper layers (1-10 bar) of the dayside atmosphere compared to transmission spectroscopy which measures the upper layers ($\sim$1mbar) of the planetary limbs. The difference in the radiative transfer path geometry means the emission spectrum is much more sensitive to the vertical thermal structure of the planet. Depending on the atmospheric chemical composition and the temperature-pressure profile, molecular absorption or emission features can be imprinted onto the emission spectrum. Notably, water features at 1.4 $\mu m$ \citep{fu_statistical_2017} has been the most robustly detected in both absorption \citep{kreidberg_precise_2014} and emission \citep{fu_strong_2022,evans_ultrahot_2017}. The presence of H- continuum opacity was inferred in multiple ultra-hot Jupiters \citep{arcangeli_climate_2019, fu_hubble_2021-1}. In the longer wavelength, the CO/CO$_2$ features were indicated by the deviations of Spitzer photometric points in 3.6 $\mu m$ and 4.5 $\mu m$ channels from the blackbody approximation \citep{garhart_statistical_2020, fu_hubble_2021}.

HAT-P-41b is an inflated hot Jupiter (R=$1.685^{+0.076}_{-0.051}$R$_{Jup}$ M=$0.8^{+0.1}_{-0.1}$M$_{Jup}$ T$_{eq}$=$1940^{+38}_{-38}$K) \citep{hartman_hat-p-39bhat-p-41b_2012} with detailed atmospheric characterization via transmission spectroscopy \citep{sheppard_hubble_2021, wakeford_into_2020, lewis_into_2020}. The planet shows a high metallicity \citep{wakeford_into_2020, sheppard_hubble_2021} atmosphere with increased H- opacity abundance \citep{lewis_into_2020}. With a dayside temperature of $\sim$2300K, HAT-P-41b sits between hot and ultra-hot Jupiters where physical processes such as molecular dissociation and H- opacity are starting to become important \citep{parmentier_thermal_2018}. HAT-P-41b allows us to pinpoint the transitional parameter space of various atmospheric processes which makes it a valuable data point in our understanding of hot Jupiter atmospheres on a population level \citep{mansfield_unique_2021, baxter_evidence_2021}. Here we present the 1.1 to 4.5 $\mu m$ emission spectrum of HAT-P-41b with combined data from HST/WFC3 G141 and $Spitzer$ channel 1 and 2.

\begin{figure*}
\centering
  \includegraphics[width=\textwidth,keepaspectratio]{transit_wfc3_8pix_r1.png}
  \caption{HST WFC3 G141 spectral bin transit lightcurves after ramp-effect correction using \texttt{RECTE} (left) and corresponding residuals (right) with their relative ratios to the photon-limit noise ($\sigma/\sigma_\gamma$) levels.}
  \label{fig:WFC3}
\end{figure*}

\begin{figure}
\centering
  \includegraphics[width=0.5\textwidth,keepaspectratio]{mansfield_compare.png}
  \caption{The WFC3/G141 emission spectrum is in excellent agreement with \cite{mansfield_unique_2021}.}
  \label{fig:mansfield}
\end{figure}

\section{Observations}

The $Hubble$ WFC3/G141 dataset was observed as a part of the Panchromatic Exoplanet Treasury Survey (PanCET GO 14767, PIs: Sing $\&$ Lopez-Morales) on 2016-10-09 and the $Spitzer$ dataset was published in \citep{garhart_statistical_2020}. The WFC3/G141 eclipse observation was taken in spatial scan mode for five consecutive HST orbits. Each frame was taken with the 256 $\times$ 256 pixel subarray in the SPARS10 and NSAMP=12 settings. The forward scanning rate is 0.065 arcsec $s^{-1}$ and the exposure time is 81 seconds.

\section{HST WFC3 data reduction}

All the orbital parameters for the WFC3/G141 data reduction have been fixed to the same values used in \citep{sheppard_hubble_2021} which are identical to what were used in the \citep{garhart_statistical_2020} $Spitzer$ analysis.

The data reduction process starts with applying the standard flat field, background subtraction and bad pixels, cosmic rays removal on the ima frames \citep{fu_hubble_2021-1}. Then we extract the non-destructive reads from each frame \citep{deming_infrared_2013}. There is a spatially resolved companion star \citep{sheppard_hubble_2021} located 3.6 arcseconds away. Due to the large spatial separation, the non-destructive reads of spatial scan from the two stars do not overlap which allows for clean removal of the companion star spectra. The companion star removed reads are then combined to the complete spatial scan frames. Next we summed each frame in the vertical direction for the 1D spectrum and normalized it by its own median flux to calculate the wavelength shifts. We then used $scipy.interpolate.interp1d$ function to interpolate the 1D spectrum in the wavelength direction and calculated the relative sub-pixel level shifts of each frame. The largest shift between any two frames is under 0.1 pixel which does not induce any excessive systematics \citep{stevenson_analyzing_2019}. Wavelength shifts corrected non-normalized frames were then summed in all wavelength channels to form the whitelight eclipse lightcurve. The lightcurve is then fitted with a combination of \texttt{BATMAN} \citep{kreidberg_batman_2015} with the \texttt{RECTE} charge trapping systematics model \citep{zhou_physical_2017}, 2nd order polynomial of the HST orbital phase and the wavelength shifts. Each wavelength channel is then fitted using the same method but with the mid-transit time fixed to the best-fit whitelight value. The best-fit lightcurves of each wavelength channel and corresponding residuals are shown in Figure \ref{fig:WFC3}.

\begin{figure*}
\centering
  \includegraphics[width=\textwidth,keepaspectratio]{atmo_retrieval.pdf}
  \caption{ATMO retrieval of HAT-P-41b emission spectrum. The featureless black-like spectrum (left) is best fitted with an isothermal TP profile (right).}
  \label{fig:atmo_retrieval}
\end{figure*}

\begin{figure}
\centering
  \includegraphics[width=0.45\textwidth,keepaspectratio]{arc_monitor_nested_corner.pdf}
  \caption{Posterior distribution of the ATMO retrieval. The metallicity is mostly unconstrained due to a lack of features within the emission spectrum.}
  \label{fig:atmo_corner}
\end{figure}

\begin{figure*}
\centering
  \includegraphics[width=\textwidth,keepaspectratio]{compare.png}
  \caption{HAT-P-41b WFC3/G141 emission spectrum compared with other hot Jupiters. The water feature strength index (S$_{H_2O}$) \citep{mansfield_unique_2021} measures the 1.4$\mu m$ water band feature size relative to the blackbody model based on the out-of-band spectrum. Centered at zero being the blackbody-like featureless spectrum. Increasing negative values (redder color) represent stronger emission, features and increasing positive values (bluer color) represent stronger absorption features. The circle size scales with the host star temperature (T$_{eff}$). The grey vertical dash lines are condensation lines for various metals at 100 mbar assuming solar metallicity. We see mostly featureless spectra at higher dayside temperatures ($>$2700K) where H$^-$ continuum opacity expects to dominate \citep{parmentier_thermal_2018}. At the lower temperature side, we see mostly water absorption features where the atmosphere is too cool for gaseous heavy metal absorbers such as TiO to form and drive thermal inversion. HAT-P-41b sits at this in-between transitional region where planets can have absorption, emission or featureless spectra. Further observations of HAT-P-41b in the longer wavelength with JWST will help us to precisely determine its thermal structure and understand what causes thermal inversion in hot Jupiters atmospheres.}
  \label{fig:compare}
\end{figure*}

\section{Retrieval analysis}

We performed the retrieval analysis on the reduced HAT-P-41b emission spectrum (Table. \ref{eclipse_spectrum}) with \texttt{ATMO} \citep{amundsen_accuracy_2014, drummond_effects_2016, goyal_library_2018, tremblin_fingering_2015, tremblin_cloudless_2016, wakeford_hat-p-26b_2017}. The retrieval setup is equilibrium chemistry with a fixed solar C/O ratio. There is a total of 4 free fitting parameters (Table. \ref{eclipse_priors}) including metallicity (Z) and three parameters $\kappa_{IR}$, $\gamma_{O/IR}$ and $\beta$ for the temperature-pressure (TP) profile as defined in \cite{line_systematic_2013}. We fixed the C/O ratio to the solar value due to the lack of constrain this dataset has on the carbon-bearing species. Also, the transit spectra retrieval \citep{sheppard_hubble_2021, lewis_into_2020} results were all consistent with the solar C/O values.

\begin{center}
\begin{deluxetable}{ccc}
\tablecaption{\textbf{ATMO eclipse retrieval priors $\&$ posteriors}}
\tablehead{\colhead{Parameter} & \colhead{Priors} & \colhead{Posteriors}}
\startdata
\hline\hline
log(Z/Z$_{\Sun}$)	&	$\mathcal{U}$(-2.8, 2.8)	&	$-0.342^{+1.677}_{-1.277}$	\\
log(K$_{IR}$)	&	$\mathcal{U}$(-5, -0.5)	&	$-2.887^{+0.982}_{-0.907}$	\\
log($\gamma$/IR)	&	$\mathcal{U}$(-4, 1.5)	&	$-0.907^{+0.672}_{-1.716}$	\\
beta	&	$\mathcal{U}$(0, 2)	&	$1.326^{+0.068}_{-0.168}$	\\
\hline
\enddata
\end{deluxetable}
\end{center}
\label{eclipse_priors}

The best-fit model to the emission spectrum from the \texttt{ATMO} retrieval (Figure. \ref{fig:atmo_retrieval}) gives a $\chi^2_{\nu}$ of 1.24 (4 degrees of freedom) with a near-isothermal TP profile. The best-fit blackbody temperature using the PHOENIX stellar model \citep{husser_new_2013} grid (logg = 4.5 and logZ = 0) interpolated to T$_{eff}$ = 6390 K is 2283$\pm$64K with a $\chi^2_{\nu}$ of 1.23 (1 degree of freedom). Based on the formalism described in \citep{cowan_statistics_2011}, this dayside temperature would suggest a circulation efficiency $\varepsilon$ of 0.44 assuming zero Bond albedo and the predicted nightside temperature would be 1572K. $\varepsilon$=0 represents the no-circulation limit where the nightside temperature is 0K and $\varepsilon$=1 represents the full-circulation scenario where nightside temperature is the same as the dayside. The measured dayside temperature is consistent with the dayside-only heat redistribution scenario. The retrieved metallicity is consistent with the solar value (Figure. \ref{fig:atmo_corner}). While \cite{sheppard_hubble_2021} retrieved metallicity $\sim$1$\sigma$ higher than solar, \cite{lewis_into_2020} reported values consistent with solar metallicity. Considering the large uncertainties on metallicity from all three studies, the current datasets can not well constrain the models.

\cite{lewis_into_2020} found H- abundance several orders of magnitude larger than equilibrium chemistry is needed to best fit the transmission spectrum. We did not find evidence for significant H$^{-}$ opacity from the emission spectrum and we believe this could be due to: (1) The blackbody-like emission spectrum and near-isothermal TP profile make the abundances largely unconstrained. (2) The H- constrain in \cite{lewis_into_2020} comes from a combination of WFC3 G280 and G141 data. \cite{wakeford_into_2020} did not find evidence for H- solely based on the G280 optical dataset. Therefore any uncorrected offsets between G280 and G141 spectra can potentially lead to high H- abundance. (3) \cite{sheppard_hubble_2021} did not retrieve significant H- abundance with both G280 and G141 datasets using a different retrieval code. So it could also be due to different forward model assumptions.

Heavy metals such as Fe I, Fe II, and molecule like TiO are considered major optical absorbers that can induce temperature inversion in hot Jupiter atmospheres \citep{fortney_unified_2008, lothringer_uv_2020}, and they have been detected in the optical transit spectrum with HST/STIS G430L and G750L \citep{fu_hubble_2021-1}. There are no significant NUV heavy metals nor optical TiO absorption detected in both STIS \citep{sheppard_hubble_2021} and UVIS \citep{wakeford_into_2020, lewis_into_2020} transmission spectra for HAT-P-41b. The non-inverted TP profile retrieved in this HAT-P-41b emission spectrum is consistent with the non-detection of NUV/optical absorbers in the transit spectra.

\section{Compared with other hot Jupiter atmospheres}

To quantify the WFC3/G141 emission spectra from various hot Jupiters and compare HAT-P-41b to them, we took the 1.4$\mu m$ water feature strength index (S$_{H_2O}$) reported in \cite{mansfield_unique_2021} and plotted them in the dayside temperature versus gravity parameter space (Figure \ref{fig:compare}) first introduced in \cite{parmentier_thermal_2018}. The important physical transitions including TiO, H$^-$, and molecular dissociation are also the same as shown in \cite{parmentier_thermal_2018}. The larger negative S$_{H_2O}$ values (redder) represent stronger water emission features while the larger positive values (bluer) indicate more prominent water absorption features. The index is centered at zero for a blackbody-like featureless spectrum. In addition, we also scale the size of each point based on the host star temperature (T$_{eff}$) to reflect the effect of increased FUV/UV radiations from earlier type stars \citep{lothringer_influence_2019, fu_strong_2022}.

We have only seen strong water absorption features among cooler hot Jupiter (T$_{day}<$2200K) atmospheres \citep{kreidberg_precise_2014, line_no_2016} driven by their decreasing TP profiles. This is due to heavy metal absorbers condensing out of the atmosphere at low temperatures. On the other hand, when the atmosphere becomes too hot (T$_{day}>$2700K), water molecule dissociates and H$^-$ continuum opacity starts to dominate. As a result, we have only observed featureless emission spectra among the hottest planets. In between the two regions ($\sim$2200K to 2700K) we see a transitional parameter space where atmospheres can be inverted \citep{fu_strong_2022, evans_ultrahot_2017} or isothermal \citep{fu_hubble_2021, mansfield_unique_2021}. HAT-P-41b sits in the middle of this transitional space. Its mostly featureless WFC3/G141 emission spectrum is similar to that of WASP-76b, WASP-74b, WASP-19b and KELT-7b. However, it is different from WASP-121b and KELT-20b where we saw evidence for water emissions. Although previous models have indicated multiple important physical transitions \citep{parmentier_thermal_2018} happening at this temperature and gravity range, it is not yet well understood what exact physical processes drive the different thermal structures and emission spectra of these planets. In addition, host star type could be another important determining factor since strong FUV/UV radiations from the host star can strengthen atmospheric thermal inversion \citep{fu_strong_2022, lothringer_influence_2019} as they deposit significant energy into the upper layers of the planetary atmospheres. This is supported by the large emission features saw in KELT-20b which orbits an A-type star. The lack of FUV/UV flux from F-type host star of HAT-P-41b could be another cause for its isothermal atmosphere.

\section{Conclusion}

We present the most complete emission spectrum for inflated hot Jupiter HAT-P-41b. The spectrum is close to blackbody-like with no significant molecular absorption or emission features. The best-fit \texttt{ATMO} model shows an isothermal TP profile agreeing with the dayside heat redistribution scenario and a metallicity value consistent with the solar value. The non-inverted TP profile is consistent with the non-detection of NUV/optical absorbers in the transit spectra. Significant H$^-$ opacity suggested in \cite{lewis_into_2020} is not required in the model to adequately fit the emission spectrum. We also do not retrieve a metal-rich atmosphere as indicated in \citep{sheppard_hubble_2021}. However, the emission spectrum does not well constrain the atmosphere metallicity due to the limited wavelength coverage. The featureless emission spectrum of HAT-P-41b indicates planets with a dayside temperature around 2300K may have relatively isothermal TP profiles in the absence of heavy metal absorbers and strong host star FUV/UV radiations to drive thermal inversion. The comparison of HAT-P-41b to other similar hot Jupiters paints a murky picture of how atmospheric physical properties transition from cooler to hotter planets. The combined effect of surface gravity, thermal dissociation, H$^-$ opacity, heavy metal absorbers, and host star type is yet to be disentangled. Future observations of more similar planets and follow-up JWST infrared measurements of HAT-P-41b will be the key to solving the mystery of hot Jupiter atmospheres.

\startlongtable
\begin{center}
\begin{deluxetable*}{cccc}
\tablecaption{\textbf{HAT-P-41b eclipse spectrum}}
\tablehead{\colhead{Wavelength midpoint ($\mu$m)} & \colhead{Bin width ($\mu$m)} & \colhead{Occultation Depth (ppm)} & \colhead{Uncertainty (ppm)}}
\startdata
1.1137	&	0.0186	&	186	&	108	\\
1.1509	&	0.0186	&	378	&	94	\\
1.1881	&	0.0186	&	171	&	86	\\
1.2253	&	0.0186	&	526	&	102	\\
1.2625	&	0.0186	&	427	&	87	\\
1.2997	&	0.0186	&	346	&	86	\\
1.3369	&	0.0186	&	435	&	91	\\
1.3741	&	0.0186	&	587	&	89	\\
1.4113	&	0.0186	&	419	&	91	\\
1.4485	&	0.0186	&	499	&	95	\\
1.4857	&	0.0186	&	383	&	94	\\
1.5229	&	0.0186	&	446	&	99	\\
1.5601	&	0.0186	&	505	&	117	\\
1.5973	&	0.0186	&	643	&	105	\\
1.6345	&	0.0186	&	733	&	102	\\
3.6000	&	0.3800	&	1842	&	319	\\
4.5000	&	0.5600	&	2303	&	177	\\
\hline
\enddata
\label{eclipse_spectrum}
\end{deluxetable*}
\end{center}