<EMAIL_ADDRESS>
# Modular Object-Oriented Games: A Task Framework for Reinforcement Learning,
Psychology, and Neuroscience

Nicholas Watters  Department of Brain and Cognitive Sciences, Massachusetts
Institute of Technology  Joshua Tenenbaum  Department of Brain and Cognitive
Sciences, Massachusetts Institute of Technology  Center for Brains, Minds and
Machines, MIT  Mehrdad Jazayeri  Department of Brain and Cognitive Sciences,
Massachusetts Institute of Technology  McGovern Institute of Brain Research,
MIT

###### Abstract

In recent years, trends towards studying simulated games have gained momentum
in the fields of artificial intelligence, cognitive science, psychology, and
neuroscience. The intersections of these fields have also grown recently, as
researchers increasing study such games using both artificial agents and human
or animal subjects. However, implementing games can be a time-consuming
endeavor and may require a researcher to grapple with complex codebases that
are not easily customized. Furthermore, interdisciplinary researchers studying
some combination of artificial intelligence, human psychology, and animal
neurophysiology face additional challenges, because existing platforms are
designed for only one of these domains. Here we introduce Modular Object-
Oriented Games, a Python task framework that is lightweight, flexible,
customizable, and designed for use by machine learning, psychology, and
neurophysiology researchers.

##  1 Introduction

In recent years, trends towards studying object-based games have gained
momentum in the fields of artificial intelligence, cognitive science,
psychology, and neuroscience. In artificial intelligence, interactive physical
games are now a common testbed for reinforcement learning (Mnih et al., 2013;
Sutton and Barto, 2018; Leike et al., 2017; François-Lavet et al., 2018) and
object representations are of particular interest for sample efficient and
generalizable AI (Battaglia et al., 2018; van Steenkiste et al., 2019; Greff
et al., 2020). In cognitive science and psychology, object-based games are
used to study a variety of cognitive capacities, such as planning, intuitive
physics, and intuitive psychology (Ullman et al., 2017; Chabris, 2017).
Developmental psychologists also use object-based visual stimuli to probe
questions about object-oriented reasoning in infants and young animals (Spelke
and Kinzler, 2007; Wood et al., 2020). In neuroscience, object-based computer
games have recently been used to study decision-making and physical reasoning
in both human and non-human primates (Fischer et al., 2016; McDonald et al.,
2019; Yoo et al., 2020; Rajalingham et al., 2021).

Furthermore, a growing number of researchers are studying tasks using a
combination of approaches from these fields. Comparing artificial agents with
humans or animals performing the same tasks can help constrain models of
human/animal behavior, generate hypotheses for neural mechanisms, and may
ultimately facilitate building more intelligent artificial agents (Lake et
al., 2017; Willke et al., 2019; Hassabis et al., 2017).

However, building a task that can be played by AI agents, humans, and animals
is a time-consuming undertaking because existing platforms are typically
designed for only one of these purposes. Professional game engines are
designed for human play and are often complex libraries that are difficult to
customize for training AI agents and animals. Reinforcement learning platforms
are designed for AI agents but are often too slow or inflexible for
neuroscience work. Existing psychology and neurophysiology platforms are too
limited to easily support complex interactive games.

In this work we offer a game engine that is highly customizable and designed
to support tasks that can be played by AI agents, humans, and animals.

##  2 Summary

The Modular Object-Oriented Games library is a general-purpose Python-based
platform for interactive games. It aims to satisfy the following criteria:

  * •

Usable for reinforcement learning, psychology, and neurophysiology. MOOG
supports DeepMind dm_env and OpenAI Gym (Brockman et al., 2016) interfaces for
RL agents and an MWorks interface for psychology and neurophysiology.

  * •

Highly customizable. Environment physics, reward structure, agent interface,
and more can be customized.

  * •

Easy to rapidly prototype tasks. Tasks can be composed in a single short file.

  * •

Light-weight and efficient. Most tasks run quickly, almost always faster than
100 frames per second on CPU and often much faster than that.

  * •

Facilitates procedural generation for randomizing task conditions every trial.

##  3 Intended Users

MOOG was designed for use by the following kinds of researchers:

  * •

Machine learning researchers studying reinforcement learning in
2.5-dimensional (2-D with occlusion) physical environments who want to quickly
implement tasks in Python.

  * •

Psychology researchers who want more flexibility than existing psychology
platforms afford.

  * •

Neurophysiology researchers who want to study interactive games yet still need
to precisely control stimulus timing.

  * •

Machine learning researchers studying unsupervised learning, particularly in
the video domain. MOOG can be used to procedurally generate video datasets
with controlled statistics.

MOOG may be particularly useful for interdisciplinary researchers studying AI
agents, humans, and animals (or some subset thereof) all playing the same
task.

##  4 Design

The core philosophy of MOOG is "one task, one file." Namely, a task should be
implemented with a single configuration file. This configuration file is a
short “recipe” for the task, containing as little substantive code as
possible, and should define a set of components to pass to the MOOG
environment. See Figure 1 for a schematic of these components.

Figure 1: Components of a MOOG environment. See main text for details.

A MOOG environment receives the following components (or callables returning
them) from the configuration file:

  * •

State. The state is a collection of sprites. Sprites are polygonal shapes with
color and physical attributes (position, velocity, angular velocity, and
mass). Sprites are 2-dimensional, and the state is 2.5-dimensional with
z-ordering for occlusion. The initial state can be procedurally generated from
a custom distribution at the beginning of each episode. The state is
structured in terms of layers, which helps hierarchical organization. See
state_initialization for procedural generation tools.

  * •

Physics. The physics component is a collection of forces that operate on the
sprites. There are a variety of forces built into MOOG (collisions, friction,
gravity, rigid tethers, …) and additional custom forces can also be used.
Forces perturb the velocity and angular velocity of sprites, and the sprite
positions and angles are updated with Newton’s method. See physics for more.

  * •

Task. The task defines the rewards and specifies when to terminate a trial.
See tasks for more.

  * •

Action Space. The action space allows the subject to control the environment.
Every environment step calls for an action from the subject. Action spaces may
impart a force on a sprite (like a joystick), move a sprite in a grid (like an
arrow-key interface), set the position of a sprite (like a touch-screen), or
be customized. The action space may also be a composite of constituent action
spaces, allowing for multi-agent tasks and multi-controller games. See
action_spaces for more.

  * •

Observers. Observers transform the environment state into a observation for
the subject/agent playing the task. Typically, the observer includes a
renderer producing an image. However, it is possible to implement a custom
observer that exposes any function of the environment state. The environment
can also have multiple observers. See observers for more.

  * •

Game Rules (optional). If provided, the game rules define dynamics or
transitions not captured by the physics. A variety of game rules are included,
including rules to modify sprites when they come in contact, conditionally
create new sprites, and control phase structure of trials (e.g. fixation phase
to stimulus phase to response phase). See game_rules for more.

Importantly, all of these components can be fully customized. If a user would
like a physics force, action space, or game rule not provided by MOOG, they
can implement a custom one, inheriting from the abstract base class for that
component. This can typically be done with only a few lines of code.

The modularity of MOOG facilitates code re-use across experimental paradigms.
For example, if a user would like to both collect behavior data from humans
using a continuous joystick and train RL agents with discrete action spaces on
the same task, they can re-use all other components in the task configuration,
only changing the action space.

For users interested in doing psychology or neurophysiology, we include an
example of how to run MOOG through MWorks, a platform with precise timing
control and interfaces for eye trackers, HID devices, electrophysiology
software, and more.

##  5 Example Tasks

Figure 2: Example tasks. Time-lapse images of four example tasks. Left-to-
right: (i) Pong - The subject aims to catch the yellow ball with the green
paddle, (ii) Red-Green - The subject tries to predict whether the blue ball
with contact the red square or the green square, (iii) Pac-Man - The subject
moves the green agent to catch yellow pellets while avoiding the red ghosts,
(iv) Collisions - the green agent avoids touching the bouncing polygons.

See the example_configs for a variety of task config files. Four of those are
shown in Figure 2. See the demo documentation for videos of them all and
instructions for how to run them with a Python gui.

The MOOG codebase contains libraries of options for each of the components in
Section 4, so implementing a task involves only combining the desired
ingredients and feeding them to the environment. For an example, the following
code fully implements a navigate-to-goal task, where the subject must move an
agent via a joystick action space to a goal location:

⬇

1"""Navigate-to-goal task."""

2

3import collections

4from moog import action_spaces, environment, observers, physics, sprite,
tasks

5

6# Initial state is a green agent in the center and red goal in the corner

7def state_initializer():

8 # Goal is a red square with side-length 0.1 at position (0.1, 0.1)

9 # Color channels are [c0, c1, c2] arguments

10 goal = sprite.Sprite(x=0.1, y=0.1, shape=’square’, scale=0.1, c0=255)

11 # Agent is a green circle at position (0.5, 0.5)

12 agent = sprite.Sprite(x=0.5, y=0.5, shape=’circle’, scale=0.1, c1=255)

13 state = collections.OrderedDict([

14 (’goal’, [goal]),

15 (’agent’, [agent]),

16 ])

17 return state

18

19# Physics is a drag force on the agent to limit velocity

20phys = physics.Physics((physics.Drag(coeff_friction=0.25), ’agent’))

21

22# Task gives a reward 1 when the agent reaches the goal, and a new trial
begins

23# 5 timesteps later.

24task = tasks.ContactReward(1., ’agent’, ’goal’, reset_steps_after_contact=5)

25

26# Action space is a joystick with maximum force 0.01 arena widths per
timestep^2

27action_space = action_spaces.Joystick(0.01, ’agent’)

28

29# Observer is an image renderer

30obs = observers.PILRenderer(image_size=(256, 256))

31

32# Create the environment, ready to play!

33env = environment.Environment(

34 state_initializer=state_initializer,

35 physics=phys,

36 task=task,

37 action_space=action_space,

38 observers={’image’: obs},

39}

This is an extremely simple task, but by complexifying the state initializer
and adding additional forces and game rules, a wide range of complex tasks can
be implemented with few lines of code.

##  6 Limitations

Users should be aware of the following limitations of MOOG before choosing to
use it for their research:

  * •

Not 3D. MOOG environments are 2.5-dimensional, meaning that they render in
2-dimensions with z-ordering for occlusion. MOOG does not support 3D sprites.

  * •

Very simple graphics. MOOG sprites are monochromatic polygons. There are no
textures, shadows, or other visual effects. Composite sprites can be
implemented by creating multiple overlapping sprites, but still the graphics
complexity is very limited. This has the benefit of a small and easily
parameterizable set of factors of variation of the sprites, but does make MOOG
environments visually unrealistic.

  * •

Imperfect collisions. MOOG’s collision module implements Newtonian rotational
mechanics, but it is not as robust as professional physics engines (e.g. can
be unstable if object are moving very quickly and many collisions occur
simultaneously).

##  7 Related Software

Professional game engines (e.g. Unity and Unreal) and visual reinforcement
learning platforms (e.g. DeepMind Lab (Beattie et al., 2016), Mujoco (Todorov
et al., 2012), and VizDoom) are commonly used in the machine learning field
for task implementation. While MOOG has some limitations compared to these
(see above), it does also offer some advantages:

  * •

Python. MOOG tasks are written purely in Python, so users who are most
comfortable with Python will find MOOG easy to use.

  * •

Procedural Generation. MOOG facilitates procedural generation, with a library
of compositional distributions to randomize conditions across trials.

  * •

Online Simulation. MOOG supports online model-based RL, with a ground truth
simulator for tree search.

  * •

Psychophysics. MOOG can be run with MWorks, a psychophysics platform.

  * •

Speed. MOOG is fast on CPU. While the speed depends on the task and rendering
resolution, MOOG typically runs at 200fps with 512x512 resolution on a CPU,
much faster than DeepMind Lab and Mujoco and at least as fast as Unity and
Unreal.

Python-based physics simulators, such PyBullet (Coumans and Bai, 2016–2019)
and Pymunk, are sometimes used in the psychology literature. While these offer
highly accurate collision simulation, MOOG offers the following advantages:

  * •

Customization. Custom forces and game rules can be easily implemented in MOOG.

  * •

Psychophysics, Procedural Generation, and Online Simulation, as described
above.

  * •

RL Interface. A task implemented in MOOG can be used out-of-the-box to train
RL agents, since MOOG is Python-based and has DeepMind dm_env and OpenAI Gym
interfaces.

Psychology and neurophysiology researchers often use platforms such as
PsychoPy (Peirce et al., 2019), PsychToolbox (Kleiner et al., 2007), and
MWorks. These allow precise timing control and coordination with eye trackers
and other controllers. MOOG can interface with MWorks to leverage all of those
features and offers the following additional benefits:

  * •

Flexibility. MOOG offers a large scope of interactive tasks. Existing
psychophysics platforms are not easily customized for game-like tasks, action
interfaces, and arbitrary object shapes.

  * •

Physics. Existing psychophysics platforms do not have built-in physics, such
as forces, collisions, etc.

  * •

RL Interface, as described above.

### Acknowledgments

We thank Chris Stawarz and Erica Chiu for their contributions to the codebase.
We also thank Ruidong Chen, Setayesh Radkani, and Michael Yoo for their
feedback as early users of OOG.

## References

  * Battaglia et al. (2018) P. W. Battaglia, J. B. Hamrick, V. Bapst, A. Sanchez-Gonzalez, V. Zambaldi, M. Malinowski, A. Tacchetti, D. Raposo, A. Santoro, R. Faulkner, C. Gulcehre, F. Song, A. Ballard, J. Gilmer, G. Dahl, A. Vaswani, K. Allen, C. Nash, V. Langston, C. Dyer, N. Heess, D. Wierstra, P. Kohli, M. Botvinick, O. Vinyals, Y. Li, and R. Pascanu.  Relational inductive biases, deep learning, and graph networks, 2018.
  * Beattie et al. (2016) C. Beattie, J. Z. Leibo, D. Teplyashin, T. Ward, M. Wainwright, H. Küttler, A. Lefrancq, S. Green, V. Valdés, A. Sadik, J. Schrittwieser, K. Anderson, S. York, M. Cant, A. Cain, A. Bolton, S. Gaffney, H. King, D. Hassabis, S. Legg, and S. Petersen.  Deepmind lab.  _CoRR_ , abs/1612.03801, 2016.  URL http://arxiv.org/abs/1612.03801.
  * Brockman et al. (2016) G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba.  Openai gym.  _CoRR_ , abs/1606.01540, 2016.  URL http://arxiv.org/abs/1606.01540.
  * Chabris (2017) C. F. Chabris.  Six suggestions for research on games in cognitive science, 2017.  URL https://doi.org/10.1111/tops.12267.
  * Coumans and Bai (2016–2019) E. Coumans and Y. Bai.  Pybullet, a python module for physics simulation for games, robotics and machine learning.  http://pybullet.org, 2016–2019.
  * Fischer et al. (2016) J. Fischer, J. G. Mikhael, J. B. Tenenbaum, and N. Kanwisher.  Functional neuroanatomy of intuitive physical inference.  _Proceedings of the national academy of sciences_ , 113(34):E5072–E5081, 2016.
  * François-Lavet et al. (2018) V. François-Lavet, P. Henderson, R. Islam, M. G. Bellemare, and J. Pineau.  An introduction to deep reinforcement learning.  _Foundations and Trends® in Machine Learning_ , 11(3-4):219–354, 2018.  ISSN 1935-8245.  10.1561/2200000071.  URL http://dx.doi.org/10.1561/2200000071.
  * Greff et al. (2020) K. Greff, S. van Steenkiste, and J. Schmidhuber.  On the binding problem in artificial neural networks.  _arXiv preprint arXiv:2012.05208_ , 2020.
  * Hassabis et al. (2017) D. Hassabis, D. Kumaran, C. Summerfield, and M. Botvinick.  Neuroscience-inspired artificial intelligence.  _Neuron_ , 95(2):245–258, 2017.
  * Kleiner et al. (2007) M. Kleiner, D. Brainard, D. Pelli, A. Ingling, R. Murray, and C. Broussard.  What’s new in psychtoolbox-3, 2007.
  * Lake et al. (2017) B. M. Lake, T. D. Ullman, J. B. Tenenbaum, and S. J. Gershman.  Building machines that learn and think like people.  _Behavioral and brain sciences_ , 40, 2017.
  * Leike et al. (2017) J. Leike, M. Martic, V. Krakovna, P. A. Ortega, T. Everitt, A. Lefrancq, L. Orseau, and S. Legg.  Ai safety gridworlds, 2017.
  * McDonald et al. (2019) K. R. McDonald, W. F. Broderick, S. A. Huettel, and J. M. Pearson.  Bayesian nonparametric models characterize instantaneous strategies in a competitive dynamic game.  _Nature communications_ , 10(1):1–12, 2019.
  * Mnih et al. (2013) V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller.  Playing atari with deep reinforcement learning.  _arXiv preprint arXiv:1312.5602_ , 2013.
  * Peirce et al. (2019) J. Peirce, J. R. Gray, S. Simpson, M. MacAskill, R. Höchenberger, H. Sogo, E. Kastman, and J. K. Lindeløv.  PsychoPy2: Experiments in behavior made easy.  _Behavior Research Methods_ , 51(1):195–203, 2019\.  ISSN 1554-3528.  10.3758/s13428-018-01193-y.  URL https://doi.org/10.3758/s13428-018-01193-y.
  * Rajalingham et al. (2021) R. Rajalingham, A. Piccato, and M. Jazayeri.  The role of mental simulation in primate physical inference abilities.  _bioRxiv_ , 2021.  10.1101/2021.01.14.426741.  URL https://www.biorxiv.org/content/early/2021/01/17/2021.01.14.426741.
  * Spelke and Kinzler (2007) E. S. Spelke and K. D. Kinzler.  Core knowledge.  _Developmental science_ , 10(1):89–96, 2007.
  * Sutton and Barto (2018) R. S. Sutton and A. G. Barto.  _Reinforcement learning: An introduction_.  MIT press, 2018.
  * Todorov et al. (2012) E. Todorov, T. Erez, and Y. Tassa.  Mujoco: A physics engine for model-based control.  In _IROS_ , pages 5026–5033. IEEE, 2012.  ISBN 978-1-4673-1737-5.  URL http://dblp.uni-trier.de/db/conf/iros/iros2012.html#TodorovET12.
  * Ullman et al. (2017) T. D. Ullman, E. Spelke, P. Battaglia, and J. B. Tenenbaum.  Mind games: Game engines as an architecture for intuitive physics.  _Trends in cognitive sciences_ , 21(9):649–665, 2017.
  * van Steenkiste et al. (2019) S. van Steenkiste, K. Greff, and J. Schmidhuber.  A perspective on objects and systematic generalization in model-based rl.  _arXiv preprint arXiv:1906.01035_ , 2019.
  * Willke et al. (2019) T. L. Willke, S. B. M. Yoo, M. Capotă, S. Musslick, B. Y. Hayden, and J. D. Cohen.  A comparison of non-human primate and deep reinforcement learning agent performance in a virtual pursuit-avoidance task.  _bioRxiv_ , page 567545, 2019.
  * Wood et al. (2020) J. N. Wood, D. Lee, B. W. Wood, and S. M. Wood.  Reverse engineering the origins of visual intelligence.  2020\.
  * Yoo et al. (2020) S. B. M. Yoo, J. C. Tu, S. T. Piantadosi, and B. Y. Hayden.  The neural basis of predictive pursuit.  _Nature neuroscience_ , 23(2):252–259, 2020.


# Manually driven harmonic oscillator

M N S Silva and J T Carvalho-Neto  Departamento de Ciências da Natureza,
Matemática e Educação, Universidade Federal de São Carlos, Caixa Postal 153,
13600-970, Araras, SP, Brasil<EMAIL_ADDRESS>
###### Abstract

Oscillations and resonance are essential topics in physics that can be
explored theoretically and experimentally in the classroom or teaching
laboratory environments. However, one of the main challenges concerning the
experimental study of resonance phenomena via forced oscillations is the
control of the oscillation frequency, which demands an electronic circuit or a
fine tuned coupled mechanical system. In this work, we demonstrate that, in
what concerns the physics teaching, such demanding accessories are not
necessary. The forced oscillations can be implemented by the teacher’s hand
guided by an oscillating circle displayed in a web application loaded in a
smartphone. The oscillations are applied to an ordinary spiral toy.
Qualitative, as well quantitative, proposals are explored in this work with
excellent results.

  * November 2019

##  1 Introduction

Oscillations and waves are essential contents for understanding various
natural phenomena ranging from basic physics to applied sciences.
Specifically, the harmonic oscillators provide great explanatory power as they
have analytical solutions in various damping and excitation regimes [1].

In addition to the theoretical description of these systems, didactic and
illustrative experiments can be performed in the classroom environment,
providing a fairly complete description of oscillatory phenomena [2].
Regarding the damped and driven harmonic oscillator, the electromagnetic and
electromechanical systems are the most versatile to work with and have many
technological applications. Consequently, they are greatly explored for the
teaching of forced oscillations and resonance [2, 3, 4]. Nowadays,
technologies such as smartphone accelerometers have been widely used to study
mechanical oscillators [5, 6, 7, 8]. Optical sensors and counters are also
employed, some of which use reasonably complex electronic circuits [9, 10,
11].

There are few purely mechanical implementations for the teaching of the driven
oscillator. As an example, in [12] the author uses the oscillations of a heavy
pendulum to drive oscillations on a lighter one. Usually, electrical motors
are coupled to spring-mass or pendulum systems in order to drive the
oscillations [13, 14]. These devices, while ingenious and some easy to
construct, require time-consuming experimental preparations that can inhibit
widespread use in the classroom.

In this work, our purpose is to demonstrate that the excitation and reading of
a purely mechanical driven harmonic oscillator aimed at teaching oscillations
and resonance can be fairly well performed with a simple spiral toy [15, 16]
driven by the own teacher’s hand. Moreover, with just few more add-ons, it is
possible to quantitatively explore the oscillator amplitude response as
function of the external force frequency in different damping conditions. The
only necessary technological resource is a smartphone loaded with our open
source oscillator web application, which serves as a reference in swinging the
spiral toy. A multimedia projector is useful just in case of performing
quantitative analysis for the whole classroom.

##  2  Materials

The materials used in this work consisted of a spiral toy and a web
application called harmetronome [17] that we developed specifically for this
work (the name makes reference to the metronome that helps musicians to keep
the pace at a chosen beat rate). The web application consists of a circle that
executes simple harmonic oscillation. The oscillation frequency or period can
be set to a specific value or can be swiped in a predefined range.

The spiral toy consisted of a plastic helical spring with 7.6 cm in external
diameter and weighting 0.84 g per loop. Optionally, metal nuts and plastic
discs can be attached to the end of the spiral toy in order to change its mass
and damping factor, respectively.

##  3 Applications

###  3.1  Qualitative demonstration of resonance phenomenon

This is the simplest propose. The mode of operation is illustrated in figure 1
[18]. The teacher holds the smartphone (loaded with the harmetronome app) with
his hand and one end of the spiral toy with the other hand. Both the
smartphone display and the spiral toy must be well seen by the students. Then,
the hand that holds the spiral toy tries to keep in pace with the oscillating
circle being displayed in the smartphone. With such scheme we want to
implement a driven harmonic oscillator in order to explore the qualitative
aspects related to harmonic oscillations and resonance. For that, we propose
the following didactic sequence:

  1. 1.

The teacher presents the concepts associated with the natural frequency of
oscillation of a system and resonance phenomenon.

  2. 2.

In order to give a simple example of such a system, he keeps one end of the
spiral toy motionless with his own hand and puts it to oscillate in the
vertical position by displacing and releasing the bottom end.

  3. 3.

The students are asked to count a number of oscillations of the spiral toy
while the teacher measures the elapsed time with the aid of the smartphone
chronometer. With that measurement, they find the empirical value for the
system natural period of oscillation.

  4. 4.

In what follows, the teacher chooses a series of periods of oscillation to be
displayed in the harmetronome app: some above, some bellow, and one at the
just measured natural period. For each case, the teacher applies the forced
oscillations on the spiral toy as described in the earlier paragraph.

Figure 1:  Snapshots of the forced oscillations on the spiral toy (a) below
resonance at $\omega=2.6$ rad/s, (b) on resonance at $\omega=4.6$ rad/s, and
(c) above resonance at $\omega=6.6$ rad/s.

Using the procedure 3 described above, we measured a natural frequency of
oscillation of 4.6 rad/s. Therefore, we choose the following frequencies for
the manually driven oscillations: 2.6 rad/s, 4.6 rad/s, and 6.6 rad/s. Figures
1(a) to 1(c) show snapshots took at the moment that the spring bottom end was
at its lowest position for each oscillation frequency. The red line was drew
to indicate the spring equilibrium position. That is, the spiral toy bottom
end oscillates symmetrically below and above the red line. An unloaded spiral
toy with 26 and one fourth loops was used in this demonstration.

As can be seen in figure 1, the expected behavior is very well reproduced,
despite the inherently inaccurate movement performed by the teacher’s hand.
Bellow and above the natural frequency, the spring oscillation amplitudes are
smaller than the amplitude at the natural frequency. Additionally, at
$\omega=2.6$ rad/s we can see that the lowest spring position coincides with
the harmetronome lowest postion, which corresponds to a zero phase difference
between them. On the contrary, at $\omega=6.6$ rad/s, when the spring bottom
end is in the lowest position, the harmetronome is in its highest position,
which corresponds to a 180∘ phase difference between them. On resonance, at
$\omega=4.6$ rad/s, the harmetronome is in its middle position, which
corresponds to a 90∘ phase difference. It is important to note that since the
damping factor is significantly lower than the natural frequency, the
resonance frequency and the natural frequency are very close to each other.

###  3.2  Determination of the resonance curve

It is possible to quantitatively study the oscillator response to an external
force by running the harmetronome web application on a computer and projecting
it on the wall by means of a multimedia projector. Besides the oscillating
circle, there is the option to show a vertical scale against which the shadow
of a hanging mass is project and its position can be read by the students.
Therefore, as the teacher swings the top end of the spiral toy in pace with
the harmetronome circle, the students take note of the minimum and maximum
position reached by the hanging mass for each frequency of oscillation choose
by the teacher. At the end, the amplitude of oscillation for each frequency of
the external force is calculated by taking the difference between the maximum
and minimum positions.

In figure 2 [19] is shown a picture of the teacher executing the experimental
demonstration where the harmetronome application is projected on the
blackboard. In order to make the experiment realization more comfortable for
the teacher, the spiral toy is attached to the top of a toy car, which in turn
is moved by the teacher against the blackboard and a straight piece of wood.
The piece of wood helps to constrain the spring movement along a straight
line.

Figure 2:  Quantitative demonstration of the driven harmonic oscillator. With
the left hand, the teacher holds a white piece of wood to vertically constrain
the movement executed by its right hand. The shadow of the hanging washer can
be seen against the projected scale and the oscillating circle to the left of
the scale.

With the purpose to demonstrate the effectiveness of this propose, we
performed this demonstration to 34 students of a first year class of a
Brazilian undergraduate physics teaching course. A piece of six and a half
loops of the spiral toy was used. Additionally, a mass of 16.5 g and a 78.5
cm2 plastic disc were attached to the spring end in order to lower the
resonance frequency and the maximum oscillation amplitude, respectively. This
procedure was important to avoid the spring reaching its minimum elongation.
The plastic disc was cut from the cover of a polypropylene folder.

Figure 3(a) shows the graphic containing the resonance curves recorded by all
students. Figure 3(b) shows the graphic of the resonance curve obtained by the
average of all curves shown in figure 3(a). The amplitude of oscillation $A$
of the harmonic oscillator as function of the external force frequency
$\omega$ is given by [1]:

$A(\omega)=\frac{A_{0}}{\sqrt{\left(\omega_{0}^{2}-\omega^{2}\right)^{2}+4\omega^{2}\beta^{2}}}\;,$
(1)

where $A_{0}=F_{0}/m$, $\beta=b/2m$, $\omega_{0}=\sqrt{k/m}$, $F_{0}$ is the
external force amplitude, $m$ is the harmonic oscillator mass, and $b$ is the
drag coefficient. The continuous curve in figure 3(b) corresponds to the
fitting of equation (1) to the students measurements.

Figure 3: (a) Experimental measurements of the spiral toy resonance
demonstration shown in figure 2. Each continuous curve was measured by one of
the 34 students. (b) The average resonance curve taken from the data of the
individual students. The dots correspond to the average values and the error
bars to the standard deviations.

As result, the following parameters were obtained: $\omega_{0}=(5.97\pm 0.03)$
rad/s, $\beta=(0.53\pm 0.03)$ rad/s, and $A_{0}=(11.3\pm 0.4)$ m/s2. The
correlation coefficient of the model to the experimental data was $r=0.991$.

###  3.3 Application in the physics lab

Our proposal can also be applied to the physics teaching laboratory, where
groups of two or three students can realize the resonance experiment
separately. In such case, the spiral toy can oscillate around a rigid tube
(e.g. PVC water tube) to which a cloth tape measure is attached. The top end
of the spiral toy is fixed to the tube. While one student oscillates the tube
in pace with the harmetronome, the other students take note of the minimum and
maximum amplitudes reached by the spring as read on the cloth tape. A picture
of such configuration is shown in figure 4. The purpose of the rigid tube is
to constrain the coil movement along a straight line and to hold the scale
given by the cloth tape.

Figure 4:  Experiment to be performed in groups of students. The spiral toy is
fixed to the top of the guiding white tube. The oscillations are read against
the cloth tape in red. The pink plastic ring and the metal washers attached to
it change the oscillator damping factor and mass, respectively.

The damping forces and masses can be changed by attaching different plastic
rings and metal washers or nuts to the bottom end of the spiral toy. To
illustrate the use of this setup, we performed the resonance experiment for
the six configurations showed in table 1. The same piece of six and a half
loops spiral toy of section 3.2 was used. The resulting resonance curves are
shown in figure 5(a) to 5(f). The experimental data is more noisy than the
data of figure 3(b) because we did just one measurement for each
configuration. However, it could be greatly improved by performing more
repetitions and taking their average values.

Figure 5:  Experimental data obtained with the setup of figure 4 for different
values of effective mass $m_{\mathrm{eff}}$ and ring area $r_{\mathrm{A}}$.
The continuous curves correspond to the fitting using equation (1). The fitted
parameters are listed at table 1.

We also measured the spring elastic constant, $k=0.57\pm 0.01$ N/m, which was
obtained using the same apparatus by reading each new spring equilibrium
position as function of the attached masses.

In table 1 we compare the undamped natural frequencies $\omega_{0}$ obtained
by the relation $w_{0}=\sqrt{k/m_{\mathrm{eff}}}$ with this same parameter
obtained from the fitting of equation (1) to the experimental data in figures
5(a) to 5(f). Since the attached masses, $m$, are in the order of the spiral
toy mass ($m_{\mathrm{s}}=5.9$ g), it was necessary to add a mass correction
factor equals to $0.33m_{\mathrm{s}}$, such that the effective mass was
calculated as $m_{\mathrm{eff}}=m+0.33m_{\mathrm{s}}$ [20, 21]. As can be seen
in table 1, the natural frequencies obtained from the two approaches were
reasonably close to each other, presenting deviations from 2% to 7%. Another
correlation that can be observed is between the drag coefficient $b$ and the
plastic ring area $r_{\mathrm{A}}$.

Figure 5 | $m_{\mathrm{eff}}$ (g) |  ring area (cm2) |  $\omega_{0}=\sqrt{k/m_{\mathrm{eff}}}$ (rad/s) | From fitting of equation (1)
---|---|---|---|---
| $\omega_{0}$
---
(rad/s)
| $b$
---
(g/s)
(a) | 13.55 | 79 | 6.50 | 6.34 | 9.7
(b) | 20.05 | 79 | 5.34 | 5.67 | 13
(c) | 9.25 | 79 | 7.86 | 7.38 | 9.5
(d) | 11.05 | 119 | 7.19 | 6.67 | 20
(e) | 15.35 | 119 | 6.10 | 6.01 | 15
(f) | 21.85 | 119 | 5.12 | 5.42 | 18
Table 1: Driven harmonic oscillator parameters used in the experimental
configuration depicted in figure 4.

##  4 Conclusions

We have shown that a manually driven harmonic oscillator can be explored in
various levels of physics teaching.

In the qualitative approach, the amplitude of oscillation and phase difference
behavior as one approaches the system’s natural frequency was clearly
observed.

A quantitative proposal involving the construction of a resonance curve was
applied to a real class of 34 students and the experimental data thus obtained
was fitted by the theoretical model with a correlation coefficient of 0.991.
One of the main difficulties in this approach is the realization of many
different oscillations by a single person (e.g. the teacher) whose arm can get
tired. In this case, the teacher may take turns with other students.

At last, an implementation using a rigid tube as a guide to the spiral toy was
proposed to be applied to a physics teaching laboratory, where each group of
students can make its own quantitative analysis probing the forced damped
harmonic oscillator in different conditions.

## References

## References

  * [1] Thornton S T and Marion J B 2003 Classical Dynamics of Particles and Systems (Cengage Learning) ISBN 0534408966
  * [2] MacLeod A 1969 Physics Education 4 272
  * [3] Jarvis W 1975 Physics Education 10 370
  * [4] Jones C C 1995 American Journal of Physics 63 232–236
  * [5] Kuhn J and Vogt P 2012 The Physics Teacher 50 504–505
  * [6] Carlos Castro-Palacio J, Velázquez-Abad L, Giménez M H and Monsoriu J A 2013 American Journal of Physics 81 472–475
  * [7] Tuset-Sanchis L, Castro-Palacio J C, Gómez-Tejedor J A, Manjón F J and Monsoriu J A 2015 Physics Education 50 580
  * [8] Giménez M H, Salinas I, Monsoriu J A and Castro-Palacio J C 2017 The Physics Teacher 55 424–425
  * [9] Toyotomi S 1994 Physics Education 29 382
  * [10] Ng T and Ang K 2005 American Journal of Physics 73 793–795
  * [11] Wadhwa A 2009 Physics Education 44 70
  * [12] Ong P 1972 Physics Education 7 391
  * [13] Boving R et al. 1983 Physics Education 18 275–76
  * [14] Goncalves A, Cena C and Bozano D 2017 Physics Education 52 043002
  * [15] Pendrill A M and Williams G 2005 Physics Education 40 527–533 URL https://doi.org/10.1088/0031-9120/40/6/003
  * [16] Pendrill A M and Eager D 2014 Physics Education 50 64
  * [17] Silva M N S and Carvalho-Neto J T 2019 Harmetronome web application https://harmetronome.glitch.me/
  * [18] Silva M N S and Carvalho-Neto J T 2019 Qualitative manually driven harmonic oscillator demonstration https://youtu.be/Q9A-8cXByHk
  * [19] Silva M N S and Carvalho-Neto J T 2019 Quantitative manually driven harmonic oscillator demonstration https://youtu.be/74YtezSDFRo
  * [20] Fox J and Mahanty J 1970 American Journal of Physics 38 98–100
  * [21] Rodriguez E E and Gesnouin G A 2007 The Physics Teacher 45 100–103


# Estimating the Feasibility of 21cm-Ly$\alpha$ Synergies using the Hydrogen
Epoch of Reionization Array

Tyler A. Cox,1,2 Daniel C. Jacobs,2 Steven G. Murray,2
1Department of Astronomy, University of California, Berkeley, CA
2School of Earth and Space Exploration, Arizona State University, Tempe, AZ
E-mail<EMAIL_ADDRESS>
(Accepted 2022 February 15. Received 2022 February 14; in original form 2021
March 30)

###### Abstract

Cross-correlating 21cm and Ly$\alpha$ intensity maps of the Epoch of
Reionization (EoR) promises to be a powerful tool for exploring the properties
of the first galaxies. Next-generation intensity mapping experiments such as
the Hydrogen Epoch of Reionization Array (HERA) and SPHEREx will individually
probe reionization through the power spectra of the 21cm and Ly$\alpha$ lines
respectively, but will be limited by bright foregrounds and instrumental
systematics. Cross-correlating these measurements could reduce systematics,
potentially tightening constraints on the inferred astrophysical parameters.
In this study, we present forecasts of cross-correlation taking into account
the effects of exact uv-sampling and foreground filtering to estimate the
feasibility of HERAxSPHEREx making a detection of the 21cm-Ly$\alpha$ cross-
power spectrum. We also project the sensitivity of a cross-power spectrum
between HERA and the proposed next-generation Cosmic Dawn Intensity Mapper. By
isolating the sources of uncertainty, we explore the impacts of experimental
limitations such as foreground filtering and Ly$\alpha$ thermal noise
uncertainty have on making a detection of the cross-power spectrum. We then
implement this strategy in a simulation of the cross-power spectrum and
observational error to identify redshifts where fiducial 21cmFAST models
predict the highest signal-to-noise detection ($z\sim 8$). We conclude that
detection of the SPHEREx-HERA cross-correlation will require an optimistic
level of 21cm foreground filtering, as well as deeper thermal noise
integrations due to a lack of overlapping sensitive modes but for CDIM with
its larger range of scales and lower noise forecast detection levels, may be
possible even with stricter 21cm foreground filtering.

###### keywords:

cosmology: dark ages, reionization, first stars, galaxies: high-redshift,
instrumentation: interferometers

††pubyear: 2021

##  1 Introduction

The Epoch of Reionization (EoR) marks a major phase transition in the history
of the Universe. During this time period, the first stars and galaxies, formed
from density fluctuations seeded by inflation, produced X-ray and ultraviolet
photons, which heated and ionized the neutral gas around them. Gradually,
these first luminous sources ionized all of the neutral hydrogen gas in the
intergalactic medium (IGM) around them, transitioning the ionization state of
the Universe from completely neutral to fully ionized (see Barkana & Loeb 2001
for a review). Only a limited number of observations have been made which
constrain reionization timing and the properties of the sources driving its
progression. Observations of the Gunn-Peterson trough (Gunn & Peterson, 1965)
in the spectra of high-redshift quasars have been used to place the conclusion
of reionization at $z\sim 6$ (Fan et al., 2006). Assuming a model in which
reionization took place instantaneously, the integrated Thomson optical depth
of cosmic microwave background photons has also been used to imply a
reionization redshift of $z\sim 7.68$ (Planck Collaboration et al. 2020a, b).
These constraints suggest that the bulk of reionization took place from $z\sim
6-10$, yet little is known about the topology of reionization, the properties
of the first structures, or the spatial or temporal evolution of the ionized
bubbles surrounding the first luminous sources.

New techniques are currently being developed to help shed light on this period
of time. One of the most promising techniques for directly observing the EoR
is intensity mapping of the 21cm hyperfine transition of neutral hydrogen (for
a review of 21cm cosmology, see Furlanetto et al. 2006). The potential for
understanding reionization through the 21cm line has spawned a number of
experiments over the past decade, including PAPER111Precision Array for
Probing the Epoch of Reionization (Kolopanis et al., 2019), LOFAR222Low
Frequency Array(Gehlot et al., 2019), the GMRT333Giant Metre Wave Radio
Telescope (Paciga et al., 2013), and MWA444Murchison Widefield Array (Li et
al., 2019), which have set increasingly stringent upper limits on the
amplitude of the 21cm power spectrum. Building on techniques developed in
previous experiments, the Hydrogen Epoch of Reionization Array (HERA; DeBoer
et al. 2017) —currently being built in South Africa—is expected to detect and
characterize the 21cm power spectrum with high significance within the coming
decade. However, while improvements have been made in collecting area and
observing technique, HERA is still expected to be limited by the ability of
the instrument to separate foregrounds and background at a level of 1 part in
$\sim 10^{5}$.

Complementary to the 21cm line is intensity mapping of highly-redshift
Lyman-$\alpha$ (Ly$\alpha$) emission which directly measures high-redshift
galaxies and, indirectly, their effect on the ionized IGM (Silva et al.,
2013). Previous work has suggested that wide-field intensity mapping of the
Ly$\alpha$ line may soon be possible with the future infrared satellites
SPHEREx (Doré et al., 2014) and Cosmic Dawn Intensity Mapper (CDIM; Cooray et
al. 2016).

Cosmological 21cm emission from the IGM is expected to be anti-correlated with
galactic Ly$\alpha$ at large spatial scales enabling study of the interface
between ionized regions dominated by Ly$\alpha$ and neutral regions emitting
21cm. Additionally, due to the fact that radio and infrared foregrounds are
largely uncorrelated, cross-correlating these two measurements may help reduce
problems associated with foreground removal, potentially leading to higher
significance detections of the EoR.

Early studies of cross-correlation between 21cm and Ly$\alpha$ emitters (LAEs)
suggest at its potential for understanding the formation of the first galaxies
(Furlanetto & Lidz, 2007; Lidz et al., 2009). More recent work has shown the
feasibility of cross-correlation between the Square Kilometer Array (SKA) and
the Subaru Hyper Supreme Cam (Hutter et al. 2017; Kubota et al. 2018; Yoshiura
et al. 2018; Weinberger et al. 2020; Kubota et al. 2020), and SKA and CDIM
(Feng et al. 2017; Heneka et al. 2017). In particular, these latter papers
demonstrated through both analytical (Feng et al. 2017) and semi-numerical
modeling (Heneka et al. 2017) that the cross-power spectrum is sensitive to
astrophysical parameters associated with reionization, such as the mean free
path of ionizing photons, and can potentially be used to set tighter
constraints on those parameters than 21cm or Ly$\alpha$ observations could
provide independently. Similarly, Sobacchi et al. 2016 showed that cross-
correlation between 21cm and LAEs improves constraints on the inferred volume-
averaged neutral fraction.

While these results are highly encouraging, many of these estimates use
instruments that may be a decade or more away from construction, delaying a
detection of the 21cm-Ly$\alpha$ until the mid-2030’s. In this paper, we
explore the feasibility of probing the EoR using the 21cm-Ly$\alpha$ cross-
power spectrum by combining HERA observations with the next-generation
infrared probes SPHEREx and CDIM. HERA and SPHEREx are expected to deliver
observations within the next decade. Given their overlapping survey area and
sensitivity of each of these experiments, HERA and SPHEREx may offer the first
opportunity to cross-correlate reionization-era intensity mapping
measurements, which should help tighten constraints on the neutral fraction
and independently confirm each detection. While a cross-power spectrum
measurement made by HERA and SPHEREx should provide some constraints on the
neutral fraction, the total signal-to-noise will likely be low. In order to
fully explore HERA’s potential to cross-correlate with other infrared
instruments, we also provide projections for cross-correlation between HERA
and CDIM. To make this estimate, we use the wide-field sensitivity for CDIM
and imagine a scenario in which a CDIM wide-field survey overlaps with HERA’s
field of view.

This paper is organized as follows. In §2, we discuss the model used for
simulating intensities of line fluctuations, including a model for attenuation
of Ly$\alpha$ by the neutral IGM. We then establish the notation and formally
describe the cross-power spectrum and cross-correlation coefficient in §3. In
§4, we describe the uncertainty associated with a measurement of the cross-
power spectrum, while including a treatment of the 21cm foregrounds, the exact
layout of HERA, and thermal noise contributions from each experiment. Finally,
we discuss the implications of our results and conclude in §5. Throughout this
work, we assume a standard flat $\Lambda$CDM cosmology with the following
parameters: $\Omega_{\Lambda}=0.69$, $\Omega_{\rm m}=0.31$, $\Omega_{\rm
b}=0.049$, $h=0.68$, $n_{s}=0.97$, and $\sigma_{8}=0.81$, which is consistent
with the latest CMB measurements (Planck Collaboration et al., 2020b).

##  2 Modeling Line Intensities

Our goal is to establish the detectability and recoverable information content
given a reasonable simulated prediction of 21cm and Ly$\alpha$. Significant
progress has been made towards accurately modeling the conditions of the
neutral IGM and Ly$\alpha$ emitter population during reionization (McQuinn et
al., 2006; Dijkstra et al., 2007; Santos et al., 2010; Kulkarni et al., 2016).

The major challenge in modeling reionization is accurately capturing the
physical processes that occur at vastly different spatial scales. While
N-body/radiative transfer codes most accurately capture the physical
properties and the evolution of reionization from relatively small scales in
galaxies and large scales in the IGM, they are computationally expensive and
are difficult to extend to the larger cosmological volumes that next-
generation intensity mapping experiments will attempt to observe. Semi-
numerical simulators are much more computationally efficient at modeling the
evolution of the IGM and strongly agree with more numerically motivated
simulators at the large scales probed by intensity mapping experiments within
a few percent. Additionally, computational efficiency is not only more
convenient, but necessary for the efficient exploration of parameter space
required for robust parameter inference.

Here use the semi-numerical code,
21cmFAST555https://github.com/21cmFAST/21cmFAST (Mesinger et al., 2011; Murray
et al., 2020), to simulate 21cm emission and to generate the halo catalogue,
IGM density, and ionization fields necessary for the calculation of Ly$\alpha$
emission in the following subsections.

###  2.1 21cm Brightness Temperature

Figure 1: Slices through the simulated 21cm brightness temperature (top), and
Ly$\alpha$ emission from galaxies (middle) and the ionized IGM (bottom) at
redshifts $z=8.06\left(x_{\rm HI}=0.42\right)$, $z=7.04\left(x_{\rm
HI}=0.05\right)$, and $z=6.0\left(x_{\rm HI}=0.001\right)$. The box size
depicted is 300 Mpc in length with a resolution of 200 cells on each side.
Further details of these simulations are described in §2.1, 2.2.1, and 2.2.2
respectively.

21cmFAST is a semi-numerical simulator of the brightness temperature field as
it evolves with cosmic time. It uses second-order Lagrangian perturbation
theory (2LPT) to evolve a set of initial dark matter density perturbations to
a given redshift. It then uses the excursion set formalism (Bond et al., 1991;
Furlanetto et al., 2004) to identify dark matter haloes in this field, and
prescribes a hydrogen “neutral fraction” $x_{\rm HI}$ per-cell via a series of
theoretically- and empirically-motivated relations. For more details of this
procedure see Mesinger et al. (2011).

The 21cm brightness temperature is observed as a delta on the background CMB

$\displaystyle\delta T_{b}\left(z\right)$
$\displaystyle=\frac{T_{s}-T_{\gamma}}{1+z}\left(1-e^{-\tau_{\nu_{0}}}\right)$
$\displaystyle\approx 27x_{\rm HI}\left(1+\delta_{\rm
nl}\right)\left(\frac{H}{dv_{r}/dv+H}\right)\left(1-\frac{T_{\gamma}}{T_{s}}\right)$
$\displaystyle\hskip
10.00002pt\times\left(\frac{1+z}{10}\frac{0.15}{\Omega_{\rm
M}h^{2}}\right)^{1/2}\left(\frac{\Omega_{\rm b}h^{2}}{0.023}\right)\textrm{
mK.}$ (1)

where $T_{S}$ is the gas spin temperature, $T_{\gamma}$ is the CMB
temperature, $\tau_{\nu_{0}}$ is the optical depth at 21cm frequency,
$\delta_{\rm nl}=\rho/\bar{\rho}_{0}-1$ is the non-linear density contrast,
$H\left(z\right)$ is the Hubble parameter, $dv_{r}/dr$ is the comoving
gradient of the line of sight component of the comoving velocity, where all
quantities are evaluated at redshift $z=\nu_{0}/\nu-1$. The approximation in
the second line of Eq. 1 makes the assumption that reionization has reached
the post-heating regime and that the CMB temperature, $T_{\gamma}$, is much
smaller than the spin gas temperature, $T_{S}$. This approximation allows us
to neglect the full spin gas temperature evolution through reionization when
calculating $\delta T_{b}$, which is much more computationally efficient.
Using the brightness temperature offset, we then calculate the 21cm
fluctuation field, which will be used for cross-correlation in later sections,
with the equation below,

$\delta_{21}\left(\mathbf{x},z\right)=\frac{\delta
T_{b}\left(\mathbf{x},z\right)}{\overline{\delta T}_{b}\left(z\right)}-1.$ (2)

Here $\overline{\delta T}_{b}$ is the spatial average of the 21cm brightness
temperature.

For our simulations we used the fiducial model parameters utilized in 21cmFAST
v3 with a box size of $\left(300\rm Mpc\right)^{3}$. The simulated 21cm
brightness temperature offset field defined in Eq. 1 can be seen in the top
row of Figure 1. Ly$\alpha$ emission from galaxies (§2.2.1) and the ionized
IGM (§2.2.2) can be seen in the following two rows of the same figure for
reference.

###  2.2 Near-Infrared Background

Armed with a 21cm brightness temperature model we turn to generation of a
self-consistent simulated Ly$\alpha$ field which we do by post-processing the
outputs of 21cmFAST. We adopt the technique developed in Silva et al. (2013)
and first applied to 21cmFAST simulations in Heneka et al. (2017) for
simulating Ly$\alpha$ intensity mapping measurements. While those works
already provide in-depth descriptions of the methods used to simulate
Ly$\alpha$ emission, we will restate them here for clarity. In this procedure,
Ly$\alpha$ fluctuations can be modeled from two distinct sources. Their
origins are:

  * •

Ly$\alpha$ Emitters: emission from within the virial radius of dark matter
halos. The dominant components are hydrogen recombinations and collisional
excitation in galaxies.

  * •

Ionized IGM: emission from the bubble of ionized gas that surround
Ly$\alpha$-emitting galaxies. Here, recombinations of ionized hydrogen are the
dominant contributor to Ly$\alpha$ emission.

####  2.2.1 Ly$\alpha$ Emitters

As mentioned above, the emission of Ly$\alpha$ photons in LAE’s is a result of
two dominant processes: the recombination of ionized hydrogen and collisional
excitations of neutral hydrogen within the virial radius of the halos. Silva
et al. (2013) modeled two additional contributions to Ly$\alpha$ emission but
found them to be subdominant, so for this work we just focus on these two main
sources. We predict both types of emission by post-processing the outputs from
21cmFAST which models galaxies and the neutral hydrogen IGM.

Both of these sources of emission are closely related to star formation and
are therefore dependent on the star formation rate (SFR) of the LAEs which
depends in turn on halo mass. We model SFR following Silva et al. (2013) who
extrapolated down from the observed SFR of higher mass ( $M>10^{11}M_{\odot}$)
halos following the empirical model

$\mathrm{SFR}=A\\!\left(\frac{M}{M_{\odot}}\right)^{a}\\!\left(1+\frac{M}{c_{1}}\right)^{b}\\!\left(1+\frac{M}{c_{2}}\right)^{d}M_{\odot}\
{\rm yr}^{-1}$ (3)

where $A=2.8\times 10^{-28}$, $a=2.8$, $b=-0.94$, $c_{1}=10^{9}M_{\odot}$,
$c_{2}=7\times 10^{10}M_{\odot}$, and $d=-1.7$.

The dominant component to Ly$\alpha$ emission in galaxies is the recombination
of ionized hydrogen. As an electron cascades down energy levels during
recombination, it has some probability of emitting a Ly$\alpha$ photon. The
number of Ly$\alpha$ photons being emitted per second through these
recombinations can be estimated using the relationship,

$\dot{N}_{\rm Ly\alpha}\left(M,z\right)=A_{\rm He}f_{\rm rec}f_{\rm
Ly\alpha}\left[1-f_{\rm esc}\left(M,z\right)\right]\dot{N}_{\rm ion}.$ (4)

Here $A_{\rm He}=\left(1-Y_{\rm He}\right)/\left(1-3Y_{\rm He}/4\right)$ for a
helium mass fraction $Y_{\rm He}=0.249$, $f_{\rm rec}$ is the fraction of
hydrogen recombinations that result in a Ly$\alpha$ photon emission, $f_{\rm
Ly\alpha}$ is the fraction of Ly$\alpha$ photons not absorbed by dust, $f_{\rm
esc}\left(M,z\right)$ is the fraction of ionizing photons that escape the halo
dependent on the halo mass following the relationship

$f_{\rm
esc}\left(M,z\right)=\exp\left[-\alpha\left(z\right)M^{\beta\left(z\right)}\right]$
(5)

where $\alpha\left(z\right)$ and $\beta\left(z\right)$ are fitted values as a
function of redshift Table 1, and $\dot{N}_{\rm ion}=Q_{\rm
ion}\times\mathrm{SFR}$ is the rate of ionizing photons emitted by stars. The
value for the average number of photons produced per solar mass of star
formation, $Q_{\rm ion}\approx 5.8\times 10^{60}\ M_{\odot}^{-1}$, was found
by modeling population II stellar lifetimes and estimating the number of
ionizing photons per unit time (Schaerer, 2002).

We estimate a fraction $f_{\rm rec}\approx 66\%$ of the hydrogen
recombinations result in the emission of a Ly$\alpha$ photon by making the
assumption that clouds of interstellar gas are roughly spherical, that the gas
temperature is of the $10^{4}\rm K$, and an ionizing equilibrium in the gas is
present as was calculated in Gould & Weinberg (1996). The Ly$\alpha$ escape
fraction, $f_{\rm Ly\alpha}$, is one that is more challenging to estimate as
the value changes from galaxy to galaxy. For this study, we assume a redshift
parameterization of the Ly$\alpha$ escape fraction that was found in Hayes et
al. (2011),

$f_{\rm Ly\alpha}=C_{\rm dust}\times 10^{-3}\left(1+z\right)^{\xi},$ (6)

where $C_{\rm dust}=3.34$ and $\xi=2.57$. This parameterization was devised
such that $f_{\rm Ly\alpha}$ accounts for the difference between Ly$\alpha$
luminosities found by scaling star formation rates and observed Ly$\alpha$
luminosities assuming that the Ly$\alpha$ photons emitted are a result of
recombinations exclusively. It has been indicated in previous work that the
$f_{\rm Ly\alpha}$ is not only a strong function of redshift, but is dependent
on the halo mass as well, with $f_{\rm Ly\alpha}$ decreasing with increasing
halo mass (Forero-Romero et al., 2011). As has been done in previous studies
(Silva et al., 2013; Heneka et al., 2017; Feng et al., 2017), we stick to the
redshift dependent power law parameterization of the escape fraction as we do
not expect it to the be dominant source of modeling error. We leave further
parameterizations incorporating Ly$\alpha$ escape fraction dependence on halo
mass to future work.

$z$ | $\alpha\left(z\right)$ | $\beta\left(z\right)$ | $f_{\rm esc}\left(M=10^{10}M_{\odot},z\right)$ |
---|---|---|---|---
10.4 | $2.78\times 10^{-2}$ | 0.105 | 0.73 |
8.2 | $1.3\times 10^{-2}$ | 0.179 | 0.45 |
6.7 | $5.18\times 10^{-3}$ | 0.244 | 0.24 |
5.7 | $3.42\times 10^{-3}$ | 0.262 | 0.24 |
Table 1: Escape fraction of UV radiation as a function of redshift fitted in
Razoumov & Sommer-Larsen (2010). For simulation cubes whose redshift falls
between these fitted redshift bins, we interpolate $\alpha\left(z\right)$ and
$\beta\left(z\right)$ to estimate the UV escape fraction as a function of halo
mass.

To calculate the luminosity of Ly$\alpha$ emission from galaxies, we multiply
the rate of Ly$\alpha$ photons being emitted by the energy of Ly$\alpha$
photons to given us an expression for Ly$\alpha$ recombination luminosity of a
galaxy with a given mass, $M$, at redshift, $z$,

$L_{\mathrm{rec}}^{\mathrm{gal}}\left(M,z\right)=E_{\rm Ly\alpha}\dot{N}_{\rm
Ly\alpha}\left(M,z\right),$ (7)

where we assume the emission of Ly$\alpha$ radiation at rest-frame frequency
with energy $E_{\rm Ly\alpha}=13.6\ {\rm eV}$.

The other dominant contributor of Ly$\alpha$ emission in galaxies is
excitation of neutral hydrogen. The Ly$\alpha$ luminosity in the interstellar
medium due to these excitations is defined as,

$L_{\mathrm{exc}}^{\mathrm{gal}}\left(M,z\right)=A_{\rm He}f_{\rm
Ly\alpha}\left[1-f_{\rm esc}\left(M,z\right)\right]E_{\rm exc}\dot{N}_{\rm
ion}\left(M,z\right),$ (8)

where all terms have been previously defined, with the exception of $E_{\rm
exc}\approx 2.14\ {\rm eV}$, which was determined by estimating the average
ionizing photon energy for thermal equilibrium, $E_{\nu}=21.4$ eV (Silva et
al., 2013), and relating it to the energy emitted as Ly$\alpha$ radiation due
to excitations, $E_{\rm exc}/E_{\nu}\approx 0.1$ (Gould & Weinberg, 1996).

With both dominant contributions to Ly$\alpha$ emission modeled, we can then
combine the contributions from excitation and recombination in galaxies to
find a total Ly$\alpha$ luminosity for galaxies dependent on halo mass and
redshift,
$L^{\mathrm{gal}}=L_{\mathrm{exc}}^{\mathrm{gal}}+L_{\mathrm{rec}}^{\mathrm{gal}}$.
Using this expression and the halo catalogue generated by 21cmFAST, we assign
a Ly$\alpha$ luminosity to each one of the halos in the catalogue. We then
create a cube of Ly$\alpha$ luminosity due to the contribution from galaxies
that matches the voxel resolution of our $\delta T_{b}$ simulation cubes by
adding the contribution of each halo to their corresponding voxel position, x.
The result is a Ly$\alpha$ luminosity cube whose emission distribution is
naturally dependent on the spatial distribution of the halos and whose
amplitude is influenced by the mass and clustering of the halos.

This luminosity cube is then converted to a luminosity density by dividing the
galactic contribution to the Ly$\alpha$ luminosity by the comoving voxel size
of the simulation cube, $\ell^{\rm gal}=L^{\rm gal}/V_{\rm vox}$. This can be
converted to a surface brightness using the expression

$I_{\nu}^{\rm
gal}\left(\textbf{x},z\right)=y\left(z\right)D_{A}^{2}\left(z\right)\frac{\ell^{\rm
gal}\left(\textbf{x},z\right)}{4\pi D_{L}^{2}}.$ (9)

Here $D_{A}$ is the comoving angular diameter distance, $D_{L}$ is the
luminosity distance, and the conversion factor from frequency to comoving
distance is $y\left(z\right)=\lambda_{0}\left(1+z\right)^{2}/H\left(z\right)$
(for the rest-frame wavelength of Ly$\alpha$ radiation, $\lambda_{0}=1216$ Å).
Slices of these simulation cubes across the redshift range of interest can be
found in Figure 1.

Figure 2: Ly$\alpha$ power spectra across the redshift range of interest. Both
the LAE and ionized IGM contributions are included in the calculation.

####  2.2.2 Ionized IGM

In this subsection, we describe the Ly$\alpha$ emission from the ionized IGM.
As mentioned at the beginning of this section, we will focus on Ly$\alpha$
emission in the ionized IGM due to hydrogen recombinations. As with Ly$\alpha$
emitting galaxies, the ionized bubbles around galaxies also emit Ly$\alpha$
photons through recombinations of ionized hydrogen. The luminosity density of
Ly$\alpha$ emission in a voxel of the simulation can be defined as

$\ell_{\mathrm{rec}}^{\mathrm{IGM}}\left(\textbf{x},z\right)=n_{\mathrm{rec}}\left(\textbf{x},z\right)f_{\mathrm{rec}}E_{\mathrm{Ly}\alpha},$
(10)

where $E_{\rm Ly\alpha}$ is the rest-frame energy of Ly$\alpha$ photons,
$f_{\rm rec}$ is the fraction of hydrogen recombinations that result in the
emission of a Ly$\alpha$ photon, and $n_{\rm rec}$ is the comoving number
density of recombinations occurring in the ionized IGM. The expression for the
number density of recombinations is

$n_{\mathrm{rec}}\left(\textbf{x},z\right)=\alpha_{\mathrm{A}}n_{e}\left(\textbf{x},z\right)n_{\textsc{HII}}\left(\textbf{x},z\right).$
(11)

Here $\alpha_{\mathrm{A}}$ is the case A comoving recombination coefficient,

$\alpha_{\mathrm{A}}\approx 4.2\times
10^{-13}\left(T_{K}/10^{4}\right)^{-0.7}\left(1+z\right)^{3}{\rm
cm^{3}s^{-1}},$ (12)

$n_{e}\left(\textbf{x},z\right)=x_{i}n_{b}$ is the free electron density, and
$n_{\textsc{HII}}\left(\textbf{x},z\right)=n_{e}A_{\rm He}$ is the comoving
number density of ionized hydrogen. With the luminosity density calculated,
Eq. 9 can be used to calculate the surface brightness of Ly$\alpha$ in the
IGM. The slices through the simulated 21cm and the various components of
Ly$\alpha$ emission cubes are shown in Figure 1. The corresponding power
spectra are found in Figure 2.

In reality, Ly$\alpha$ emission from the ionized IGM also includes the
scattered IGM Ly$\alpha$ background whose main contributors are X-ray and UV
heating, as well as the scattering of Lyman-n photons emitted from galaxies by
residual neutral hydrogen in the ionized IGM (Pritchard & Furlanetto (2007)).
For this work, we chose to neglect the contribution from this Ly$\alpha$
background as its contribution is subdominant to hydrogen recombination in the
ionized IGM (roughly half the mean surface brightness) and the galactic
Ly$\alpha$ contribution (about an order of magnitude lower) (Silva et al.,
2013; Heneka et al., 2017).

####  2.2.3 Ly$\alpha$ Attenuation

Of course the Ly$\alpha$ emission does not propagate unimpeded. We must also
estimate attenuation of Ly$\alpha$ by the neutral IGM. The emitted galactic
Ly$\alpha$ experiences an exponential attenuation

$L^{\textrm{gal}}_{\textrm{obs}}=L^{\textrm{gal}}\exp\left[-\tau_{\textrm{Ly}\alpha}\right],$
(13)

where $L^{\textrm{gal}}$ is the luminosity of Ly$\alpha$ emitting galaxies
defined in the previous section and $\tau_{\textrm{Ly}\alpha}$ is the optical
depth of Ly$\alpha$ emission at some redshift.

The geography of the ionized bubbles around the emitters determines the extent
of the absorption effect. Ly$\alpha$ radiation is emitted by some source and
as that radiation travels to the edge of the ionized bubble, it is redshifted
out of resonance and into the line damping wings, where it has a lower
probability of being absorbed by the neutral medium. As reionization
progresses and the ionized bubbles around halos grow, the probability of
Ly$\alpha$ radiation being attenuated decreases, as UV photons redshift the
further they travel away from the halo before reaching the neutral hydrogen.

To simulate this behavior, we use a model for the optical depth of Ly$\alpha$
emission that is defined in Mesinger & Furlanetto (2008). In this model, the
optical depth of Ly$\alpha$ is related to the neutral fraction of IGM that
photons encounter along the line-of-sight and the amount by which Ly$\alpha$
emission is redshifted as it propagates from galactic halos through the
neutral IGM. Using our simulation cubes, we estimate a value for the optical
depth for Ly$\alpha$ emission, $\tau_{Ly\alpha}$, by tracing skewers from
halos in the simulation through to the edge of the simulation box calculating
the redshift on the near, $z_{bi}$, and far, $z_{ei}$, sides of each neutral
patch and recording the neutral fraction in that particular voxel,
$x_{\textrm{HI}}\left(i\right)$. By tracking these quantities, we calculate
the contribution from each neutral hydrogen patch encountered along the line-
of-sight of the halo using the approximation (Miralda-Escudé, 1998)

$\displaystyle\tau_{\textrm{Ly}\alpha}$
$\displaystyle=\tau_{\textrm{s}}\sum_{i}x_{\textrm{HI}}\left(i\right)\left(\frac{2.02\times
10^{-8}}{\pi}\right)\left(\frac{1+z_{bi}}{1+z_{\textrm{s}}}\right)^{3/2}$
$\displaystyle\hskip
10.00002pt\times\left[I\left(\frac{1+z_{bi}}{1+z_{\textrm{s}}}\right)-I\left(\frac{1+z_{ei}}{1+z_{\textrm{s}}}\right)\right].$
(14)

where $z_{s}$ is the redshift of the Ly$\alpha$ emitting source. This assumes
that the optical line depth at Ly$\alpha$ line resonance, $\tau_{s}$, can be
approximated as

$\tau_{\textrm{s}}\approx 6.45\times
10^{5}\left(\frac{\Omega_{b}h}{0.03}\right)\left(\frac{\Omega_{m}}{0.3}\right)^{-0.5}\left(\frac{1+z_{\textrm{s}}}{10}\right)^{3/2}$
(15)

at high-redshifts for a source at some redshift, $z_{s}$, given present-day
$\Omega_{b}$, and $\Omega_{m}$ (Gunn & Peterson 1965; Barkana & Loeb 2001). In
the expression above, $I\left(z\right)$ is the helper function,

$I\left(x\right)=\frac{x^{4.5}}{1-x}+\frac{9}{7}x^{3.5}+\frac{9}{5}x^{2.5}+3x^{1.5}+9x^{0.5}-\ln\left(\frac{1+x^{0.5}}{1-x^{0.5}}\right),$
(16)

derived in Miralda-Escudé (1998). This expression is only valid for
frequencies far from the line center, but is used in this case because the
optical depth is so large at the line center that the emission becomes
attenuated to the point of being unobservable and therefore gives a fairly
accurate approximation for the optical depth at high redshifts.

Once calculated for each halo in the simulation cube, $\tau_{\rm Ly\alpha}$ is
applied as an attenuation factor to each halo’s intrinsic luminosity through
Eq. 13 before constructing the cubes and calculating the Ly$\alpha$ power
spectrum and 21cm-Ly$\alpha$ cross-power spectrum for sensitivity calculations
in Section 4.

##  3 Cross-Correlation Statistics

As direct observation of 21cm image cubes will require next generation
sensitivity and precision, current experiments are opting to focus on the
power spectrum, a measurement of statistical fluctuations on various spatial
scales. The power spectrum, $P\left({\bf k}\right)$, is formally defined as

$\langle\widetilde{\delta T}_{b}({\bf k})\widetilde{\nu
I}_{\mathrm{Ly}\alpha}({\bf
k^{\prime}})\rangle\equiv(2\pi)^{3}\delta^{D}\left({\bf k}-{\bf
k^{\prime}}\right)P_{21,\mathrm{Ly}\alpha}\left({\bf k}\right),$ (17)

where $\widetilde{\delta}$ is the Fourier transform of some fluctuation field,
in this case either 21cm or Ly$\alpha$, $\delta^{D}$ is the Dirac delta
function. More commonly used in the literature is the dimensionless power
spectrum,

$\Delta^{2}_{21,\mathrm{Ly}\alpha}\left(k\right)=\frac{k^{3}}{2\pi^{2}}P_{21,\mathrm{Ly}\alpha}\left(k\right),$
(18)

which represents the contribution to the co-variance in bins of $k$, and in
which the $k^{3}$ cancels the spatial units from the power spectrum. Figure 3
shows the cross-power spectrum for a range of redshifts from $z\approx 6-8$.

Figure 3: The dimensionless 21cm-Ly$\alpha$ cross-power spectrum. The solid
lines on each of the curves represent positive values in the cross-power
spectrum, while dotted line on the same curves represent negative values. The
cross-power spectrum is expected to turn over from positive to negative on the
scale of the mean ionized bubble size at that redshift. We find the cross-
power spectrum turns over at increasingly large spatial scales (small
$k$-modes) as reionization progresses, tracing the growth of ionized bubbles.

In the power spectrum it is difficult to decouple the total power of the two
fields from their correlation. Though we can’t directly measure the
correlation it is useful to inspect it in simulation as a way to understand
the dependence of the correlation on size scale. We expect the correlation
coefficient

$r_{21,\rm Ly\alpha}\left(k\right)=\frac{P_{21,\rm
Ly\alpha}\left(k\right)}{\sqrt{P_{21}\left(k\right)P_{\rm
Ly\alpha}\left(k\right)}},$ (19)

to have an anti-correlation $(r<0)$ on large scales dominated by bubbles which
approaches zero at a scale roughly corresponding to the average size of
ionized regions.

Figure 4: The cross-correlation coefficient as a function of spatial scale.
Here, a cross-correlation coefficient value of $-1$ indicates that 21cm and
Ly$\alpha$ emission are completely anti-correlated, while a value of 0
indicates no correlation. At a fixed spatial scale, $r\left(k\right)$ becomes
less negative between $z\approx 8$ and $z\approx 6$ as reionization
progresses, which tracks the growth of ionized bubbles around galaxies. Above
redshift of $z\sim 8$ we see an increase in the correlation at large scales
which is explained by few galaxies forming in overdense regions that haven’t
yet been ionized leading to a correlation between 21cm and Ly$\alpha$
fields(Lidz et al., 2009).

The correlation (shown in Figure 4) behaves as expected, transitioning from
uncorrelated on small-scales to anti-correlated on large-scales, where the
fluctuations from the two fields are not overlapping. It is also interesting
to note that the cross-correlation coefficient progresses from generally anti-
correlated at high-redshifts to generally uncorrelated towards the end of
reionization. This change traces the growth of ionized bubbles through the
process of reionization. The scales at which the cross-correlation coefficient
transitions from uncorrelated to anti-correlated represents the size of
typical ionized bubbles around galaxies, which matches results seen in
previous literature (Lidz et al., 2009; Heneka et al., 2017; Kubota et al.,
2018).

##  4 Detectability of the Cross-Power Spectrum

###  4.1 Thermal Noise Contribution

The feasibility of a 21cm-Ly$\alpha$ cross-power spectrum measurement is
governed by sources of uncertainty associated with both measurements
independently and their overlap. Calculation of thermal error bars captures
both. In this work, we use HERA as our 21cm instrument, and SPHEREx and CDIM
individually for Ly$\alpha$.

The noise on one particular $k$-mode in the cross-power spectrum depends on
noise from both the Ly$\alpha$ and 21cm measurements, and therefore the
uncertainty from both measurements must be calculated. The variance on the
cross-power spectrum due to contributions from both instruments and sample
variance from the cross power spectrum is,

$\displaystyle\sigma^{2}$ ${}_{21,\rm
Ly\alpha}\left(\textbf{k}\right)=\frac{1}{2}\left[P^{2}_{21,\rm
Ly\alpha}\left(\textbf{k}\right)+\sigma_{21}\left(\textbf{k}\right)\sigma_{\rm
Ly\alpha}\left(\textbf{k}\right)\right]$ $\displaystyle\propto P^{2}_{21,\rm
Ly\alpha}+\left(P_{21}+P_{21,\mathrm{N}}\right)\left(P_{\rm Ly\alpha}+P_{\rm
Ly\alpha,\mathrm{N}}\right).$ (20)

Here, $P_{21,\rm Ly\alpha}$ is sample variance due to the cross-power
spectrum, $P_{21}$ is sample variance from the 21cm signal, $P_{\rm Ly\alpha}$
is sample variance from the $\rm Ly\alpha$ measurement, and $P_{21,N}$ and
$P_{{\rm Ly\alpha},N}$ are thermal noise uncertainty terms from 21cm and
$Ly\alpha$ measurements respectively.

To calculate the thermal noise for the 21cm observation, we use the method
described by Pober et al. (2013a)666We used the updated version of the code at
https://github.com/steven-murray/21cmSense. This version adds support for
modern configuration formats like YAML, significant unit-testing,
documentation and tutorials, and greatly modularises the code.. In this
method, uv-coverage of the observation is taken into account using the exact
layout of HERA and applying Earth-rotation synthesis to simulate changing uv-
bins sampled by each pair of antennas. This uv-coverage then dictates the
exact $k_{\perp}$ resolution of the instrument while its spectral resolution
sets the $k_{\parallel}$ resolution. We then perform a spherical average over
the $k_{\parallel}$ and $k_{\perp}$ bins to identify the observation time,
$t_{\rm int}$, associated with each $k$-bin. As described in Pober et al.
(2014), the power spectrum error depends on the observation time
$t_{\textrm{int}}$ and the system temperature $T_{\rm sys}$ according to

$P_{21,\textrm{N}}=X^{2}Y\frac{T^{2}_{\textrm{sys}}}{2t_{\mathrm{int}}}\frac{\Omega_{\textrm{p}}^{2}}{\Omega_{\textrm{pp}}}.$
(21)

Here, $X^{2}Y$ converts bandwidth and solid angle to their relevant
cosmological scale equivalents. The field of view also factors in via
$\Omega_{p}$, the solid angle of the primary beam, and $\Omega_{\rm pp}$, the
solid angle of the primary beam squared. For the thermal noise estimates
above, the fiducial observing parameters for HERA ($t_{\rm int}=1000\ \rm hr$,
$T_{\rm sys}=100+120\left(\nu/150\right)^{-2.55}$, and $B=8\ \rm MHz$) are
assumed (DeBoer et al., 2017).

The thermal noise associated with infrared intensity mapping experiments
(SPHEREx and CDIM in this work) can be written as

$P_{\mathrm{Ly}\alpha,\mathrm{N}}=\left(\nu_{\rm obs}\sigma_{\rm
N}\right)^{2}V_{\mathrm{vox}}\textrm{W}_{\rm Ly\alpha},$ (22)

where $\nu_{\rm obs}$ is the observed frequency of Ly$\alpha$,
$V_{\mathrm{vox}}$ is the comoving voxel volume and $\textrm{W}_{\rm
Ly\alpha}$ is the window function defined in Lidz et al. (2011), which
accounts for limitations in the spectral and spatial instrumental resolution
of the instrument. In the equation above, we take $\sigma_{\rm N}=10^{-18}\
{\rm erg\ s^{-1}\ cm^{-2}\ sr^{-1}\ Hz}$, which is consistent with the
5$\sigma$ value reported by Doré et al. (2016) for SPHEREx, and $\sigma_{\rm
N}=1.5\times 10^{-19}\ {\rm erg\ s^{-1}\ cm^{-2}\ sr^{-1}\ Hz}$ for CDIM
Cooray et al. (2019). In addition to these values, we also explore more
optimistic sensitivity values assuming deeper integrations than the minimum
system requirements for both SPHEREx and CDIM as an exploratory measure for
the requirements of a detection of the cross-power spectrum. For these
optimistic values, we take $\sigma_{\rm N}=3\times 10^{-20}\ {\rm erg\ s^{-1}\
cm^{-2}\ sr^{-1}\ Hz^{-1}}$777https://github.com/SPHEREx/Public-products/ for
SPHEREx, and $\sigma_{\rm N}=1.5\times 10^{-21}\ {\rm erg\ s^{-1}\ cm^{-2}\
sr^{-1}\ Hz^{-1}}$ for CDIM Heneka et al. (2017). We do expect the nominal
values to be the more realistic values measured, but better thermal noise
values are possible given certain instrument assumptions (Doré et al., 2018;
Symons et al., 2021).

Eq. 20 is defined for the unaveraged power spectra. In practice, we are more
interested in the spherically averaged noise power spectrum. To obtain the
variance on the spherically averaged cross-power spectra, we use

$\frac{1}{\sigma^{2}\left(k\right)}=\sum_{\textbf{k}\in
k}\frac{N_{m}}{\sigma^{2}\left(\textbf{k}\right)},$ (23)

where $N_{m}$ is the number of modes within a particular $k$-bin, which are
explicitly counted when averaging. The results of these cuts can be found in
Figure 5. A discussion of the results of this figure coupled with the
different foregrounds strategies in the following section.

###  4.2 Foreground Contamination

In addition to the thermal noise and limited spectral and spatial resolution
of each instrument, we also like to explore the effect of foregrounds on our
ability to effectively measure the cross-power spectrum. Cross-correlation of
reionization-era 21cm observations with Ly$\alpha$ intensity mapping surveys
have the advantage that 21cm foregrounds are expected to have no correlation
with low redshift interloper lines that affect high-redshift Ly$\alpha$
measurements. Because the two are uncorrelated, power from the each of these
foregrounds should not be directly added to the cross-power spectrum. However,
while bright foregrounds are not expected to contribute to the amplitude of
the cross-power spectrum, they do contribute to the total variance on the
cross-power spectrum if not removed. To truly be confident in a detection of
the cross-power spectrum, foregrounds must be accounted for.

While Ly$\alpha$ intensity mapping experiments do have to contend with
interloper lines at lower redshifts, such as H$\alpha$, Oii, and Oiii that are
orders of magnitude brighter than the infrared background, much work has gone
into effectively removing these lines. The technique for removing foreground
sources simply involves applying a cutoff flux and removing all pixels whose
amplitude falls above that threshold. Previous work has shown that for a
SPHEREx-like infrared satellite with spectral resolution $R\approx 41.5$ only
3% of pixels would need to be removed to bypass all foreground interlopers
(Pullen et al., 2014; Gong et al., 2014). Similarly, Feng et al. (2017) found
that for a CDIM-like experiment $R\approx 300$ only 0.1% of pixels would need
to be removed to significantly lower the amplitude of the interloper to orders
of magnitude below the Ly$\alpha$ power spectrum. For this reason, in this
work we ignore the effect infrared foreground removal would have on decreasing
the amplitude of the cross-power spectrum and instead focus on the effect of
21cm foreground removal on the cross-power spectrum.

Figure 5: Error budget of the sensitivity of the cross-power spectrum at
$z=7.04$ between the HERA 21cm experiment and Ly$\alpha$ observations by
SPHEREx and CDIM for the deep integration case. Here, the cross terms defined
in Eq. 20 are plotted individually to isolate the sources of error on a
measurement of cross-power spectrum. The relevant sensitivity threshold is the
quadradic sum of all terms which is dominated by the whichever term is
largest. The simulated cross-power spectrum is shown in grey with negative
values in dashes. Moderate 21cm foreground filtering in shown at top and
optimistic in the bottom row. See section 4.2 for a complete discussion.

The foreground wedge is a well-documented feature in 21cm literature (eg.
Datta et al., 2010; Morales et al., 2012; Liu et al., 2014) and filtering
matched to this shape has been identified as a potential method of removing
bright 21cm foregrounds. The foreground wedge is a feature that appears in the
cylindrically-averaged 21cm power spectrum, arising as a product of spectrally
smooth foregrounds and the chromatic response of the interferometer. Because
of their smooth spectral structure, bright foregrounds are confined to low-
order Fourier modes, thus constraining their power to low $k_{\parallel}$
values. Experiments such as HERA have leveraged the fact that the edge of the
foreground wedge is dependent on the baseline length of two antennas by
building densely packed arrays that sample lower $k_{\perp}$ values, thus
increasing the EoR window. The approximate relationship defining the
$k_{\parallel}$-edge of the wedge can be written as

$k_{\parallel,\mathrm{max}}=\frac{D_{M}}{D_{H}}\frac{E\left(z\right)}{\left(1+z\right)}\sin{\theta_{0}}k_{\perp}$
(24)

where the characteristic angle, $\theta_{0}$, is the instrument field of view,
$k_{\perp}$ is the Fourier mode dependent on the distance between two dishes,
and $k_{\parallel,\mathrm{max}}$ corresponds to the maximum $k_{\parallel}$
value dominated by bright foregrounds.

Typically, the safest assumption to make is that $\theta_{0}=\pi/2$, which
corresponds to a field of view that includes bright foregrounds at the
horizon, far from the pointing center. In practice, the wedge can extend even
beyond the horizon given imperfectly calibrated chromaticity and internal
instrument systematics, such as cable reflections and cross-coupling effects
(Pober et al., 2013b; Kern et al., 2020). However, work by (Pober et al.,
2014) has argued that outside the field of view, foreground contamination is
sufficiently attenuated by the primary beam so as not to corrupt the
cosmological 21cm signal. If this is achievable in practice, it would increase
the size of the EoR window and provide a significant sensitivity boost to a
cross-power spectrum measurement.

To investigate the effect of the foreground wedge on the ability to measure
the cross-power spectrum, we adopt two treatments of 21cm foregrounds
described in Pober et al. (2014): a moderate foreground treatment, where the
foreground wedge extends to the horizon with a horizon buffer added to account
for improper calibration, and an optimistic foreground treatment, where the
wedge is confined to the first null in the primary beam of the instrument. In
both of these treatments, we remove all $\mathbf{k}$-modes that fall within
the foreground wedge of the cylindrically-averaged power spectra before
averaging down to the spherically-averaged power spectra. We compute each of
the cross terms defined in Eq, 20 shown in Figure 5. There are two signals for
four cross terms. The total error is the quadrature sum of these terms which
is well approximated by the largest error on the plot. For SPHEREx the two
largest components are the $\rm Ly\alpha,N$ noise-noise term followed closely
by the $\rm 21cm-Ly\alpha$ cross term. This suggests the sensitivity of the
Ly$\alpha$ measurement is the limiting factor. For CDIM, however the reverse
is true. At the most sensitive $k$-modes the dominant source of error is the
correlation between $\rm Ly\alpha$ sample variance and 21cm thermal noise
uncertainty, suggesting that with CDIM’s greater range of available modes the
21cm sensitivity becomes the limiting factor.

The other strong factor at work is the 21cm foreground treatment; filtering
modes measured by both instruments reduces sensitivity with a strong
dependence on $k$. In the moderate treatment there is a sharp rise in the
uncertainty towards lower $k$s. This is due to the complete loss of large
scales in the 21cm wedge filter. We can also see that sensitivity improves
across the board when the filter is relaxed reflecting the fact that a range
of $k_{\perp}$ are being included. For SPHEREx, the filter has the largest
impact on the cross term between the 21cm signal and $\rm Ly\alpha$ noise.
With the limited number of modes overlapping between HERA and SPHEREx, the
loss of modes to the foreground filter is keenly felt. With its much wider
range of available modes wedge filtering has a much smaller impact on the
correlation with CDIM.

###  4.3 Sensitivity Estimates

With the thermal noise, spectral and spatial resolution effects, and
foregrounds taken into account, we can examine the sensitivity these
instrument pairs have to the cross-power spectrum. Using the cross-power
spectrum and the noise cross-power spectrum calculated previously, we can
calculate the total signal-to-noise ratio across the cross-power spectrum by
summing across the $k$-bins using the expression

$\mathrm{SNR^{2}_{total}}=\sum_{i}\mathrm{SNR}^{2}_{i}=\sum_{i}\left(\frac{P_{21,\mathrm{Ly}\alpha}\left(k_{i}\right)}{\sigma_{21,\mathrm{Ly}\alpha}\left(k_{i}\right)}\right)^{2},$
(25)

where the index, $i$, iterates through each of the $k$-bins. This signal-to-
noise calculation was done for each redshift bin, for both instruments, and
both foreground treatments. These calculated signal-to-noise ratios for
HERA/SPHEREx and HERA/CDIM cross-power spectra are found in Figure 6.

Figure 6: Integrated signal-to-noise on the cross-power spectrum as a function
of redshift. For each panel, the blue line represents SNR estimates on the
cross-power spectrum in the moderate foreground case, while the purple line
represents the optimistic foreground case. In addition to different foreground
strategies, we also differentiate between the minimum thermal noise
requirement for each instrument (solid lines) and a more optimal deep
integration (dotted-dashed lines) for each Ly$\alpha$ experiment

The overall SNR prediction (Figure 6) tells us that cross-correlation with
SPHEREx requires an optimistic treatment of the 21cm foregrounds and deeper
integrations for SPHEREx than the minimum requirements to make a detection of
the cross-power spectrum; which here means making noise-limited measurements
at delay bins up to the first null of the beam. SPHEREx and HERA sample
different $\mathbf{k}$-mode ranges which overlap best at small $k$. Meanwhile,
the correlation power spectrum, like the auto spectrum, is roughly flat in
$k$. Looking further into the future, a HERA/CDIM-like cross-correlation may
be possible even if the entire wedge is excluded from the 21cm data. This is
made possible by the much larger sensitivity of CDIM at 21cm $k$-modes which
are foreground free.

##  5 Summary

We have tested the feasibility of detecting large scale structure during
reionization by cross-correlating HERA with two future infrared intensity
mapping satellites, SPHEREx and CDIM. In the near future, using cross-
correlations between HERA and SPHEREx, we find that the cross-power spectrum
may be detectable from $z\approx 7-8.5$ in the fiducial 21cm model, but only
with aggressive removal of 21cm foregrounds and deeper thermal noise
integrations than the minimum system requirements. This is due to a lack of
overlapping sensitive modes between the two instruments that is only remedied
by aggressive foreground removal or deep integrations from either instrument.
A HERAxSPHEREx cross-correlation will likely set upper limits on the intensity
of the cross-power spectrum, as well as some constraints on astrophysical
parameters, but a detection of the cross-power spectrum will be challenging
without better control of systematics and improved foreground subtraction.
More forward-looking, we also show that a HERAxCDIM cross-power spectrum
measurement ought to have sufficient signal-to-noise to detect the fiducial
model across a significant portion of reionization ($z\approx
6-9;x_{\textrm{HI}}\approx 0.01-0.75$), assuming a fairly aggressive
foreground treatment and slightly deeper integrations than the minimum system
requirements.

More work is needed to extend this initial study. We treated foregrounds by
removing likely contaminated modes. One could imagine that small residual
foregrounds could cancel in cross-correlation, thus opening up more modes.
While it is likely that residual foregrounds from even the best levels of
foreground modeling and subtraction currently in use would dominate the error
budget of the cross-power spectrum, the relatively low noise contribution from
CDIM may eventually allow the cross-power spectrum to be detectable even in
cases of imperfect 21cm foreground removal at large scales.

Potentially the most exciting aspect of cross-correlating 21cm and Ly$\alpha$
measurements will be determining its ability to constrain astrophysical and
cosmological model parameters using the cross-power spectrum. Synergies
between line intensity mapping experiments will help drive constraints on
parameter estimates by breaking down expected degeneracies. We leave a study
of astrophysical and cosmological parameter estimation for future work.

## Acknowledgements

The authors would like to thank Alexander van Engelen and Judd Bowman for
their helpful comments on an early draft. This material is based upon work
supported by the National Science Foundation under Grant #1636646 and #1836019
and institutional support from the HERA collaboration partners. This research
is funded in part by the Gordon and Betty Moore Foundation. HERA is hosted by
the South African Radio Astronomy Observatory, which is a facility of the
National Research Foundation, an agency of the Department of Science and
Innovation.

### Software

This work was enabled by a number of software packages including 21cmSense
(Pober et al., 2013a), astropy (Astropy Collaboration et al., 2013),
matplotlib (Hunter, 2007), numpy (Oliphant, 2006), powerbox (Murray, 2018),
and scipy (Virtanen et al., 2020).

### Data Availability

The 21cm and Ly$\alpha$ simulation cubes and noise curves used in this study
can be provided upon reasonable request to the main author.

## References

  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Barkana & Loeb (2001) Barkana, R., & Loeb, A. 2001, Phys. Rep., 349, 125, doi: 10.1016/S0370-1573(01)00019-9
  * Bond et al. (1991) Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440, doi: 10.1086/170520
  * Cooray et al. (2016) Cooray, A., Bock, J., Burgarella, D., et al. 2016, arXiv e-prints, arXiv:1602.05178.  https://arxiv.org/abs/1602.05178
  * Cooray et al. (2019) Cooray, A., Chang, T.-C., Unwin, S., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 23.  https://arxiv.org/abs/1903.03144
  * Datta et al. (2010) Datta, A., Bowman, J. D., & Carilli, C. L. 2010, ApJ, 724, 526, doi: 10.1088/0004-637X/724/1/526
  * DeBoer et al. (2017) DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001, doi: 10.1088/1538-3873/129/974/045001
  * Dijkstra et al. (2007) Dijkstra, M., Lidz, A., & Wyithe, J. S. B. 2007, MNRAS, 377, 1175, doi: 10.1111/j.1365-2966.2007.11666.x
  * Doré et al. (2014) Doré, O., Bock, J., Ashby, M., et al. 2014, arXiv e-prints, arXiv:1412.4872.  https://arxiv.org/abs/1412.4872
  * Doré et al. (2016) Doré, O., Werner, M. W., Ashby, M., et al. 2016, arXiv e-prints, arXiv:1606.07039.  https://arxiv.org/abs/1606.07039
  * Doré et al. (2018) Doré, O., Werner, M. W., Ashby, M. L. N., et al. 2018, arXiv e-prints, arXiv:1805.05489.  https://arxiv.org/abs/1805.05489
  * Fan et al. (2006) Fan, X., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415, doi: 10.1146/annurev.astro.44.051905.092514
  * Feng et al. (2017) Feng, C., Cooray, A., & Keating, B. 2017, ApJ, 846, 21, doi: 10.3847/1538-4357/aa7ff1
  * Forero-Romero et al. (2011) Forero-Romero, J. E., Yepes, G., Gottlöber, S., et al. 2011, MNRAS, 415, 3666, doi: 10.1111/j.1365-2966.2011.18983.x
  * Furlanetto & Lidz (2007) Furlanetto, S. R., & Lidz, A. 2007, ApJ, 660, 1030, doi: 10.1086/513009
  * Furlanetto et al. (2006) Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, Phys. Rep., 433, 181, doi: 10.1016/j.physrep.2006.08.002
  * Furlanetto et al. (2004) Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 1, doi: 10.1086/423025
  * Gehlot et al. (2019) Gehlot, B. K., Mertens, F. G., Koopmans, L. V. E., et al. 2019, MNRAS, 488, 4271, doi: 10.1093/mnras/stz1937
  * Gong et al. (2014) Gong, Y., Silva, M., Cooray, A., & Santos, M. G. 2014, ApJ, 785, 72, doi: 10.1088/0004-637X/785/1/72
  * Gould & Weinberg (1996) Gould, A., & Weinberg, D. H. 1996, ApJ, 468, 462, doi: 10.1086/177707
  * Gunn & Peterson (1965) Gunn, J. E., & Peterson, B. A. 1965, ApJ, 142, 1633, doi: 10.1086/148444
  * Hayes et al. (2011) Hayes, M., Schaerer, D., Östlin, G., et al. 2011, ApJ, 730, 8, doi: 10.1088/0004-637X/730/1/8
  * Heneka et al. (2017) Heneka, C., Cooray, A., & Feng, C. 2017, ApJ, 848, 52, doi: 10.3847/1538-4357/aa8eed
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Hutter et al. (2017) Hutter, A., Dayal, P., Müller, V., & Trott, C. M. 2017, ApJ, 836, 176, doi: 10.3847/1538-4357/836/2/176
  * Kern et al. (2020) Kern, N. S., Parsons, A. R., Dillon, J. S., et al. 2020, ApJ, 888, 70, doi: 10.3847/1538-4357/ab5e8a
  * Kolopanis et al. (2019) Kolopanis, M., Jacobs, D. C., Cheng, C., et al. 2019, ApJ, 883, 133, doi: 10.3847/1538-4357/ab3e3a
  * Kubota et al. (2020) Kubota, K., Inoue, A. K., Hasegawa, K., & Takahashi, K. 2020, MNRAS, 494, 3131, doi: 10.1093/mnras/staa979
  * Kubota et al. (2018) Kubota, K., Yoshiura, S., Takahashi, K., et al. 2018, MNRAS, 479, 2754, doi: 10.1093/mnras/sty1471
  * Kulkarni et al. (2016) Kulkarni, G., Choudhury, T. R., Puchwein, E., & Haehnelt, M. G. 2016, MNRAS, 463, 2583, doi: 10.1093/mnras/stw2168
  * Li et al. (2019) Li, W., Pober, J. C., Barry, N., et al. 2019, ApJ, 887, 141, doi: 10.3847/1538-4357/ab55e4
  * Lidz et al. (2011) Lidz, A., Furlanetto, S. R., Oh, S. P., et al. 2011, ApJ, 741, 70, doi: 10.1088/0004-637X/741/2/70
  * Lidz et al. (2009) Lidz, A., Zahn, O., Furlanetto, S. R., et al. 2009, ApJ, 690, 252, doi: 10.1088/0004-637X/690/1/252
  * Liu et al. (2014) Liu, A., Parsons, A. R., & Trott, C. M. 2014, Phys. Rev. D, 90, 023018, doi: 10.1103/PhysRevD.90.023018
  * McQuinn et al. (2006) McQuinn, M., Zahn, O., Zaldarriaga, M., Hernquist, L., & Furlanetto, S. R. 2006, ApJ, 653, 815, doi: 10.1086/505167
  * Mesinger et al. (2011) Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955, doi: 10.1111/j.1365-2966.2010.17731.x
  * Mesinger & Furlanetto (2008) Mesinger, A., & Furlanetto, S. R. 2008, MNRAS, 385, 1348, doi: 10.1111/j.1365-2966.2007.12836.x
  * Miralda-Escudé (1998) Miralda-Escudé, J. 1998, ApJ, 501, 15, doi: 10.1086/305799
  * Morales et al. (2012) Morales, M. F., Hazelton, B., Sullivan, I., & Beardsley, A. 2012, ApJ, 752, 137, doi: 10.1088/0004-637X/752/2/137
  * Murray (2018) Murray, S. G. 2018, The Journal of Open Source Software, 3, 850, doi: 10.21105/joss.00850
  * Murray et al. (2020) Murray, S. G., Greig, B., Mesinger, A., et al. 2020, Journal of Open Source Software, 5, 2582, doi: 10.21105/joss.02582
  * Oliphant (2006) Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (Trelgol Publishing USA)
  * Paciga et al. (2013) Paciga, G., Albert, J. G., Bandura, K., et al. 2013, MNRAS, 433, 639, doi: 10.1093/mnras/stt753
  * Planck Collaboration et al. (2020a) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020a, A&A, 641, A5, doi: 10.1051/0004-6361/201936386
  * Planck Collaboration et al. (2020b) —. 2020b, A&A, 641, A6, doi: 10.1051/0004-6361/201833910
  * Pober et al. (2013a) Pober, J. C., Parsons, A. R., DeBoer, D. R., et al. 2013a, AJ, 145, 65, doi: 10.1088/0004-6256/145/3/65
  * Pober et al. (2013b) Pober, J. C., Parsons, A. R., Aguirre, J. E., et al. 2013b, ApJ, 768, L36, doi: 10.1088/2041-8205/768/2/L36
  * Pober et al. (2014) Pober, J. C., Liu, A., Dillon, J. S., et al. 2014, ApJ, 782, 66, doi: 10.1088/0004-637X/782/2/66
  * Pritchard & Furlanetto (2007) Pritchard, J. R., & Furlanetto, S. R. 2007, MNRAS, 376, 1680, doi: 10.1111/j.1365-2966.2007.11519.x
  * Pullen et al. (2014) Pullen, A. R., Doré, O., & Bock, J. 2014, ApJ, 786, 111, doi: 10.1088/0004-637X/786/2/111
  * Razoumov & Sommer-Larsen (2010) Razoumov, A. O., & Sommer-Larsen, J. 2010, ApJ, 710, 1239, doi: 10.1088/0004-637X/710/2/1239
  * Santos et al. (2010) Santos, M., Ferramacho, L., Silva, M., Amblard, A., & Cooray, A. 2010, SimFast21: Simulation of the Cosmological 21cm Signal.  http://ascl.net/1010.025
  * Schaerer (2002) Schaerer, D. 2002, A&A, 382, 28, doi: 10.1051/0004-6361:20011619
  * Silva et al. (2013) Silva, M. B., Santos, M. G., Gong, Y., Cooray, A., & Bock, J. 2013, ApJ, 763, 132, doi: 10.1088/0004-637X/763/2/132
  * Sobacchi et al. (2016) Sobacchi, E., Mesinger, A., & Greig, B. 2016, MNRAS, 459, 2741, doi: 10.1093/mnras/stw811
  * Symons et al. (2021) Symons, T., Zemcov, M., Bock, J., et al. 2021, ApJS, 252, 24, doi: 10.3847/1538-4365/abcaa5
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Weinberger et al. (2020) Weinberger, L. H., Kulkarni, G., & Haehnelt, M. G. 2020, MNRAS, 494, 703, doi: 10.1093/mnras/staa749
  * Yoshiura et al. (2018) Yoshiura, S., Line, J. L. B., Kubota, K., Hasegawa, K., & Takahashi, K. 2018, MNRAS, 479, 2767, doi: 10.1093/mnras/sty1472