\section{Introduction}

The hippocampal formation has been shown to be involved in spatial mapping \citep{OKeefe1971, Taube1990, Hafting2005}. Analogies to geographic maps or computer graphics suggest that spatial mapping should involve associating a Euclidean space of unique location representations, for example Cartesian coordinates, with a set of locations in the real world. In neuroscience, such a Euclidean space is often theorized to be provided by entorhinal grid cells \citep{Hafting2005}. A single grid cell module can represent Euclidean space ambiguously, while an array of grid cells can theoretically create unique location representations \citep{Fiete2008}. In models that simulate spatial mapping of the brain at a low level, the population activity of multiple grid cell modules is often used in this way \citep{Lewis2019, Whittington2020}, or they use a single grid cell module to complement traditional Cartesian coordinates \citep{Milford2008}. These models imply that the brain maps space in a way that is fundamentally similar to metrical geographic maps or computer graphics. On the other hand, in cognitive-inspired robotics \citep{Kuipers2000}, spatial mapping is less about creating a single coherent Cartesian map, and more about representing relationships between environment segments, with a patchwork of small Cartesian maps optionally learned on top of this system.

In this paper I argue that the multi-module grid cell code is actually not ideal as an analog of long-range Cartesian coordinates for spatial mapping. Instead, I suggest that the hippocampal formation approaches spatial mapping as graph-learning, quickly learning graphs which store relations between environment parts, and I suggest that grid cells provide additional spatial information that is not captured by these relations. This alternate approach to spatial mapping loses some generalization benefits that come naturally with Cartesian coordinates, but this approach also has more potential for flexibility. I argue that fast graph learning, with a learned repertoire of relations, is a more likely fundamental algorithm for the hippocampal formation than Cartesian-like mapping.

\section{Model}

In this section I present the theoretical and empirical motivations for graph-based spatial mapping, then I describe the general model at a high level. Finally, I present a possible low-level implementation built mainly from empirically observed neuron types.

\subsection{Problems with using only grid cells for spatial mapping}

There are theoretical and empirical reasons to doubt that the brain uses the multi-module grid cell code to create an analog of Cartesian coordinates.

Theoretically, the multi-module grid cell code is initially appealing if you approach mapping as a process of moving from location to location and associating sensory input with those locations \citep{Lewis2019, Whittington2020}. However, this approach to mapping is inefficient in that it requires visiting every location in the environment. A more practical and less expensive mapping algorithm would represent environments as arrangements of parts, using sensory input to quickly infer that arrangement of parts \textbf{(Figure \ref{fig:efficient-mapping})}. This requires the agent to be able to detect locations of these ``parts'' from a distance. Thus, it becomes important for grid cells to support the two computational operations of adding a displacement to a location and inferring the displacement between two locations. It is controversial whether grid cells support inferring or applying spatial displacements; some show that it is theoretically possible \citep{Bush2015}, while others argue that it is too difficult \citep{Fiete2008}, requiring an impractically large lookup table. Without such an operation, grid cells lose much of their desirability as a foundation for a map, especially given the fact that this also eliminates the ability to infer large navigation vectors. In previous work \citep{Hawkins2019}, my collaborators and I attempted to solve efficient mapping with the multi-module grid cell code by posing spatial mapping as composition of grid cell maps, associating one map with each environment ``part'' and quickly detecting the relationships of these maps. The two downsides of this approach were that it could only translate the maps relative to each other, not rotate or scale them, and that it had no notion of distance between the agent and the sensed ``part'', eliminating opportunities for incorporating this distance into sensory processing (for example, scaling the image on the agent's retina according to the distance).

\begin{figure}[]
\begin{center}
\includegraphics[width=0.9\textwidth]{figures/efficient-mapping.pdf}
\end{center}
\caption{\label{fig:efficient-mapping}\textbf{Efficient mapping.} An agent (a rodent) observes an environment and learns a map. One mapping strategy is to visit each location in the environment and record what you sense. This is inefficient, as indicated by the number of dots in the second panel. It is much more efficient to split the environment into parts or chunks and somehow store their relative locations. One strategy, common in computer graphics, is to store Cartesian coordinates for each part. Another strategy is to instead store a set of parts and a set of relations between them.}
\end{figure}

Empirically, there are multiple reasons to doubt that grid cell modules are working closely together to create unique location codes. They are anatomically separated \citep{Stensola2012}, and anatomical studies have not yet found any significant population of neurons receiving projections from more than one or two modules. In electrophysiological recordings, grid cell modules appear to scale and distort independent of one another \citep{Stensola2012}. The anatomical separation issue has been addressed by models by relying instead on horizontal connections in the hippocampus \citep{Whittington2020}, but the independent scaling and distortion would not be expected in a population location code. In particular, independent distortion of grids eliminates even the theoretical possibility of detecting spatial displacements between locations.

Thus, grid cells likely play a role in spatial mapping, but there are reasons to doubt that the brain uses multiple grid cell modules to create a Cartesian-like map.

\subsection{Using relations for spatial mapping}

Instead of associating parts in the environment with unique Cartesian-like locations, the hippocampal formation could represent environments as a set of relations between parts. In this view, an environment is not represented as ``\texttt{Part 1 at (x1, y1), Part 2 at (x2, y2), Part 3 at (x3, y3), ...}'', but instead as ``\texttt{Part 1 related to Part 2 with (relation1), Part 3 related to Part 2 with (relation2), ...}''. This representation is fundamentally a graph \textbf{(Figure \ref{fig:graphs})}. Processing input by building up graph-like data structures was an early idea in the history of AI \citep{Minsky1974} and is highly general. ``Scene graphs'' are an active area of research in computer vision \citep{Chang2021}.

\begin{figure}[]
\begin{center}
\includegraphics[width=0.9\textwidth]{figures/fast-learned-graphs.pdf}
\end{center}
\caption{\label{fig:graphs}\textbf{Fast-learned graphs as a general data structure.} The system stores in memory a set of nodes, associations with those nodes, a set of edges, and associations with those edges. This figure shows a graph-based map of the environment from Figure \ref{fig:efficient-mapping}. A neural system could build graph-based spatial maps by attending to different parts, assigning nodes to each part, and detecting relations between attended parts and associating those relations with edges. After a graph has been learned, it can be used for prediction (via the down arrows).}
\end{figure}

During spatial mapping, these graphs could be quickly built up from a few glances. \textbf{Figure \ref{fig:graphs}} shows a high-level mechanism. The agent uses its sensory input to infer a description of the sensed environment part as well as the spatial relationship between the agent and the part. A graph node is allocated to this part, and it is associated with the part's description. The agent attends to another environment part, activating a different description and spatial relationship, and allocating a second graph node. From this pair of agent-part spatial relationships, the agent infers the spatial relationship between the two parts, and this relation is associated with this ordered pair of graph nodes. Now that this graph has been learned, it can be used to guide attention shifts and to predict what will be sensed after each attention shift. The size of the map scales with the number of environment parts. In environments with many parts, the number of parts can be reduced by chunking them into aggregate parts.

The upside of this approach is that it is flexible and can easily apply to both spatial and non-spatial domains. The set of relations between nodes can be fully learned; in spatial domains, they could simply be spatial displacement vectors, while in other domains they could be more similar to the relations that are expressed in language.

The downside of this approach, relative to the Cartesian map, is that it does not provide easy direct access to relations of arbitrary pairs of parts. When parts are assigned Cartesian locations, the displacement between any two parts can be trivially computed, but if instead a graph of relations between parts is learned, detecting the relation between a pair of parts requires the relation to have been previously learned, or it must be inferred by finding a path between the parts on the graph and aggregating the relations. The ``Tolman-Eichenbaum Machine'' \citep{Whittington2020} aims to get the best of both of these worlds, applying to non-spatial domains while working similarly to Cartesian coordinates, essentially associating information with individual graph nodes. In this paper, I instead propose that the hippocampal formation associates information with directed edges between nodes, and this admittedly makes such generalization require more advanced computation over graphs. This is further discussed in section \ref{section:graph-processing}.

In the domain of spatial mapping, the set of relations will tend to leave some information ambiguous. When the relations are simple displacement vectors, those displacement vectors will tend to be coarse for two reasons. First, these graphs should be built up very quickly, and given the uncertainty inherent to sensory input this is much easier if the displacements are coarse. Second, using neurons to represent and decode fine-grained displacement vectors would require a lot of neural tissue, but with coarse displacements a small number of neurons and synapses could represent displacements over medium-to-large distances. The ambiguity of these relations is where grid cells will often become useful. Grid cells themselves are ambiguous, but in a different way. They can be viewed as filling in the information that was is left ambiguous by the relation graph. This principle is elaborated in in section \ref{section:combined-building-blocks}.

Thus in this theory, the hippocampus quickly learns graphs, associating information with nodes and edges of the graph. Relations are associated with edges. Grid cell locations, among other things, are associated with nodes providing extra information that was not provided by the relation. Alternately, we can invert this narrative. Grid cells solve one part of representing location, but they do it ambiguously and in a way that does not provide good support for representing displacements, and those aspects are supported by fast-learned relation graphs.

\subsection{Detailed example low-level model}

Here I describe a neural implementation of a system that learns environments as relation graphs, mainly using cell types empirically observed in the hippocampal formation. This section walks through the building blocks of this model then describes how they are combined. The model is human-interpretable; once the model is understood, the reader will be able to mentally simulate it.

Multiple design decisions of this model were informed by empirical evidence. For example, how should the spatial relationship between the self and the sensed part be represented? Another example: should the conversion from egocentric to allocentric locations occur before detecting relations between parts, or after? From a theoretical standpoint these questions do not have single obviously-correct answers. Based on the empirical existence of object-vector cells (described below), this model uses object-vector cells and performs the coordinate transformation before detecting relations. The actual cortex may use a hybrid approach, and the approach may vary by species.

\subsubsection{A short review of existing neural building blocks}

This model is built up from a set of components that have each been empirically observed and/or used previously in other models. These components are reviewed below.

\textbf{Head-direction cells:} These cells fire when an animal is facing a particular direction relative to the environment \citep{Taube1990}. They are analogous to the compass rose of a map. Head direction cells are mapped onto environment directions using sensory heuristics, for example detecting the main axis of the environment and anchoring relative to that axis \citep{Julian2018b}.

\textbf{Object-vector cells:} These cells fire when an animal is at a particular distance and environment-centric direction from any object \citep{Hoydal2019}. They can be thought of as a reusable set of location representations that can be bound to any object, but they use head direction cells as their compass rose rather than anchoring to the orientation of the object. Object-vector cells and grid cells both rely on the brain performing a coordinate transformation based on head direction cells. An animal's sensory and movement information are naturally encoded in a viewer-centric or egocentric reference frame, and head-direction cells are needed to convert these sensory vectors and movement vectors into an environment-centric or allocentric reference frame \citep{Bicanski2018}. This coordinate transformation in this model largely mirrors the model of \cite{Bicanski2018}. Head direction cells provide a routing signal, specifying how vectors must be rotated to convert between viewer-centric and environment-centric directions. Object-vector cells are just one of many types of ``-vector'' cells observed in the hippocampal formation \citep{Bicanski2020}. For simplicity, this model only focuses on object-vector cells, but many other vector cells could be substituted.

\textbf{Grid cells:} These cells fire when an animal is at particular locations in the environment. The set of firing locations forms a triangular or hexagonal ``grid''. The firing patterns of different grid cells in a local area of neural tissue are similar, shifted / translated versions of the same grid. A cluster of grid cells with a complete set of translations of the same grid are known as a ``grid cell module''. This model uses a single grid cell module \citep{Stensola2012}. Multiple computational models of grid cells exist that show how they might perform path integration or dead reckoning \citep{Giocomo2011}. This model does not depend on which of these implementations grid cells are using, it only depends on grid cells being able to update given movement information.

\textbf{Grid cells that represent location of object in environment:} This model uses a second population of grid-like cells that tracks locations of objects rather than self-location. Following basic intuition of mapping systems, it is theoretically likely that objects are associated with their locations on the grid, but it is an open question whether any grid-like cells would need to be dedicated to representing this location. A set of experiments has shown a population of grid cells that track viewed location \citep{Killian2012, Nau2018, Julian2018a} and covertly attended locations \citep{Wilming2018}. One possible interpretation of these empirical results is that the recorded grid cells were in fact these theorized grid cells. In this model, this location acts as a routing signal between object-vector cells and grid cells, and it is computed by detecting which object-vector cells fire with which grid cells.

\textbf{Displacement cells:} Based on theoretical considerations, it has been proposed that there are likely neurons that detect the spatial displacement between two metric location representations \citep{Bush2015, Hawkins2019}. This model uses the same mechanism as previously proposed, with a population of displacement cells responding to or invoking sudden changes in a metric location representation. One noteworthy difference in this model is that the displacement cells detect sudden changes in object-vector cells, rather than sudden changes in grid cells.

\textbf{Memory nodes:} Index theory of hippocampus suggests that some subset of hippocampal neurons form representations whose job is simply to learn associations with multiple other populations \citep{Teyler1986, Josselyn2020}. The representations themselves act as an index (similar to an index of an array) for recalling all attributes of a memory. Neurons in these representations, sometimes referred to as engram cells, are often hypothesized to carry no meaning of their own. Rather, they are given meaning via learned associations to populations that carry the content of the memory. For example, in this model these cells represent instances of environment parts. If an agent sees an environment part in Environment 1, it will create a memory node for this instance of this environment part, and it associates with it a description of the part, and it also associates its location and its displacement relative to other nearby parts. This content is all stored by learning associations with memory node representations. If a similar environment part is encountered in Environment 2, a different memory node will be allocated to that part and its attributes.

\subsubsection{How this model combines these building blocks}
\label{section:combined-building-blocks}

This model can be understood by first considering what it can do without grid cells, then adding grid cells \textbf{(Figure \ref{fig:circuit})}.

\begin{figure}[]
\begin{center}
\includegraphics[width=0.9\textwidth]{figures/circuit.pdf}
\end{center}
\caption{\label{fig:circuit}\textbf{Neural mechanism for quickly building memory graphs.} Spatial mapping can work without grid cells. Adding grid cells makes the spatial map more precise and enables path integration (also known as ``dead reckoning''). \textit{(Left)} A minimal circuit. Information from sensory input goes through a coordinate transformation, resulting in object-vector cells and an allocentric description of the sensed part. Spatial relations between object-vector cells are simply coarse displacement vectors. Graph nodes are represented by random sparse hippocampal activity patterns. Graph nodes are associated with other temporally adjacent graph nodes, enabling the system to either oscillate between nodes or allow the attention system to know which nodes can be attended to. After relation information has been learned, object-vector cells can be coarsely predicted during attention shifts, but they cannot be precisely predicted. \textit{(Right)} The same model, now with grid cells added. Movement information goes through a similar coordinate transformation and is used to update grid cells. Using grid cells and object-vector cells, the location of the sensed part on the grid is detected and associated with the graph node. The coarse relation between parts and the grid location of each part can together be used to accurately predict object-vector cells during movement and during attention shifts.}
\end{figure}

In this model, self-location is primarily represented by object-vector cells in conjunction with a graph node. The network represents one environment part at a time, quickly switching between representing different environment parts. This switching can be interpreted as attention shifts or as fast oscillations between nodes, enabled by learned associations between temporally adjacent nodes. Each time a node becomes active, it drives the part description to become active. Additionally, each ordered pair of nodes drives a learned relation between parts to become active. During attention shifts, these learned relations can be used to coarsely predict which object-vector cell should become active.

Now we add grid cells. Similar to before, self-location is primarily represented coarsely by object-vector cells and a graph node, but grid cells provide a complementary uniform resolution location, augmenting the agent's location representation with extra precision. In this extended model, the grid location of each environment part is associated with graph nodes. This information enables the model to more precisely predict object-vector cells during attention shifts. Grid cells receive movement information and perform path integration, and these updates are propagated to the object-vector cells.

An important design principle of this model is that it uses different spatial resolutions in grid-like vs.\@ vector-like populations \textbf{(Figure \ref{fig:ambiguity})}. In this model, grid cells use a fixed spatial resolution so that they can perform path integration, whereas object-vector cells use varying resolution (which is speculative -- not empirically verified) so that they can be distributed over a larger area using a relatively small number of cells. This division of labor enables the object-vector cells to be path-integrated by the grid cells, and it allows precise coding of locations that are far from environment parts via the combined representations of the coarse-but-unambiguous object-vector cells and the precise-but-ambiguous grid cells. Similarly, displacement vectors between objects are stored using varying resolution while the grid locations of objects are stored in fixed resolution, and so the two combined learn the precise relative locations of objects. Considering the fundamentals of sensory processing, this design is a natural fit in multiple ways. First, as a viewer moves away from an object, the image of the object on the retina changes more slowly, so it is natural for sensory-driven cells like object-vector cells to have larger and larger fields as distance increases. Second, displacement vectors between objects can be easily identified coarsely, but often not precisely, so it is pragmatic for the system to have one set of cells that quickly learn the coarse displacements between objects, while another set of cells learn the precise relative locations more slowly by learning the objects' grid locations, similar to how robots iteratively update locations of objects when performing Simultaneous Localization and Mapping (SLAM).

\begin{figure}[]
\begin{center}
\includegraphics[width=0.9\textwidth]{figures/ambiguity.pdf}
\end{center}
\caption{\label{fig:ambiguity}\textbf{Grid cells' ambiguity complements the ambiguity of relations.} \textit{(Top-left)} A 2x2 taxonomy of this model's representations of location. Vector-like cells use varying resolution, making them long-range but coarse. Grid-like representations are ambiguous over long ranges, but they use uniform resolution. \textit{(Top-right)} Depiction of an example ambiguous scenario. The blue circle indicates how the object-vector cells are mapped onto space, while the black rhombus indicates the grid cells, which tile the space. The \textit{Location of part on grid} describes where the center of the blue circle lands on the black rhombus, giving us this picture. As shown here, the grid cells alone do not uniquely determine which object-vector cells should be active. Different sets of object-vector cells that are consistent with the active grid cells are shown as dotted circles. \textit{(Bottom-left)} This problem is solved by using the memory graph's coarse displacement information. During path integration, the object-vector cells give preference to cells that are near the previously active cells. When switching between nodes on the memory graph, the orange coarse displacement information is used to determine which subset of the object-vector cells should be given preference. Thus, the \textit{Location of part on grid} and \textit{Coarse displacement} work together to enable mapping spaces much larger than one grid cell period and potentially much larger than the object-vector cell range. \textit{(Bottom-right)} This operation of predicting object-vector cells requires all of the other populations, and the computation also depends on the current active object-vector cells, as indicated by the recurrent arrow.}
\end{figure}

A core capability of this model is its ability to update object-vector cells during movement and during attention shifts. These operations can be performed precisely when the model uses grid cells. This relies on being able to route between grid cells and object-vector cells given the grid location of an object, shown as a three-way arrow in \textbf{Figure \ref{fig:circuit}}. This is built on the intuitive property that if an agent knows its location on the grid and it also knows an object's distance and direction from the agent, it can infer the object's location on the grid. The required connectivity for this three-way circuit is similar to that of cells that detect and apply displacements to grid cells \citep{Bush2015, Hawkins2019} or rotations to vector cells \citep{Bicanski2018}. Each cell receives projections from pairs of cells in the other populations so that they receive strong excitatory input when that pair of cells becomes active. However, in this model there as an additional complication, due to the ambiguity of grid cells \textbf{(Figure \ref{fig:ambiguity})}. Because this model's object-vector cells can cover a much larger area than one grid cell spatial period, it is not possible for this three-way circuit to unambiguously infer object-vector cells from grid cells. This problem is solved by making the object-vector cells' activations competitive and using learned coarse displacements to supply them with additional information. During path integration, the grid cells update in response to movement, and this update is propagated to the object-vector cells. Because no node switches occur, no displacement cells are activated, so object-vector cells that represent locations near that of previously active cells are given preference, solving the ambiguity problem during path integration. This leaves the problem of ambiguity when rapidly switching between objects. This rapid switch causes displacement cells to activate, providing excitatory input to a certain subset of the object-vector cells, depending on which ones are currently active. This excitatory input combined with the excitatory input from the three-way circuit will activate the appropriate object-vector cell. This mechanism demonstrates how object-vector cells can be updated by grid cells, and it demonstrates how object-vector cells can rapidly switch between representing different objects if the system has learned those objects' grid locations and the coarse displacements between those objects.

Putting these pieces together, we get a human-interpretable mapping system that readers will be able to mentally simulate. Consider the scenario of an animal arriving in a novel environment. The agent already has a grid cell and head direction cell representation active (or, if not, it can randomly activate one). The agent sees an object, inferring some description of the object and the agent's location relative to that object, and it uses head-direction cells to perform a coordinate transformation, inferring allocentric object-vector cells. Using the current active grid cells and object-vector cells, the agent infers the object's location on the grid. Because the location of the object on the grid is now known, updates to grid cells will drive updates to object-vector cells without the agent needing to sense the object again, enabling prediction during self-motion. We continue our mental simulation, now focusing on how the map is learned. The agent forms a new memory, allocating a memory node and associating it with this object's description and its location on the grid. The agent now senses a second object in the environment, going through the same above process. This drives a different object-vector cell to suddenly become active, driven by this sensory input. The displacement between the vectors (i.e.\@ between this object-vector cell and the previously active object-vector cell) is detected by a set of displacement cells, and this displacement is associated with this pair of memory nodes. Finally, the agent attends to the first object. The map ought to enable the agent to anticipate which object-vector cells should become active. Using the combination of the currently active object-vector cells, the learned coarse displacement information, the currently active grid cells, and the learned location of the object on the grid, the correct newly active object-vector cells can be anticipated. This can feed backward into the sensory system, enabling prediction. Thus, as the agent attends to different objects, grid cells and head direction cells remain stable while all the other populations switch between objects. Continuing this learning process, the system can learn arbitrarily large maps as long as the environment-parts are not long distances apart (i.e. longer than the object-vector cell range).

\subsection{Graph processing}
\label{section:graph-processing}

By moving from a system built on Cartesian coordinates analogs to a system built on graphs, we gain generality, but we also lose some conveniences of Cartesian coordinates. With Cartesian coordinates, after each environment part is attended to once, the system can trivially compute displacements between any pair of parts. In a graph-based system, computing a displacement between two parts either requires having previously observed those parts in that order, or it requires finding a path between them in the graph and aggregating over that path. We even see this issue coming up in the mental simulation from the previous paragraph. In that example, when the agent attended back to the first object, it had never previously experienced the ordered transition from Object 2 to Object 1, it had only experienced Object 1 to Object 2. It is trivial on a computer to infer this reversed displacement, but can the brain's neural circuitry do it?

It an open interesting question what types of computations the brain performs on these graphs. There are three broad possibilities. In one extreme, the system may simply require all relevant transitions to be learned by brute force. If, in practice, transitions tend to happen in certain orders, this may be sufficient. In the other extreme, the system may be able to learn a minimal graph and perform advanced computations on it during inference. This may involve a form of replay \citep{Foster2017}, with different readout neurons performing different computations on the replayed sequence. Finally, in a middle-ground possibility, the brain could quickly learn a minimal graph, then later during rest or sleep fill in missing edges, possibly relying on replay.

\section{Predictions of the model}

In this paper I have argued that the hippocampal formation learns environments as arrangements of familiar parts, and in particular that these arrangements of parts are learned as relation graphs. One high-level implication of this idea is the hippocampal formation does not need to form an association at each location in an environment. In the model of \cite{Whittington2020}, every place cell population code is learned as an attractor state, so every place cell representation that occurs in the environment is memorized. In this model, only the graph of environment parts is memorized, so learning the environment does not involve learning these place cell attractor states. Depending on how ``place cells'' are defined, it is possible that the memory index / engram cells representing nodes in this graph will be classified as place cells. These cells do not in fact encode self-location, rather they become active when the agent is attending to a particular environment part, though this may correlate strongly with certain locations.

Diving deeper into the particulars of this model, further predictions can be made. This model suggests that grid cells tracking eye position \citep{Killian2012, Nau2018, Julian2018a} and covertly attended locations \citep{Wilming2018} are actually tracking locations of attended objects, or the locations of the subject of any other type of ``vector'' cells \citep{Bicanski2020}, e.g.\@ goals. Object-vector cells have been seen to respond at multiple locations in environments that have multiple objects \citep{Hoydal2019}, and this model suggests that object-vector cells will either respond to different objects at different times in the theta band, quickly oscillating between multiple objects, or will shift with attention -- either way, there will be time-based separation. It is important to note that this paper focused on interpretable models and toy datasets. Actual neural tuning curves will be more complicated.

\section{Discussion}

The hippocampal formation is often described as a spatial mapping system. It is worth asking whether spatial mapping is the base algorithm of the system, or if there is some other base algorithm and spatial mapping is just one of its applications. With this model, I aimed to frame spatial mapping as just one instantiation of a more general fast-learned-relation-graph algorithm. By allowing the repertoire of relations to be learned, this system becomes highly general, applicable to spatial and non-spatial problem domains. In this paper I did the preliminary work of framing spatial mapping as a graph-learning process. Future work will exploit this gained generality.

This point of view is compatible with the high-level point of view of TEM \citep{Whittington2020}, but it differs in some important details. Using the language of this paper, TEM attempts to capture all relations by associating location information with graph nodes, not by associating information with graph edges. I argue that, though node-only associations are desirable when possible, it will not always be possible, and the ability to associate arbitrary relation information with graph edges is important. It is likely that a slight variation on TEM could learn exactly the model that I have presented here.

\section*{Acknowledgements}

Many thanks to Mirko Klukas, Jeff Hawkins, and Subutai Ahmad for helpful discussions.

\bibliographystyle{apa}


\section{Introduction}

Intensity-modulated radiation therapy (IMRT) has become one of the predominant forms of radiation therapy (RT) currently utilized in clinical treatments. This prevalence is the result of ongoing research efforts that continue to improve the conformity and accuracy with which IMRT plans can be delivered. As the precision of equipment and accuracy of treatment planning algorithms improves, understanding the radiobiology involved becomes increasingly important \cite{rosenthal_mouse_2007,liu_lessons_2013,butterworth_small_2015,coleman_improving_2016}. As such, preclinical studies are essential for fully investigating novel RT treatment techniques and technologies before they can be implemented in the clinic \cite{butterworth_small_2015}. Accurate delivery of radiation to the site of disease in humans is a critical part of achieving a favorable outcome, so it would reasonably follow that the analogously accurate delivery of radiation to small animals desireable in order to facilitate translation of preclinical findings to clinical improvements. Given the small scales involved in mice, the tolerance of spatial deviation is on the submillimeter scale, rather than the multiple millimeter tolerance acceptable in humans \cite{desrosiers_use_2003}. For these reasons, the need for IMRT and image guided RT in preclinical studies at this submillimeter scale is important for enhancing the clinical relevance of associated results.

Currently, small animal IMRT methods include the use of millimeter scale secondary collimators with standard geometric shapes, and 3D-printed compensators to create a specific fluence map.
Secondary collimators has been implemented in two forms, the first of which utilizes a series of motorized jaw positions to create small rectangular segments delivered in a serial fashion and the second uses a Sparse Orthogonal Collimator (SOC), which consists of 4 sets of orthogonally positioned leaves \cite{stewart2013two,nguyen2016novel,reinhart_uimrt_2018,woods2019sparse1,woods2019sparse2}. Aperture-based methods are limited by a minimum size of the apertures, as well as potentially long treatment times requiring numerous exposures through different jaw/leaf positions. Another approch is through compensator-based methods which utilize 3D printed patterns through which the beam is delivered at specified angles. In order to achieve their attenuation, a shaped polylactic acid (PLA) compensator can either be packed with sodium iodide powder, or the compensator itself can be printed out of PLA impregnated with copper or tungsten \cite{yoon2020method,gage2020small}. Some advantages of the compensator method are the submillimeter resolution achievable by modern 3D printers, a reduction in the number of segments that need to be delivered in order to achieve a conformal dose, and the absence of additional mechanical parts which require maintenance and quality control \cite{stewart2013two,gage2020small}. However, 3D-printed compensators are specific to each plan and beam angle, thus a new set needs to be made for each plan, consuming additional material and printing time. The printing process, while capable of the aforementioned submillimeter accuracy, is more prone to error in abutting regions of greatly differing attenuation (and therefore thickness). Furthermore, adjacent regions with a large difference in intensity will lose lateral electronic equilibrium, as electrons scatter out of the copper in the highly attenuated area into the column of air above low-attenuation segments.

Total variation regularization (TVR) is most often used in digital image processing for noise removal. Total variation refers to the integral of the absolute gradient of the signal. TVR is based on the principle that noisy signals have high total variation \cite{rudin1992nonlinear}.TVR has been adopted in the field of radiation oncology in the form of compressed sensing used during inverse treatment planning\cite{zhu2009search}. Compressed sensing facilitates planning by reducing the complexity of the optimized fluence map, reducing the number of multi-leaf collimator segments for a given treatment plan without sacrificing plan quality. A TVR term can be added to the cost function in order to minimize the complexity of fluence patterns for IMRT segments, or in this case, modulation between adjacent beamlets.

In this work we develop a TVR-based inverse planning method for IMRT in small animals based on a novel micrometer-scale dose deposition kernel designed for accurate treatment planning at the sub-millimeter scale. This approach is designed to reduce the complexity of the modulation pattern, overall beam-on time, and the amount of material used to create 3D printed compensators, while maintaining or improving the DVH that can been achieved with previous compensator-based IMRT methods.

\section{Methods and materials}

\subsection{TVR-based inverse planning}

In TVR based inverse planning for small animal, a gradient operator is defined as
\be\label{equ_xuv}
\nabla x(u,v) = \bigl|x_{u,v} - x_{u-1,v}\bigl| \;+  \; \bigl|x_{u,v} - x_{u,v-1}\bigl|,
\ee
where $x$ is beamlet intensity, and $u$ and $v$ are its row and column indices, respectively.
The total variation (TV) of intensity is defined as the
summation of the gradients in beam view of all fields as the following  $\ell_1$ norm term,
\be\label{intensity_TV}
\phi=\sum_{u,v}  \nabla x(u,v)=\|Gx\|_1.
\ee
The complexity of the fluence map can be evaluated by TV of intensity.
The TVR based IMRT includes this extra gradient term, and so it is given by
\be\label{form_compressed}
\begin{array}{rl}
\mbox{minimize} & h(y) + \beta \|Gx\|_1 \\
\mbox{subject to}
 & y=Dx,\;\; 0\preceq x \preceq u_x,
\end{array}
\ee
where
$D$ is the dose matrix,
$y$ is the dose to the voxels,
and $u_x$ is the upper bounds of beamlet intensity,
$\beta$ is the weight of variation.
The function $h(y)$ is the cost function for conventional IMRT, and it can be a linear, quadratic or other nonlinear functions of the voxel doses of the targets and organs at risk (OAR).
Due to the  $\ell_1$ norm term, the TVR-IMRT cost function in (\ref{form_compressed}) is now convex but not differentiable.
To apply the efficient
Limited-memory Broyden–Fletcher–Goldfarb–Shanno optimization algorithm (LBFGS), in which the derivative of the objective function is used, it is therefore necessary to approximate $|x|$ by a smooth function. The function $|x|$ can be approximated by $\sqrt{x^2+\nu^2}$ as in \cite{ramirez2014x2}, or by parabolic functions at the neigherhood of $x=0$. The differentials involved in optimization can be derived accordingly.

\subsection{Development and 3D Printing of Compensators}

A photon beam passing through a sample of
thickness $z$ is attenuated via numerous types of interactions with the surrounding medium. The photons are transmitted according to
Beer–Lambert’s law \cite{kucuk2012determining}.
The parameters of polychromatic transmission through both copper- and tungsten-PLA filaments were calculated by measuring the radiation transmitted through material thicknesses ranging from 0 to 20 mm with an ionization chamber. See Figure~\ref{SARRP-transmission} in the next section for the copper-PLA transmission. The radiation attenuation of tungsten-PLA is about four times larger than copper-PLA. The copper-PLA was used to shape the fluence map inside the beams-eye-view region containing the target, while the tungsten-PLA was used to print the trimmers to block the radiation outside the desired field (analogous to portions of beam collimated by MLCs and jaws, respectively, with cinical treatment machines).
In the case of our copper filament compensators, a thin base layer of 0.4mm is first printed in the modulation area, and the maximum thickness of the compensators was assumed to be 20mm.
Therefore, the bound of beamlet intensity $x$ in (\ref{form_compressed}) should be replaced by $l_x\preceq x \preceq u_x,$, and the ratio of minimum intensity $l_x$ to maximum intensity $u_x$ for the optimization was set to
\be\label{l_x_u_x}
l_x/u_x=T_{\tiny\mbox{Cu}}(20)/T_{\tiny\mbox{Cu}}(0.4),
\ee
where $T_{\tiny\mbox{Cu}}(z)$ is the transmission through Copper-PLA filament of thickness $z$.
The optimization was first run without constraint on fluence intensity, then $u_x$ was set to an appropriate value accordingly.
Once obtaining the optimized bixel intensity $x$, the compensator thickness can be computed accordingly. Denote $x=\mbox{col}(x_1,\cdots,x_m)$, where $m$ is number of gantry angles. For the $i$-th each gantry angle,
the expose time $t_i=\mbox{max}(x_i/T_{\tiny\mbox{Cu}}(0.4))$. The dose in a single time interval is $x_i/t_i$, and thus the thickness $z_i$ is given by
\[
z_i=T_{\tiny\mbox{Cu}}^{-1}(x_i/t_i),
\]
The total variation (TV) of thickness is defined in the similar way as in (\ref{equ_xuv}) and (\ref{intensity_TV}),
\be\label{thickness_TV}
\varphi=\sum_{u,v}  \nabla z(u,v)=\|Gz\|_1,
\ee
which can be used to evaluate the complexity of 3D printing surface inside beamview.

The compensators used in this study were printed using a Raise3D Pro2 Plus (Raise3D Technologies Inc, Irvine, CA). The current version of these 3D-printed compensators have a tab and outer ring for mounting and alignment. The associated tungsten trimmer is positioned downstream, in contact with the compensator, in order to remove any scatter from outside the desired region. See Figure~\ref{SARRP-transmission}.

\subsection{IMRT plan for a mouse heart with lung OARs}

A comprehensive practical implementation of IMRT using 3D-printed beam compensators on the Xstrahl Small Animal Radiation Research Platform (Xstrahl Inc, Suwanee, GA) (SARRP) was used to illustrate the TVR method.
The Source-Axis-Distance (SAD) of the SARRP is 350mm, and, for the customized collimator, the source-to-compensator distance is 290mm. The bixel width was set to 1.0mm at SAD for all treatment planning, thus the bixel width of the 3D-printed compensators was 290/350=0.8286mm.

Plans were generated from a CT taken for use in a study on heart toxicity, using the built-in imaging system on the SARRP \cite{dreyfuss_novel_2021}. Hounsfield unit (HU) values were discretized to 5 categories: air, lung, fat, muscle, and bone, which were assigned an electron density by Muriplan (xStrahl) which are based on ICRU-44\cite{white_report_1989}. The CT dimensions used for optimization were $126\times 126\times 134$ voxels with a voxel width of 0.254mm. The lengths of target in $xyz$ direction are 8.9, 5.6, and 8.6mm. Relative dosimetry measurements and subsequent dose calculation using the 3D printed compensators was verified as in the previous work \cite{redler2020small}.

The heart was designated as the planning target volume (PTV) and the lungs were the only organ at risk (OAR) contoured and given dose constraints.
The plan generated by each optimization method was evaluated by comparing the $D95_{\tiny\mbox{PTV}}$[\%], $V95_{\tiny\mbox{PTV}}$[\%], $D95_{\tiny\mbox{lung}}$[\%], $V95_{\tiny\mbox{lung}}$[\%], and $V110_{\tiny\mbox{lung}}$[\%], as well as the total variation (TV) of both the fluence and compensator thicknesses, and total beam-on time.

In order to evaluate the accuracy of the delivered plan, each of the five beams, as well as the composite dose was delivered to EBT2 film placed in a $32\times 32$ mm high-density polyethylene (HDPE) phantom with the center of the film placed at isocenter. The accuracy of the dose calculation algorithm and resultant compensator was assessed both on a per-field and composite dose planar gamma analysis by using previously published software \cite{low1998technique} \cite{mark2015gamma}.

\section{Results}

\subsection{Experimental setup and copper filament transmission}
\begin{figure}[ht]
\centerline{
		\includegraphics[height=73mm]{fig/SARRP}\hspace{10mm}
\begin{minipage}[b]{0.3\textwidth}
	\centerline{
		\includegraphics[width=65mm]{fig/Cu_PLA_transmission}
		}\vspace{3mm}
	\centerline{
		\includegraphics[width=54mm]{fig/compensators_trimmer_fat2}
		}
\end{minipage}
}
	\caption{Left panel: Irradiation and measurement setup. Right upper panel: Radiation transmission of copper-PLA. Right lower panel: 3D-printed tungsten-PLA trimmer (black) and corresponding copper-PLA IMRT and TVR-IMRT compensators.}\
	\label{SARRP-transmission}
\end{figure}

For this study, a custom brass collimator and compensator-holder was machined to fit the Xstrahl SARRP. This collimator held a series of notches and alignment pins to ensure that the compensator for each beam was held at the same source to compensator distance, without allowing rotation around the beam axis. The $32\times 32$ mm HDPE phantom can be seen with the film inside placed at isocenter (Figure 1, left panel).

Transmission through a range of Cu-PLA thicknesses was measured in order to translate desired bixel intensity into a physical compensator (Figure 1, right upper panel). Cu-PLA compensators for each beam angle were designed with and without TVR. A tungsten-PLA trimmer was designed to reduce out of field scatter and aligned to the compensators with three tongue-and-groove points (Figure 1, right lower panel).

\subsection{Plan delivery statistics and compensator dimensions}
\begin{figure}[ht]
	\centerline{
		\includegraphics[width=50mm]{fig/hp_newG_TV0_CTDose}\hspace{-9.7mm}
		\includegraphics[width=50mm]{fig/hp_newG_TV1_CTDose}\hspace{6mm}
		\includegraphics[width=60mm]{fig/hp_newG_all_dvh}
	}
	\caption{Dose distribution and DVH comparison for IMRT (solid) and TVR-IMRT (dashed)
	}
	\label{fig_Dose}
\end{figure}

The dose heat maps shown on an axial slice of the mouse CT (Figure~\ref{fig_Dose}) containing the heart and lungs and the corresponding DVH were similar between both planning methods,
with the TVR-IMRT having slightly more homogenity (lower hotspot). This is also shown by the data in Table~\ref{table_D95V95}.

\begin{table}[htbp]
\caption{DVH statistics comparison of optimization results}
\centerline{\small
\begin{tabular}{l|rrr|rrr}
  \hline\hline
&  PTV $D_{95}$ [\%] &  $V_{95}$ [\%]   &$V_{110}$ [\%]   & Lung $D_{95}$ [\%] &  $V_{50}$ [\%] & $V_{95}$ [\%] \cr
\hline
      IMRT &     94.5\%&     94.2\%&  0.9\%&   0.0\%&     21.9\%&      0.1\% \cr
TVR-IMRT&     94.4\%&     93.9\%&  0.0\%&   0.1\%&     23.1\%&      0.1\% \cr
\hline\hline
\end{tabular}
}
\label{table_D95V95}
\end{table}

All dose objectives were met by both plans and were within 0.5\% for all but the PTV $V_{110}$ and the Lung $V_{50}$ (Table 1). The TVR-IMRT plan did not deliver 110\% of the dose to any portion of the PTV, while the IMRT plan delivered it to 0.9\%. The Lung $V_{50}$ was 1.2\% lower for the IMRT plan than the TVR-IMRT plan (Table~\ref{table_D95V95}).

\begin{figure}[ht]
	\centerline{\includegraphics[width=170mm]{fig/hp_newG_TV0_fluence}}
	\vspace{2mm}
	\centerline{\includegraphics[width=170mm]{fig/hp_newG_TV1_fluence}}
	\caption{Bixel intensity comparison of the compensators designed for each gantry angle. Top: IMRT. Bottom: TVR-IMRT}
	\label{fig_beamlet_intensity}
\end{figure}

\begin{figure}[ht]
	\centerline{\includegraphics[width=170mm]{fig/hp_newG_TV0_thick}}
	\vspace{2mm}
	\centerline{\includegraphics[width=170mm]{fig/hp_newG_TV1_thick}}
	\caption{Compensator thicknesses for each gantry angle.  Top: IMRT. Bottom: TVR-IMRT}
	\label{fig_compensators_thick}
\end{figure}

\begin{table}[htbp]
	\caption{IMRT vs TVR-IMRT}
	{\small
		\begin{tabular}{l|r|r|r|r}
			\hline\hline
			&  Intensity TV & Thickness TV (mm) & Total thickness (mm) & Exposure time (s) \cr
			\hline
			IMRT &      997 (   0\%)&     1620 (   0\%)&     1771 (   0\%)&   280 (   0\%) \cr
			TVR-IMRT&      517 ( -48\%)&     728 ( -54\%)&     902 ( -28\%)&   218 ( -22\%) \cr
			\hline\hline
		\end{tabular}
	}
	\label{table_compare_TV}
\end{table}

The TV of the delivered dose and physical thickness for the TVR-IMRT compensator was reduced by 48\% and 54\% respectively, when compared to the IMRT compensator as shown in Table~\ref{table_compare_TV}. TVR reduced the sum of the thicknesses of each discrete region in the beam path (corresponding to the filament used during printing) by 28\%. The time to deliver the plan was similarly decreased by 22\% for TVR-IMRT.

\subsection{Gamma analysis}

The planar dose profile (perpendicular to beam-axis) of each individual compensator for each of 5 gantry angles (0$^\circ$, 72$^\circ$, 144$^\circ$, 216$^\circ$, 288$^\circ$),
were all delivered from gantry angle 0$^\circ$ to a HDPE phantom for ease of comparison and to facilitate gamma analysis for both IMRT and TVR-IMRT (Figure~\ref{fig_gamma1}).
Composite dose delivery with the film horizontally placed for each plan was also performed (Figure~\ref{fig_gamma2}).

The films were calibrated by net-OD with actual dose delivered to the films, and fitted by a polynomial with degree of two. The measured dose distribution was compared with calculated dose distribution by dose profiles and gamma analysis.

For single field dose delivery (Figure~\ref{fig_gamma1}), the measured and calculated dose distributions match each other very well, which indicated that all dose calculation algorithm, measurement of radiation transmission, 3D printing, experimental dose delivery and film calibration worked appropriately. Compared to IMRT, the TVR-IMRT plan appears more blurred, which is evidence of the reduction in fluence TV. Gamma analysis indicated that TVR-IMRT has a slightly better fidelity than IMRT.

\begin{figure}[ht!]
	\centerline{
	\includegraphics[width=160mm]{fig/IMRT_newG_3pct_03mm_ref_calc0_local0_g3_filmno_25_disp}}
	\vspace{2mm}
	\centerline{
	\includegraphics[width=160mm]{fig/TVR_newG_3pct_03mm_ref_calc0_local0_g3_filmno_26_disp}}
	\caption{Gamma analysis of single field dose delivery. a: films. b and c: measured and calculated dose heat maps.
			d and e: profiles for measured dose (dotted line) and calculated dose (solid line).
			Green line denotes difference between measured and calculated.
			f: Gamma analysis of 3\%/0.3mm heat maps (Top: IMRT. Bottom: TVR-IMRT)
			}
	\label{fig_gamma1}
\end{figure}
\begin{figure}[ht!]
	\centerline{
	\includegraphics[width=160mm]{fig/IMRT_newG_3pct_03mm_ref_calc0_local0_g6_filmno_20_disp}}
	\vspace{2mm}
	\centerline{
	\includegraphics[width=160mm]{fig/TVR_newG_3pct_03mm_ref_calc0_local0_g6_filmno_19_disp}}
	\caption{Gamma analysis of composite dose delivery. Each column is the same as that in Figure~\ref{fig_gamma1}. (Top: IMRT. Bottom: TVR-IMRT)
			}
	\label{fig_gamma2}
\end{figure}

The gamma analysis for each plan composite dose delivery with 3\%/0.3mm criterion achieved a passing rate of 89.07\% and 95.37\% for IMRT and TVR-IMRT respectively, but the IMRT plan showed more distinct regions with significant deviation around the edges of the field (Figure~\ref{fig_gamma2}).

The gamma analysis of the composite dose meets the 90\% clinical passing threshold at 3\%/0.5mm for the IMRT plan, while TVR-IMRT passes at the more strict 3\%/0.2mm criteria (Table 3).
\begin{table}[ht]
	\caption{Gamma analysis results for IMRT and TVR-IMRT}
	{\small
		\begin{tabular}{r|r|r|r|r|r|r}
			\hline\hline
			  & beam1 & beam2 & beam3 & beam4 & beam5 & composite\cr
			\hline
			3\%/0.5mm: $\qquad$ IMRT&     99.60&    99.15&    99.69&    97.88&    96.59&    94.35\cr
			TVR-IMRT&        99.53&    98.81&    99.79&    99.52&    97.84&    97.51\cr
			\hline
			3\%/0.3mm: $\qquad$ IMRT& 99.23&    97.36&    98.33&    93.62&    94.05&    89.07    \cr
			TVR-IMRT&    99.09&    97.76&    99.04&    96.63&    95.28&    95.37         \cr
			\hline
			3\%/0.2mm: $\qquad$ IMRT&     97.51&    94.03&    93.80&    86.23&    87.76&    84.11    \cr
			TVR-IMRT&        98.09&    93.47&    96.79&    88.78&    88.96&    91.42    \cr
			\hline\hline
		\end{tabular}
	}
	\label{table_compare_gamma}
\end{table}

\section{Discussion and Conclusions}

 Preclinical research is essential to the continued development of RT technologies, fractionation schema, and knowledge of the underlying radiobiology which have all contributed to the improvement of clinical outcomes. The growing interest in small animal IMRT stems, in part, from the high degree of conformality and precision that is demanded in the clinic. In order to best predict the biological effects and outcomes of a procedure prior to use in the clinic, preclinical studies must be performed which can adhere to those same standards.

Dose painting and the associated sparing of neighboring healthy tissues has been proven to improve patient outcomes, so certain preclinical studies could benefit from these clinical-quality plans even though it would result in a smaller sample size.

We have shown that this method of planning is capable of producing the high fidelity dose distributions which make 3D-printed compensators an attractive method of small animal IMRT, with the plan generated by both IMRT and TVR-IMRT meeting all dose constraints and remaining within 1.2\% of one another for all of the plan evaluation metrics that we used. This was achieved while limiting the total variation, resulting in compensator modulation patterns that have shallow physical gradients relative to previous design methods\cite{gage2020small,yoon2020method}.

The reduction in steep gradients and variation produces a compensator design which enables higher fidelity printing, as shown by both the individual field dose heat maps and improved gamma passing rate for the TVR-IMRT plan (Figures~\ref{fig_gamma1}-\ref{fig_gamma2} and Table \ref{table_compare_gamma}). Furthermore, complex modulation patterns present in each beam angle may also result in a plan which is less robust with respect to small (sub degree, sub millimeter) variations in compensator positioning between dose deliveries, resulting in several areas of dose deviation in the composite dose films for the IMRT plan which were not present for the TVR-IMRT plan (Figure~\ref{fig_gamma2}). Both plans exhibited some deviation around the periphery, but it was more pronounced in the IMRT plan. The composite dose delivery of the TVR-IMRT plan passing the gamma analysis with 3\%/0.2mm criteria, compared with the IMRT plan only passing the 3\%/0.5mm criteria is further evidence that TVR-IMRT produces a compensator which is both more likely to print accurately and more robust to small positioning errors.

The TVR-IMRT method of 3D-printed compensator design was shown to provide an advantage over previous methods by reducing the sum of thicknesses used in the fluence map modulation region, resulting in lower cost and shorter print times, as well as reducing the total beam- on time for each exposure, while preserving or improving the dose distribution in the target as confirmed by DVH and gamma analysis.

In conclustion, here we have presented a technique utilizing total variation regularization (TVR) which improves upon the current design method for 3D-printed compensators for small animal IMRT. This is a subject which is garnering increasing interest due to the need to produce clinically analogous results in preclinical studies. TVR-IMRT limits the total variation across the modulated intensity pattern to be produced by the variable thickness compensator, resulting in more shallow gradients that facilitate high fidelity printing and can produce plans with improved gamma analysis passing rates. Furthermore, there is a significant reduction in total beam-on time to deliver a plan, total compensator thickness, 3D printer filament consumed, and print time when compared to previous methods.

\section*{References}
\addcontentsline{toc}{section}{\numberline{}References}
\vspace*{-15mm}

\bibliographystyle{ieeetr}


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

With an integrated luminosity of $150\,\mathrm{ab}^{-1}$ collected in $\approx 4$ years of running at centre-of-mass energies between 88 and 94\,GeV, FCC-ee~\cite{Benedikt:2651299} offers a unique opportunity to perform ultra-precise electroweak measurements of the Z resonance. More than $5\times 10^{12}$ Z decays, constituting the so-called TeraZ scenario, will be available for study. The statistical power of this sample is complemented by an extraordinary precision in the knowledge of the collision energy, of $\approx 100\,\mathrm{keV}$~\cite{R0polblondel}. Robust procedures to monitor other relevant beam collision parameters and the relative uncertainties between the energy points in Z lineshape scans~\cite{Blondel:2019jmp} are also an integral part of the physics program.

The increase in luminosity with respect to past LEP experiments translates into an increase of more than two orders of magnitude in statistical sensitivity. With such large potential improvement over previous measurements, a key question is how much the associated systematic uncertainties can be reduced,  and ultimately match the statistical uncertainties. Challenges arise both at the theoretical and experimental levels. New theoretical paths will have to be pursued in order to provide most precise predictions for experimental observables~\cite{Blondel:2018mad}, and detailed experimental studies will have to be performed to optimize accelerator and detector designs. Exploring new analysis strategies and observables to simultaneously reduce both theoretical and experimental uncertainties will be another essential component of the challenge.

An initial review of the FCC-ee potential regarding ultra-precise measurements at the Z pole can be found in Ref.~\cite{Abada:2019lih}. In the following we focus on some elements of the challenge that we consider relevant for success. The physics implications of the proposed program are expected to be deep. Regarding universal deviations, the new estimates of the oblique parameters $S$ and $T$ will constitute an early probe for Higgs compositeness or new interactions possibly occurring at the deca-TeV scale~\cite{Blondel:2021ema,Strategy:2019vxc}, i.e. one order of magnitude above currently explored scales. Deviations from universality related with these scales will be searched for through the precise measurement of the couplings of each fermion flavour to the Z. Separate access to the left- and right-handed components of the couplings will be available at FCC-ee even in the absence of polarised beams~\cite{Blondel:2019yqr}, as described below.

\section{General considerations for a successful TeraZ program}

Regarding cross-section measurements, relative statistical uncertainties are expected in the $1/\sqrt{N}=10^{-5}$-$10^{-6}$ range, where $N$ is the number of events selected in the decay channel under study. Similar statistical uncertainties, of the order of $\scriptstyle{\sqrt{(1-\mathcal{A}^2)/N}}$, are expected on the absolute value of measured asymmetries $\mathcal{A}$

At the theoretical front, and focusing on cross-section measurements, a limiting factor is the precision of the theoretical predictions used in luminosity measurements. Current studies using Bhabha scattering at low angle consider a relative accuracy of order $10^{-4}$ as a realistic target~\cite{Jadach:2018jjo}. The $\eegg$ process has also been suggested as an alternative channel to consolidate a $10^{-4}$ precision measurement, owing to the almost negligible size of the theoretical uncertainties of hadronic origin for this process~\cite{CarloniCalame:2019dom}. Cross-section ratio and asymmetry measurements are not affected by luminosity uncertainties. Besides the calculation of missing higher orders in these observables, reaching precisions of order $10^{-5}$ or better will probably require a deep change of philosophy in the measurement of lineshape parameters: inclusion of non-factorizable terms; interference effects between initial and final state radiation; redefinition of electroweak parameters at the amplitude level; fits of complete differential distributions using dedicated weighted Monte Carlo calculations, etc.~\cite{Blondel:2018mad} For the discussion below, and in order to better illustrate the main challenges, we nevertheless assume a simplified LEP-like strategy in the measurement of electroweak parameters~\cite{ALEPH:2005ab}.

Focusing on experimental aspects, a typical limiting factor for cross-section measurements is the systematic uncertainty on the acceptance determination. A $10^{-5}$ uncertainty, even in processes presenting a relatively smooth behaviour of the
angular distributions, implies a knowledge of the positions of the edges of sub-detectors at the $10\,\mu\mathrm{m}$ level over distances of the order of a meter. A first consequence is that detectors should be as homogeneous as possible. Such a precision is a realistic target given current tracking accuracy, but it demands dedicated efforts in terms of metrology, alignment, monitoring and designs able to ensure the stability of large detector volumes as a function of time. The challenge is even bigger for detectors located at very low polar angles and measuring differential cross sections with a $d\sigma/d\theta \propto 1/\sin\theta$ behaviour. For instance, a luminosity monitor located at 1\,m of the interaction point with an inner radius of $\approx 65$\,mm demands a $1\,\mu\mathrm{m}$ ($1\,\mu \mathrm{rad}$) precision in positioning, in order to reach $10^{-4}$ uncertainties~\cite{Benedikt:2651299}. Other requirements imposed by acceptance systematics are the uniformity in the detector response, redundant particle identification capabilities, beam stability and a detailed monitoring of the beam geometry conditions at the interaction point.

\section{Z lineshape determination and ppm/keV precision observables}
\label{section:xsections}

At FCC-ee Z lineshape scans are expected to provide a measurement of $\MZ$ with unprecedented precision, $\delta\MZ \approx 0.1\,\mathrm{MeV}$, i.e., $\approx 20$ times better than the present precision from LEP. The uncertainty is fully dominated by the uncertainty on the collision energy, which can be determined with $\approx 100\,\mathrm{keV}$ precision using resonant depolarisation of the transversely polarised beams~\cite{Blondel:2019jmp}. This method was already used at LEP and will be significantly improved at FCC-ee. There will be very frequent in-situ calibrations using non-colliding  pilot bunches, simultaneous with nominal collisions. Nevertheless, detailed studies of the differences between colliding and pilot bunches will be necessary. Also, the beam energy spread and the energy asymmetries between the two beams should be monitored via analysis of the longitudinal boosts of $\epem\to\mumu$ events~\cite{Blondel:2019jmp}.

The total Z width, $\GZ$, is directly connected with the width of the Z lineshape. A statistical precision of $4\,\mathrm{keV}$ is expected from a fit of the hadronic lineshape. The overall precision is dominated by the so-called ``point-to-point'' uncertainties, which correspond to systematic differences between their central $\rts$ values that are not 100\% positively correlated. The high statistics and the expected muon momentum resolution of tracking detectors ($\delta \left( 1/p_{\rm T} \right) \approx$ a few $10^{-5}$ for $p = 45$\,GeV at normal incidence) allow a quantification of these differences from the dimuon invariant mass distributions obtained at each collision energy. An uncertainty of $\approx 25\,\mathrm{keV}$ was obtained in preliminary studies~\cite{Blondel:2019jmp}, which translates into a precision improvement of two orders of magnitude over current LEP results.
re
The ratio between hadronic and leptonic cross sections, $\Rl$, is an essential observable for the extraction from the total width of the global and individual leptonic and hadronic partial widths. The global leptonic partial width is a direct test of new physics involving weak isospin violations, while the individual partial widths constitute a powerful test of lepton coupling universality in the neutral current. $\Rl$ also provides one of the most precise ways to measure $\alphasMZ$~\cite{Gomez-Ceballos:2013zzn,Benedikt:2651299}. Its measurement is independent of luminosity uncertainties, and therefore relative precisions below $10^{-4}$ can be contemplated. Improvements over LEP results by a factor ranging between $20$ and $100$, i.e., $\delta\Rl/\Rl = (1-5)\times 10^{-5}$ are expected. Leptonic and hadronic Z decays provide clean signatures at the Z pole, and are only affected by limited backgrounds at small visible mass (two-photon processes). According to LEP experience, the largest source of uncertainty could be the acceptance of the leptonic channels. Even with a rather hermetic detector and sufficiently redundant identification criteria, the edges of the tracking acceptance and the interplay with beam position and width parameters will have to be understood in detail. A precision in the position of these edges at the level of of $10\,\mu\mathrm{m}$ might be required.

A measurement of the relative invisible width $\Ginv/\Gl$, or equivalently of the number of neutrinos, $N_{\nu}= \Ginv/\Gn$, involves the measurement of the Z peak cross section~\cite{ALEPH:2005ab,Voutsinas:2019hwu,Janot:2019oyi}. The measurement is therefore limited by the precision in the measurement of the luminosity, as commented before. A $10^{-4}$ precision represents already one order of magnitude improvement with respect to LEP for a similar luminosity detector coverage. The study of radiative recoil ratios like $\sigma(\Pnu\APnu\gamma)/\sigma(\ell^+\ell^-\gamma)$ above the Z pole is also being considered as an alternative at a similar level of precision, dominated by statistical uncertainties~\cite{Abada:2019lih}.

\section{Precise measurements of the electroweak fermion couplings to the Z and of \texorpdfstring{$\alphaMZ$}{qed}}
\label{section:asymmetries}

One of the main targets of the electroweak program is the precise measurement of the chiral couplings to the Z for each individual fermion. The baseline FCC-ee proposal relies on a direct measurement of the Z polarisation induced by the beam particles with studies of the $\tau$ polarisation as a function of the polar angle. This approach suppresses the need of longitudinal beam polarisation, which would otherwise imply a non-negligible loss in luminosity.
The tau polarisation in Z decay is measured from the charged particle momentum distribution in the semi-leptonic decays $\tau \to {\rm e} \nu_{\rm e} \nu_\tau$, $\mu \nu_\mu \nu_\tau$ or in hadronic decays $\tau \rightarrow {\rm h} \nu_\tau$ where  $\rm h$ can be $\pi$, ${\rm K}$,  $\rho$, ${\rm K}^\ast$, ${\rm a}_1$, etc. Each channel has a different polarisation analysis power and must therefore be analysed independently. A clean separation between channels is also essential. The analysis of the $\tau$ polarisation dependence on the $\epem\to\tau^+\tau^-$ scattering angle $\theta$ gives access to both the tau and electron chiral coupling asymmetries ${\cal A}_\tau$ and ${\cal A}_{\rm e}$ independently
\begin{equation}
P(\cos\theta)   =  - \frac{  {\cal A}_\tau ( 1 + \cos^2 \theta)  + 2 {\cal A}_{\rm e} \cos\theta    } {  (1 + \cos^2 \theta) + 2 {\cal A}_{\rm e} {\cal A}_\tau \cos\theta  },
\end{equation}
The previous expression is valid at Born-improved level and in the massless lepton limit. $\Af$ is the chiral asymmetry parameter of the fermion $\Pf$ in its coupling to the Z.
The average $\tau$ polarisation, $\Poltau$, provides a direct measurement of $\Atau$: $\Poltau=-\Atau$, whereas the forward-backward polarisation asymmetry, $\AFBpol$, provides a direct measurement of $\Ae$, the induced Z polarisation in the $\epem$ collision: $\AFBpol=-\frac{3}{4}\Ae$. The remaining systematic uncertainty on $\Ae$ at LEP was originating from the limited knowledge of non-$\tau$ backgrounds~\cite{Heister:2001uh}. At FCC-ee, huge control samples will be available to reduce this component, but dedicated studies are still necessary to estimate the ultimately reachable precision.

The forward-backward asymmetry in the $\epem\to\ffbar$ process, $\AFBf = \frac{3}{4}\Ae\Af$ gives access to the chiral couplings of the fermion $\Pf$ when combined with tau polarisation studies. Chiral couplings can be trivially converted in measurements of vector and axial couplings or, alternatively, in a Born-improved spirit, in measurements of the effective $\rhof$ and $\sinf$ parameters for each individual fermion. Assuming universal deviations, $\sinf$ measurements become a measurement of the effective weak mixing angle $\sinSqEff$. The muon channel is particularly promising in this respect, aiming for a precision approaching the ultimate statistical sensitivity of $3\times 10^{-6}$ on $\AFBm$. At that level of precision, further theoretical studies of QED corrections and in particular of interference effects between initial and final-state radiation will be mandatory. Detailed studies of the beam parameters and of the polar angular resolution will be required as well~\cite{Blondel:2019jmp,Jadach:2018lwm}.

The extreme FCC-ee precision requirements demand an equivalent level of precision on the input parameters to theoretical predictions. One of these parameters is $\alphaMZ$, the electromagnetic coupling constant at the Z scale, which is the source of one of the dominant uncertainties in present fits. In practice, $\alphaMZ$ can be treated as another parameter to be determined in the Z lineshape running. This is supported by the study presented in Ref.~\cite{Janot:2015gjr}, which proposes a measurement of $\alphaMZ$ with a $3\times 10^{-5}$ relative precision, largely uncorrelated with other Z lineshape parameters. Most of the sensitivity lies in the linear dependence of the $\AFBm$ asymmetry with respect to the $\Pgamma-\PZ$ interference term around the Z peak. An appropriate combination of measurements at the energy points with the largest sensitivity, $\sqrt{s}=87.9$ and $94.3$~GeV, using one year of integrated luminosity should provide the required precision. Still, reaching the aimed precision will require the calculation of missing electroweak corrections of higher order, as well as more detailed studies of initial-final state interference effects.

\section{Heavy quark precision measurements}
\label{section:hq}

The cross-section ratios and asymmetries $\Rb\equiv \Gb/\Gh$, $\Rc\equiv \Gc/\Gh$, $\AFBb$, and $\AFBc$ are expected to be measured at FCC-ee with more than one order of magnitude better precision than at LEP/SLC. On top of the increase in the number of collected events, the statistical precision should be boosted with respect to LEP/SLC by the much higher b- or c-tagging efficiencies and purities obtained with better detectors and more powerful lifetime tagging techniques. Systematic uncertainties are expected to be significantly reduced by the use of exclusive B decays and enriched control samples in specific regions of phase space.

At FCC-ee, $\Rb$ and $\Rc$ will be likely measured with double-tagging techniques on both hemispheres of a $\bbbar$ or a $\ccbar$ event~\cite{ALEPH:2005ab}. This strategy allows a measurement independent of the knowledge of the tagging efficiency in the limit of negligible backgrounds and correlations between  hemispheres. At LEP, hemisphere correlations due to QCD effects (mostly hard gluon emission) and primary vertex determination were dominant sources of uncertainty. Gluon splitting increases the number of single tags in events with two light-flavour jets, and also constitutes a significant source of correlated uncertainty between experiments. The vertex precision of the new generation of detectors should contribute to reduce vertex correlations significantly.  Studies as a function of the acoplanarity between b-tagged jets can help reducing QCD correlations, and huge gluon splitting samples will become available for a precise understanding of this source of uncertainty. Nevertheless, coming studies are needed to quantify these improvements in more detail.

The bare forward-backward asymmetry of b quarks at the Z pole, $\AFBzerob$, is the electroweak observable that currently presents the largest deviation with respect to the standard model expectation in current fits~\cite{ALEPH:2005ab} ($\approx 3\sigma$ pull). An order-of-magnitude improved measurement at FCC-ee could thus become a clean signal of new physics if the deviation in the central value is confirmed. The world-average measurement is still dominated by statistical uncertainties ($\delta\AFBzerob({\rm stat.}) = 0.0016$), but is also affected by non-negligible systematic uncertainties ($\delta\AFBzerob({\rm syst.})=0.0007$). A fraction of it can be reduced at FCC-ee through dedicated studies on high-statistics control samples. A detailed analysis of the detector requirements to maximise flavour identification capabilities is also mandatory. Exclusive B decays can be exploited as well. For instance, about $10^8$ $\mathrm{B}^+$ decays, not affected by charm contamination or B-mixing effects, will be available at FCC-ee~\cite{Fabrizio-afbb}. The expected improvements in flavour tagging will be much more visible in the case of $\AFBzeroc$ measurements. In addition, the availability of exclusive, high purity, and large statistics D decay samples should provide a significant improvement in terms of precision compared with LEP measurements~\cite{ALEPH:2005ab}. There, most measurements were performed with inclusive or pseudo-inclusive techniques on samples with significant b-quark contamination. Dedicated studies are certainly necessary.

An irreducible source of uncertainty in the current estimate of $\AFBzerob$, fully correlated among experiments, is the presence of a QCD correction factor, of order $1-\alphas/\pi$, that accounts for the shift between the experimentally observed asymmetry and $\AFBzerob$~\cite{Djouadi-afbb,Altarelli:1992fs,ep-98}. The main role of the correction is to absorb the angular distortions due to final-state QCD radiation in the $\mathrm{Z}\to\mathrm{b\overline{b}(g)}$ decay. Recent re-evaluations of that uncertainty~\cite{dEnterria:2018jsx} based on modern parton shower tunes seem to be consistent with the initial estimates. New strategies to reduce or constrain experimentally the size of these uncertainties are being developed, for an initial target of $\delta\AFBzerob \approx 0.0001$~\cite{AlcarazMaestre:2020fmp}.

\section{Outlook}
\label{section:outlook}
Besides theory requirements, discussed in more detail in~\cite{Freitas:2019bre,Heinemeyer:2021rgq}, careful experimental studies with realistic detector descriptions are necessary to estimate the ultimate precision for FCC-ee electroweak measurements at the Z pole. In the particular case of cross section measurements, aiming for a precision of $\approx 10^{-4}$ is already imposing severe constraints on the design and tolerances of luminosity monitors. Regarding final states involving heavy flavours (tau, bottom, charm), an improvement of one order of magnitude with respect to previous LEP/SLC measurements may imply new constraints on the detector design and the development of new, more powerful tagging techniques. Let us note that the huge available statistics at the Z peak suggests that exclusive decay identification should be explored in more depth, as a complementary path to reduce systematic uncertainties in some particular cases. Addressing all these challenges is a critical step for success and, as such, one of the main objectives of the present electroweak physics program at FCC-ee.

\bibliographystyle{myutphys}