# Calcium bright knots and the formation of chromospheric anemone jets on the
Sun

Kunwar Alkendra Pratap Singh  Plasma Astrophysics Research Laboratory,
Department of Physics, Institute of Science, BHU, Varanasi-221005, India
Astronomical Observatory, Graduate School of Science, Kyoto University,
Yamashina, Kyoto 607-8471, Japan  Keisuke Nishida Astronomical Observatory,
Graduate School of Science, Kyoto University, Yamashina, Kyoto 607-8471, Japan
Kazunari Shibata Astronomical Observatory, Graduate School of Science, Kyoto
University, Yamashina, Kyoto 607-8471, Japan  School of Science and
Engineering, Doshisha University, Kyotanabe-city, Kyoto 610-0394, Japan

###### Abstract

Space-based observations show that the solar atmosphere from the solar
chromosphere to the solar corona is filled with small-scale jets and is linked
with small-scale explosions. These jets may be produced by mechanisms similar
to that of large-scale flares and such jets may be related to the heating of
corona and chromosphere as well as the acceleration of solar wind. The
chromospheric anemone jets on the Sun remain puzzling because their footpoints
(or bright knots) have not been well resolved and the formation process of
such enigmatic small-scale jets remains unclear. We propose a new model for
chromospheric jets using the three-dimensional magnetohydrodynamic (MHD)
simulations, which show that the continuous, upward rising of small-scale
twisted magnetic flux ropes in a magnetized solar chromosphere drive small-
scale magnetic reconnection and the launching of several small-scale jets
during the evolution of the chromospheric anemone jets. Our new, self-
consistent, three-dimensional computer modeling of small-scale, but ever-
changing flux rope emergence in the magnetized solar atmosphere is fully
consistent with observations and provides a universal mechanism for nanoflare
and jet formation.

Solar chromosphere (1479); Solar activity (1475); Solar magnetic fields
(1503); Solar magnetic reconnection (1504)

##  1 Introduction

The small-scale jets, associated with small-scale explosions fill the solar
atmosphere present between the solar chromosphere and to solar corona
(Schmieder 2022). There is increasing indirect evidence that these jets may be
produced by a physical mechanism similar to that of large-scale flares, i.e.,
the release of magnetic energy stored in the solar atmosphere through magnetic
reconnection (e.g. Yokoyama & Shibata 1995, Shibata & Magara 2011, Sterling
2015), and also that these jets may be related to the heating of corona and
chromosphere as well as the acceleration of solar wind (Raouafi 2016, Sterling
& Moore 2020). Among small-scale jets, the so-called chromospheric anemone
jets has a inverse Y-shape structure in active regions on the Sun remain
puzzling, because they are at such small scales that their footpoints (bright
knots) have not been well resolved (Shibata et al. 2007, Nishizuka et al.
2011, Singh et al. 2012a,b).

In addition to the chromospheric anemone jets, the solar chromosphere harbors
a variety of jet-like phenomena such as spicules (Zhang et al. 2012),
penumbral microjets (Katsukawa et al. 2007), light bridge jets (Asai, Ishii &
Kurokawa 2001; Shimizu et al. 2009), IRIS bombs (Tian et al. 2016).
Chromospheric anemone jets usually occur in solar active regions and the name
anemone is derived from the morphology of coronal X-ray jets reported earlier
(Shibata et al. 1994). The apparent velocity of the chromospheric anemone jets
ranges between 5 km s-1 to 20 km s-1, the lifetime of these jets ranges
between 100 s to 500 s, and the apparent velocity of chromospheric anemone
jets deduced from the SOT/Hinode observations lies somewhat close to the local
Alfvén speed in the lower solar chromosphere (Nishizuka et al. 2011).
Sometimes the jets are found to be invisible in the Calcium-$K_{3}$ (or
Ca-$K_{3}$) component (Morita et al. 2010). Since the Ca-$K_{3}$ component is
emitted from the top of the solar chromosphere, some of the jets that were
reported earlier correspond to the lower chromospheric phenomena (Morita et
al. 2010). However, giant chromospheric anemone jets are also reported in
observations and such jets appear to be associated with the reconnection in
the upper solar chromosphere or solar transition region (Nishizuka et al.
2008).

The chromospheric anemone jets are collimated plasma ejections, and, below the
jet’s spire, there is a presence of a bright loop or a bright knot. In
SOT/Hinode observations, the bright knots are visible in Ca $\scriptsize IIH$
broadband filtergram images. There are blob-like, multiple plasma ejections
that occur in the bright knots and develop gradually into a narrow, collimated
jet (Singh et al. 2012a), usually appearing as a single thread at least within
an arcsec (Nishizuka et al. 2011, Singh et al. 2011, Singh et al. 2012a). It
is believed that the chromospheric anemone jets are formed as a result of the
magnetic reconnection between an ambient magnetic field and an emerging
magnetic bipole (Shibata et al. 2007), and, if the magnetic reconnection
occurs below the middle solar chromosphere or solar photosphere, it indeed
leads to the formation of jets that appear to correspond to chromospheric
anemone jets (Takasao, Isobe & Shibata 2013). Jets on small- scales in the
solar chromosphere are very much intermittent in nature (Singh et al. 2012a,
Tian et al. 2014), however, there is no convincing, self-consistent picture
that establishes a link between the emergence of bright knots, the launching
of small-scale jets from the footpoint region or bright knots during the
evolution of chromospheric anemone jets. The 2D magnetohydrodynamic (MHD)
models alone cannot produce many jets from the entire body of magnetic bipole,
because, in 2D MHD, the jet ejection is localized. In case of 2D MHD, the
magnetic reconnection occurs primarily on one side of the footpoint and does
not erupt during jet activity. Therefore, we need a new 3D model in order to
explain enigmatic features linked with the chromospheric jets.

##  2 Observations

In the present work, we report a high-resolution Calcium $\scriptsize IIH$
filtergram (FG) image of a chromospheric anemone jet that was observed by
SOT/Hinode (Tsuneta et al. 2008, Suematsu et al. 2008, Ichimoto et al. 2008)
aboard Hinode (Kosugi et al. 2007), on 2006 November 20 in NOAA AR10923, in
the northern hemisphere, with a time cadence of about 8 s, and having a
$\lambda$-shaped footpoint from SOT/Hinode.

Figure 1a shows fine-scale, enigmatic structures of one of the chromospheric
anemone jet (Singh et al. 2012b) and the formation mechanism of such features
in the solar chromosphere is still unclear.

In a well-developed Calcium bright knot that is larger than usual, one can
notice fine structures along with a well developed jet spire (main jet).
Figure 1b shows running difference image of the small region indicated by the
box in Figure 1a at 12:40:22 UT -12:40:14 UT to show the small-scale jets
emanating from the bright knots.

Figure 1c shows the snapshot of the small-scale jets and a $\lambda$-shape
region above the bright knots at 12:43:17 UT. Several small-scale jets appear
from the Calcium bright knots. The small-scale jets develops further in to the
solar atmosphere and move systematically from one end of the bright knot.
Later on, these small-scale jets merge into the main jet (see e.g. Figures 1b,
1c, 1d along with online movie1). The main jet is clearly visible in the
running difference images, derived from Calcium filtergram images of
SOT/Hinode (see e.g. Figure 1b, 1d and online movie1). The unsharp masking is
applied to Fig. 1a and 1c. The The faint jets as well as multi-thread
structures are clearly visible in the running difference image (Figure 1d).
The online movie1 provides the evolution of the jets between 12:00 UT to 13:59
UT. The left panel in the online movie1 is obtained by performing the unsharp
masking, the middle panel is obtained by applying the derivative filter, and
the right panel is obtained by applying the running difference method.

The launching positions of the main jet and the small-scale jets are not the
same. The small-scale jets are launched from various locations of the bright
knot that is an extended structure. The jet activity near the bright knots
become pronounced and several small-scale jets are launched from the Calcium
bright knots, after the main jet is developed. The systematic motion of small-
scale jets are seen more clearly when the pieces of bright knots merge into
each other, and form a rather extended structure, at some stage.

The top end of the jet (or the jet spire) in the Calcium image provides the
information of the apparent velocity of the jet. The jet velocity of the main
jet, at a later stage is $\sim$ 20 km/s.

Figure 1e shows the height-time plots of the brightness of small scale jets
along the slit (long rectangle in Figure 1a). In the height-time plots, we are
interested mainly in the small-scale jets located near the footpoint region or
the bright knots. The slit is made wider to accommodate the chosen jets inside
the slit. The analysis of the height-time diagram of the jets near a selected
portion of the Calcium bright knot shows several small-scale jets emanating
from the footpoint, along with a presence of fine-scale features in the
footpoint. The very fast velocities indicated by the dashed white line, in
Figure 1e, is not real but the apparent motions due to the transversal motion
of the jets inside the slit. Nevertheless, the small-scale jets emanating near
the Calcium bright knots have higher velocities ($\sim$ 40 km/s).

##  3 Three Dimensional MHD Simulation of Chromospheric Anemone Jets

To explain enigmatic chromospheric anemone jet shown in Figure 1, Singh et al.
(2012b) proposed the hypothetical model in which many small scale jets are
generated by the successive magnetic reconnection associated with the
emergence of twisted flux tube (as shown in Fig. 4c of Singh et al. 2012b).

Here we try to develop above Singh-model using 3D MHD simulations of
reconnection associated with emergence of twisted flux tube in the solar
active regions.

###  3.1 Details of Numerical experiment

Our model is a direct extension of the 2D MHD model of the coronal jets and
surges developed by Yokoyama & Shibata (1995, 1996). In our self-consistent,
3D MHD simulations, the initial gas layer contains the sub-photosphere
(convection zone or CZ), solar photosphere, chromosphere, transition region,
solar corona. All the gas layers are in the hydrostatic equilibrium under
uniform gravity. The CZ gas layer is convectively unstable. The initial
temperature at the solar photosphere is set to 6000 K, whereas the initial
temperature in the solar corona is set to 6 $\times$ 105 K. The temperature
distribution is akin to the solar atmosphere. Figure 2 shows the initial
distribution of magnetic field and size of the computational box. Note that
$z=0$ is assumed to be the base of the solar photosphere.

The background magnetic field is uniform and parallel in the vertical
direction (Figure 2). A weakly twisted magnetic flux rope is inserted below
the solar photosphere. See Gold & Hoyle (1960) for the initial magnetic flux
rope configuration utilized in our simulations. The axis of the twisted flux
tube is in the x-direction, the positive z-direction is vertically upwards and
the y-direction is perpendicular to x- and z-directions but lies on the xy-
plane.

The initial flux rope configuration is kink unstable. The initial values for
the coronal magnetic field are taken to be 370 G and the magnetic field of the
flux rope at the centre and at $z=-1100$ km is taken to be 5600 G. The twist
in the magnetic flux rope is in the anti-clockwise (right-handed) direction.

A twisted magnetic field configuration has to be utilized for the emerging
magnetic fields in the simulations due to the fact that the S-shaped, sigmoid
structures are often associated with large-scale flux emergence. Solar jets
are often associated with spinning motions (Canfield et al. 1996), so there
always exists a possibility of twisted flux tubes at small scales. Sometimes
the twisted magnetic fields can be formed as a result of local sheared
configuration (e.g. Takasao et al. 2015). The possibility of wave excitation
due to the small-scale, reconnection-generated jet formation in relation to
the chromospheric heating cannot be ignored (Kotani & Shibata 2020).

We have examined the emergence of twisted flux tube in solar atmosphere and
interaction (reconnection) with ambient field, by solving the 3D resistive MHD
equations. We assume adiabatic energy equation and consider resistivity in the
magnetic induction equation. We assumed anomalous resistivity depending upon
the current density for causing fast reconnection (e.g., Yokoyama and Shibata
1996, Nishida et al. 2013). The heating term due to the resistivity is
neglected in the energy equation. The ratio of specific heat is 5/3.

The boundary condition is free at top and side boundaries. The bottom boundary
is assumed to be slipping wall boundary (so that magnetic field and plasma
cannot pass through there but move freely along the boundary).

We did not assume convective flow initially in the simulations. Rather, a
slight upward velocity field is added to the centre of the flux rope. The
initial velocity distribution, between -550 km $<x<$ 550 km, is: $V_{z}$ =
$v_{0}$ cos[($\pi$/10)($x/H$)] , where $v_{0}$ = 7.1 $\times$ 10-3 km s-1, and
$H$=110 km. Here $v_{0}$ can be considered as the minimum enough value to
trigger the emerging flux in the absence of solar convection zone. The
numerical computation is carried out using the CIP-MOCCT scheme (Kudoh,
Matsumoto & Shibata 1999).

###  3.2 Numerical Results

Figure 3a shows the density distribution and provides a simulated image of
jets from our 3D MHD simulations. It is quite evident that the morphology of
the entire jet structure appears differently when the structure is observed
from two different viewing directions (Figure 3a). Figure 3b provides the
temperature distribution of the jet. The 3D computer simulation provides the
clear information of the region present below the jet. An animation (movie2)
of Figures 3a and 3b is available in digital version. The movie2 provides the
jet evolution obtained using density distribution (upper panel) from the
simulation data. The region below the jet shows a clear anemone (or inverse Y)
shape when viewed along the axis of the twisted magnetic flux rope. The online
movie2 also provides the jet evolution obtained using temperature distribution
(lower panel) from the simulation data. The loop or the extended structure
(footpoint) present below the jet has higher density compared to the jet. The
temperature of the footpoint region below the jet has a cooler and a hotter
component, although the majority of the footpoint structure present below the
jet has a cooler component. Such a footpoint region having a higher density
and multi-thermal temperature structure, compared to the overlying jet
structure, could be visualized as a bright loop or bright knot that are
reported quite often in the SOT/Hinode observations.

Figures 3c and 3d provide the absolute velocity of the jets in the yz-plane
for $t=1.94\times 10^{3}$ s. Figure 3c corresponds to $x=0.028$ Mm (i.e. near
the centre of the twisted flux tube) and Figure 3d corresponds to $x=8.8$ Mm
(i.e. away from the centre). The velocity distribution of the jet near the
centre of the flux tube i.e. $x=0.028$ Mm shows that the jet velocity reaches
a value of about 100 km s-1. The velocity distribution of the jet in the outer
region i.e. $x=8.8$ Mm is in the range of 7 km s-1 to 15 km s-1.

Our 3D MHD simulations show that the continuous, upward-rising motion of the
magnetic flux rope creates multiple reconnection sites and results in the jet
formation at various locations and heights in the solar atmosphere. Many
small-scale jets are launched from the body of the flux rope, develop
gradually towards the main jet and merge with the main jet. The evolution of
the main jet and small-scale jets in the computer simulations are very much
similar to the observations of the jets. A weakly twisted magnetic flux rope
provides a guide magnetic field, and the presence of this axial component,
during the upward rising motion of the flux rope, is vital for the launching
of many small-scale jets and their systematic motions from one end of the
bright knot to the other.

The chromospheric anemone jets are embedded in the magnetic field of the
continuously emerging magnetic flux rope. Previous SOT/Hinode observations
show that after the launching of jet activity, the bright knots or anemone-
shaped footpoint start to shrink with time (Singh et al. 2012a), and it is
quite likely that a major portion of the Calcium bright knots, present below
the jets, correspond to emerging flux.

In the present observation, the bright knots do not appear at once. The
different parts of the footpoint (or bright knot) appear in stages, form an
extended structure later on, and eventually gives a $\lambda$ -shape to the
jet. Figure 4a provides the height-time plot of the high-speed, small-scale
jets formed near the bright knots in observations. The height-time plots
corresponding to Figure 4a, obtained from 3D MHD simulation, are provided in
Figure 4b and 4c. The bright parts near the footpoint appear with a speed of
about 6.1 km s-1 or slightly larger. The emerging flux region in the simulated
height-time plot shows motions of about 8.6 km s-1 which is comparable to the
upward speed ($\sim$ 6.1 km s-1) of the rising motion of the bright knot at
the footpoint of the jets. The 3D MHD simulations clearly show vertical jet
velocities as high as 80 km/s.

Fig. 5a shows the magnetic field structure at $t=1.94\times 10^{3}~{}s$. The
magnetic field lines at $t=1.94\times 10^{3}~{}s$ show the reconnected
magnetic field, along with the magnetic fields of the flux rope and plasmoid
formed by magnetic reconnection. It is clear that the reconnected field lines
are twisted. Fig. 5b shows the magnetic field and the velocity field at
$x=0.028$ Mm and $t=1.94\times 10^{3}~{}s$, and during this time, the x-point
is located around $z=1.5$ Mm. The time evolution of the magnetic field and
velocity show the presence of high-speed jets (e.g. around 40 km s-1 to 60 km
s-1) in the computer simulations. It has been found that the height of the
x-point influences the speed of the jet strongly, and the high-speed jets are
formed due to the magnetic reconnection occurring at around $z=1.5$ Mm in the
solar atmosphere. The schematic diagram i.e. Figure 5c to Figure 5j shows the
formation and ejection of jets from the body of flux rope, as viewed from two
different viewing directions.

##  4 Discussion and Conclusion

Understanding the jet velocity, especially the faster velocity components
($\sim$ 40 km s-1 \- 100 km s-1 or more) present in the jets, is fundamental.
The jet velocity ($\sim$10 km s-1) corresponding to the outer region, in our
simulations, is comparable to the velocity of chromospheric anemone jets
(Shibata et al. 2007, Nishizuka et al. 2011), though, in observations,
velocities faster than 40 km s-1 are also reported in some cases (Nishizuka et
al. 2011). It is interesting to note that velocities around 100 km s-1 are
found near the centre of the flux tube in our simulations, much faster than
those of typical chromospheric anemone jets. On the other hand, in one
observation (Nishizuka et al. 2008) where the size of the emerging bipole was
relatively large, the jet was accelerated to around 100 km s-1. Overall, a
careful estimate of the small-scale jet velocities in our particular
observations, near the Calcium bright knots, gives a higher value compared to
the usual values reported earlier during the statistical studies of
chromospheric anemone jets (e.g. Nishizuka et al 2011). The apparent velocity
of small-scale jets through eye-detection is $\sim$ 40 km s-1 \- 100 km s-1,
but higher velocity components are also possible.

The Calcium bright knots before the appearance of jetting are studied in
observations. The high-resolution Calcium filtergram observations from
SOT/Hinode have revealed that there exists a significant time lag between the
bright point or Calcium loop/knot and the appearance of the main jet, see e.g.
Table 1 of Singh et al. (2012a), where SOT/Hinode observations show that the
time- lag between the appearance of the bright loop and appearance of the main
jet ranges between 1 min- 10 min, in fact, the time-lag between 1.5 min - 3
min is more common. In our present observations where the Calcium bright knot
appears in parts that merge gradually, the time-lag is not fixed with typical
time-lag between 3 min - 4 min. It is interesting to see that in our 3D MHD
simulations, the typical time- lag between the appearance of the loop and the
main jet is $\sim 3.2$ min. The time-lag noticed in our computer modelling
should also be applicable to emerging flux and jets.

Multiple jets could be launched from a location if there is a recurrent
emergence of bipoles in a small region of the solar chromosphere. In the
present work, recurrent jets are produced in our 3D MHD simulations, and we
provide a different scenario where there would be a presence of a single,
dynamically evolving bright knots (or bipole), and many small-scale jets are
launched as a result of multiple magnetic reconnection occurring in the entire
body of the bright knots, and such a scenario fits very well with the
observations of the chromospheric anemone jets. There are jets associated with
light bridges that are observed in the umbral region of the sunspots (Shimizu
et al. 2009, Shimizu 2011, Asai, Ishii & Kurokawa 2001, Schmieder 2022). The
morphology of such jets is similar to the successive ejections or small scale
jets associated with chromospheric anemone jets, but the latter occur much
more frequently than the light bridge jets.

Interestingly, in our self-consistent simulations, our new model can describe
jets not only occur near the mixed polarity region (at the level of solar
photosphere) but elsewhere in the body of the flux rope, exhibiting a thread-
like pattern, and the footpoint continues to evolve in the upward direction
and horizontal expansion takes place during the launching of the small-scale
jets. The velocity distribution of the jet on the yz-plane at different
x-positions in the simulation domain provides different values of the jet
velocity. While developing interpretations of solar jets that are based only
on the examination of the observed morphological features, one should be
cautious because the jet formation requires a proper understanding of the
magnetic reconnection (Isobe et al. 2005; Isobe, Proctor & Weiss 2008) along
with the driving mechanisms that can prevail at various heights in the solar
atmosphere (Isobe et al. 2005; Shibata et al. 2007; Tian et al. 2014).
Observations show that the nature of magnetic reconnection is quite bursty and
time-dependent in solar chromosphere (Singh et al. 2012a). The height of the
magnetic reconnection and interaction of slow mode shocks with the solar
transition region is important (Takasao, Isobe & Shibata 2013) for
chromospheric jets and such effects should be studied in 3D MHD simulations.

The initial velocity trigger (by a small amount) is near the centre of the
flux rope, so the middle part of flux rope will start rising and appear
earlier in the higher atmosphere. The magnetic reconnection occurs also in the
solar chromosphere in our 3D MHD model. As time goes on, the magnetic
reconnection near the centre of the flux rope would take place at higher
heights in the solar atmosphere, compared to regions of the flux rope that are
located near the outer regions of the flux rope, and the main jet appears. In
this situation, the small-scale jets would move towards the main jet. If
location of the initial velocity perturbation is near the outer region of the
flux rope (initial, twisted magnetic field), the outer part will emerge first
at the higher heights, and one would then expect higher velocities of small-
scale jets in the outer region, because the magnetic reconnection will occur
at higher heights in the solar atmosphere. Even though the location of the
initial velocity perturbation is important, the underlying physics would not
be affected much.

The faster velocity jets ($\sim 100$ km s-1) could be due to the fact that the
magnetic reconnection takes place in higher layers where the local Alfvén
speed is high. The magnetic field and density both decrease with height. The
ratio of magnetic fields in the upper and middle chromosphere is $\frac{B_{\rm
upper}}{B_{\rm middle}}\sim\frac{1}{5}$ but the ratio of the number densities
would be $\frac{n_{\rm upper}}{n_{\rm middle}}\sim 10^{-3}$. So, if we take
the ratio of the local Alfvén speeds at two different heights i.e. upper
chromosphere and middle chromosphere, the local Alfvén speed in the upper
chromosphere should be 10 times higher than the local Alfvén speed in the
middle chromosphere. Therefore, one would expect jets formed in the upper
chromosphere to have higher velocities than jets formed below the middle
chromosphere.

We now discuss the effect of our assumption, especially adiabatic energy
equations in 3D MHD simulations. The main reason of this simplified assumption
is the difficulty of assuming realistic energy equations for chromospheric
plasma, since radiation cooling is controlled by non-LTE radiative transfer.
On the other hand, the gravitational stratification plays an important role in
chromospheric jet formation and subsequent acceleration (e.g. Takasao, Isobe &
Shibata 2013). Even though we do not include radiation, we consider that the
gravitational stratification provides a much stronger effect on the jet
dynamics in the solar chromosphere. This is because the characteristic
velocity of the jet is determined by the local Alfvén speed and the density
determines the local Alfvén speed. The temperature can indirectly affect
density distribution, e.g., if pressure balance is satisfied, and if the
temperature increases twice, then density decreases to half, leading to an
increase in Alfvén speed by $\sqrt{2}$. However, even if temperature changes,
if density does not change, the Alfvén speed remains the same, so there is no
influence on jet dynamics. The resolution ($\sim$ 55 km) is sufficient enough
to resolve the small-scale jets. Our model is the first step of 3D MHD
modeling of enigmatic chromospheric jets, and we believe, we can capture
essential physics of dynamics associated with magnetic reconnection and
emergence of twisted flux tube even under the assumption of adiabatic energy
equation. The inclusion of realistic energy equation will be important future
task.

The coronal physical quantities such as temperature, density, and gas pressure
are not realistic values in our initial conditions, as in the 2D model of
Yokoyama & Shibata (1995, 1996). The temperature of solar corona is lower than
typical coronal temperature in active region (i.e. 1 $\times$ 106 K to 2
$\times$ 106 K ). The height of the solar transition region is much less than
the realistic value. The density of coronal plasma is much greater than the
realistic value (see e.g. Yokoyama & Shibata 1996). Also see two-dimensional
MHD simulations of Nishizuka et al. (2008) where a realistic coronal density
and temperature is adopted compared to Yokoyama & Shibata (1996), but it bring
out similar physics. These demerits are not essential for the modeling of 3D
jets. The presence of radiation in the energy equation would lead to a more
realistic temperature distribution, though dynamics will remain unaffected,
driven as it is by the Lorentz force. In the future, more realistic 3D
modeling of jets will be made using large computational resources.

In conclusion, the present work provides insight through our new model of
chromospheric anemone jets. For the first time, it has become very much
evident, that the continuous, upward rising of the twisted magnetic flux rope
in the magnetized solar atmosphere, is the key behind the formation of the
chromospheric anemone jets, and the new MHD model has the capacity to resolve
diverse, intricate issues related with the chromospheric anemone jets and its
footpoint evolution. Under such circumstances, the jet activity in the solar
chromosphere has to be understood from the viewpoint of the driving mechanisms
and the underlying jet dynamics common to large-scale flares and mass
ejections.

## Acknowledgment

Hinode is a Japanese mission developed and launched by ISAS/JAXA, with NAOJ as
a domestic partner and NASA and STFC (UK) as international partners. It is
operated by these agencies in cooperation with ESA and NSC (Norway). This
research is supported by JSPS KAKENHI grant number 21H01131 (K. Shibata). KAPS
gratefully acknowledges the UGC Faculty Recharge Program of the Ministry of
Human Resource Development (MHRD), Govt. of India and University Grants
Commission (UGC), New Delhi, the incentive grant of the Institute of Eminence
(IoE) Program, BHU, and the Visiting Associateship Program of IUCAA, Pune. The
numerical computations were carried out on SX9 at the Center for Computational
Astrophysics (CfCA) of the National Astronomical Observatory of Japan, and the
KDK system at the Research Institute for Sustainable Humanosphere (RISH) at
Kyoto University. The authors are grateful to M. Shimizu, K. Uehara, and all
the staff, members of the Astronomical Observatory of the Graduate School of
Science, Kyoto University, Japan. We thank Prof. Marco Velli, Department of
Earth, Planetary and Space Sciences, University of California, Los Angeles
(UCLA) for constructive comments.

## References

  * Asai et al. (2001) Asai, A., Ishii, T. & Kurokawa, H. 2001, ApJ, 555, L65-L68
  * Canfield et al. (1994) Canfield, R.C. et al. 1996, ApJ, 464, 1016
  * Gold et a. (1960) Gold, T. & Hoyle, F. 1960, MNRAS, 120, 89
  * Ichimoto et al. (2008) Ichimoto, K. et al. 2008, Sol. Phys. 249, 233
  * Isobe et al. (2005) Isobe, H., Miyagoshi, T., Shibata, K. & Yokoyama, T. 2005, Nature, 434, 478
  * Isobe et al. (2008) Isobe, H., Proctor, M.R.E. & Weiss, N.O. 2008, ApJ, 679, L57
  * Katsukawa et al. (2007) Katsukawa, Y et al. 2007, Science, 318, 1594
  * Kosugi et al. (2007) Kosugi, T. et al. 2007, Sol. Phys. 243, 3
  * Kotani et al. (2020) Kotani, Y., & Shibata, K. 2020, PASJ, 72, 75 (22pp.)
  * Kudoh et al. (1999) Kudoh, T., Matsumoto, R., & Shibata, K. 1999, Computational Fluid Dynamics Journal, 8, 56
  * Morita et al. (2010) Morita, S. et al. 2010, PASJ, 62, 901
  * Nishida et al. (2013) Nishida, K., Nishizuka, N. & Shibata, K. 2013, ApJ, 775, L39 (6pp.)
  * Nishizuka et al. (2008) Nishizuka, N. et al. 2008, ApJ, 683, L83
  * Nishizuka et al. (2011) Nishizuka, N. et al. 2011, ApJ, 731, 43(11pp.)
  * Raouafi et al. (2016) Raouafi, N. E. 2016, Space Sci. Rev., 201, 1-53 (2016)
  * Schmieder et al. (2022) Schmieder, B. 2022, Frontier in Astron. and Space Sci., 9, Article 820183 (15pp.)
  * Shibata et al. (1994) Shibata, K. et al. 1994, ApJ, 431, L51
  * Shibata et al. (2007) Shibata, K. et al. 2007, Science, 318, 1591
  * Shibata et al. (2011) Shibata, K. & Magara, T. 2011, Living Rev. in Sol. Phys., 8, 6 (99pp.)
  * Shimizu et al. (2009) Shimizu, T. et al. 2009, ApJ, 696, L66-L69
  * Shimizu (2011) Shimizu, T. 2011 ApJ, 738, 83
  * Singh et al. (2011) Singh, K. A. P., Shibata, K., Nishizuka, N. & Isobe, H. 2011, Phys. Plasmas, 18, 111210 (8pp.)
  * Singh et al. (2012a) Singh, K.A.P. et al. 2012a, ApJ, 759, 33 (14pp.)
  * Singh et al. (2012) Singh, K.A.P., Isobe, H., Nishida, K. & Shibata, K. 2012b, ApJ, 760, 28 (5pp.)
  * Sterling et al. (2015) Sterling, A.C. et al. 2015, Nature, 523, 437
  * Sterling et al. (2020) Sterling, A.C. and Moore, R.L. 2020, ApJ, 896, L18 (9pp.)
  * Suematsu et al. (2008) Suematsu, Y. et al. 2008, Sol. Phys. 249, 197
  * Takasao et al. (2013) Takasao, S., Isobe, H., Shibata, K. 2013, PASJ, 62 (22pp.)
  * Takasao et al. (2015) Takasao, S., Fan, Y., Cheung, Mark C. M. & Shibata, K. 2015, ApJ, 813, 112 (14pp.)
  * Tian et al. (2014) Tian, H. et al. 2014, Science 346, 1255711 (4pp.)
  * Tian et al. (2016) Tian, H., Xu, Z., He, J. & Madsen, C. 2016, ApJ, 824, 96 (14pp.)
  * Tsuneta et al. (2008) Tsuneta, S. et al. 2008, Sol. Phys. 249, 167
  * Yokoyama et al. (21995) Yokoyama, Y., & Shibata, K. 1995, Nature, 375, 42
  * Yokoyama et al. (1996) Yokoyama, Y., & Shibata, K. 1996, PASJ, 48, 353
  * Zhang et al. (2012) Zhang, Y.Z. et al. 2012, ApJ, 750, 16 (9pp.)

Figure 1:  The small-scale jets from the bright knot in Ca $\scriptsize IIH$
broadband filtergram data obtained from high-resolution data from SOT/Hinode.
The rectangular box with a solid line in the image shows the slit position.
The rectangular slit is 10 pixels wide and it is utilized in the making of a
height-time plot later on. (a) The main jet along with a bright knot in
enigmatic solar chromosphere during high-resolution observations, as seen from
SOT aboard Hinode satellite. (b) Running difference image to show the small-
scale jets emanating from the bright knots. The running difference images are
derived from the observations of Calcium filtergram images taken from
SOT/Hinode to show the small-scale jets jutting from Calcium bright knot, also
see Singh et al. (2012b) (c) Ca $\scriptsize IIH$ image to show the small-
scale jets and $\lambda$-shape at 12:43:17 UT, (d) Running difference image to
show that, along with a main jet, there exist small-scale jets along with
multi-thread structures. The multi-thread structures at a later stage are more
clear. The field of view of Figs. 1b to 1d are shown in Fig. 1a with a dotted
rectangular box. Scales are provided for the sake of comparisons, (e) The
height-time plot, derived from Calcium filtergram images (in observations), of
a region below the observed jet that is present in the calcium bright knot (or
footpoint). The rectangular slit position is shown in Fig. 1a. Typical
velocities are shown with white ticks in the height-time plot. The very fast
velocities indicated by dashed white lines, on Fig. 1e, are not real but
apparent motion due to transversal motion of jets on the slit. Also, see
online movie1 where the Field-of-view (FOV) of the jets, obtained from the Ca
$\scriptsize IIH$ filtergram data, is shown in Fig.1a with a dotted
rectangular box, and jets arise between 12:00 UT to 13:59 UT. Figure 2: The
initial magnetic flux rope configuration along with the uniform, ambient
magnetic field utilized in 3D MHD simulation. As time goes on, the initial
magnetic flux rope configuration evolves dynamically during the upward rising
and at some point in time, the jets are formed. Figure 3: The small jets along
with the first jet (or main jet) are also noticed in 3D MHD simulations. (a)
The density distribution from 3D MHD simulation shows several such jets
launched from the footpoint. The density distribution at $t=1.94\times 10^{3}$
s is shown for two different viewing directions, (left: along the axis of the
twisted magnetic flux tube, yz-plane, x = 0.028 Mm,(right: perpendicular to
the axis of the twisted magnetic flux tube, xz-plane, y = 0.028 Mm (b) The
temperature distribution at $t=1.94\times 10^{3}$ s, $y=0.028$ Mm for viewing
direction perpendicular to the jet axis is shown. (c and d) The absolute
velocity distribution of the jets on yz-plane. The figures show two cases for
$t=1.94\times 10^{3}$ s, panel (c): $x=0.028$ Mm, panel (d): $x=8.8$ Mm. The
$x=0.028$ Mm inside the flux rope is the position near the axis of the flux
rope and the jets have higher velocity components. In the outer region i.e.
$x=8.8$ Mm the jets have relatively smaller velocity components. The white box
in Fig. 3(c) corresponds to the region shown later on in Fig. 5(b). An
animation (movie2) of Figs. 3a and 3b is available in the digital version. The
upper panel in movie2 provides the time evolution of jets using density
distribution obtained from 3D MHD simulation data. The lower panel in movie2
provides the time evolution of jets using temperature distribution obtained
from 3D MHD simulation data. Figure 4:  (a) The height-time plot, derived from
Calcium $\scriptsize IIH$ filtergram images a selected region near the Calcium
bright knot. The rectangular box is the same as the one used for Fig. 1e.
Velocities are shown with white ticks in the height-time plot. The observed
height-time plot can be utilized for comparisons with MHD simulations. The
very fast velocities indicated by dashed white lines are not real but apparent
motion due to transversal motion of jets on the slit, (b) a simulated height-
time plot, derived from 3D MHD modeling, of the jet using the density
distribution, (c) a simulated height-time plot of the jet using the velocity
distribution. The contours in the velocity distribution shows the velocities
in units of km s-1. Both (b) and (c) correspond to $x=0.58$ Mm, $y=-0.47$ Mm.
Figure 5: (a) The magnetic field structure at $t=1.94\times 10^{3}$ s obtained
from the simulation data. The iso-surface of density is shown by grey
colour.(b) The magnetic field and velocity field are shown at $x=0.028$ Mm and
$t=1.94\times 10^{3}$ s. The magnetic field and velocity field, obtained from
simulation data, are shown in the panel (b) here (i.e. Magnetic field and
velocity field at $t=1.94\times 10^{3}$ s and $x=0.028$ Mm) correspond to the
box shown in Fig. 3c. The schematic diagrams (c to j) show the origin and
evolution of the jets, the oblique view (c to f) and axial view (g to j) are
provided in the diagrams. The magnetic flux rope is shown by grey colour. The
sign ’X’ in the schematic diagram shows the reconnection point, the flux rope
magnetic field is shown in orange color and the ambient magnetic field is
shown in green color. At first, Jet1 is formed and magnetic reconnection
proceeds further leading to more jets such as Jet2 and so on. In this manner,
many small-scale jets would be launched from the entire body of the magnetic
flux rope.


††thanks: e-mail<EMAIL_ADDRESS>
# Simulating COVID19 Transmission From Observed Movement: An Agent-Based Model
of Classroom Dispersion

Yi Zhang  University of Miami, Department of Physics, Coral Gables, FL  Yudong
Tao  University of Miami, Department of Electrical and Computer Engineering,
Coral Gables, FL  Mei-Ling Shyu  University of Miami, Department of Electrical
and Computer Engineering, Coral Gables, FL  Lynn K. Perry  University of
Miami, Department of Psychology, Coral Gables, FL  Prem R. Warde  Care
Transformation, University of Miami Hospitals and Clinics, Miami, FL  Data
Analytics Research Team (DART) Research Group, University of Miami Hospitals
and Clinics, Miami, FL  Daniel S. Messinger  University of Miami, Department
of Psychology, Coral Gables, FL  Chaoming Song  University of Miami,
Department of Physics, Coral Gables, FL

###### Abstract

Current models of COVID-19 transmission predict infection from reported or
assumed interactions. Here we leverage high-resolution observations of
interaction to simulate infectious processes. Ultra-Wide Radio Frequency
Identification (RFID) systems were employed to track the real-time physical
movements and directional orientation of children and their teachers in 4
preschool classes over a total of 34 observations. An agent-based transmission
model combined observed interaction patterns (individual distance and
orientation) with CDC-published risk guidelines to estimate the transmission
impact of an infected patient zero attending class on the proportion of
overall infections, the average transmission rate, and the time lag to the
appearance of symptomatic individuals. These metrics highlighted the
prophylactic role of decreased classroom density and teacher vaccinations.
Reduction of classroom density to half capacity was associated with an 18.2%
drop in overall infection proportion while teacher vaccination receipt was
associated with a 25.3% drop. Simulation results of classroom transmission
dynamics may inform public policy in the face of COVID-19 and similar
infectious threats.

## Introduction

On March 11, 2020, the World Health Organization declared COVID-19 a pandemic
and called for coordinated mechanisms to support preparedness and response
efforts across health sectors[1, 2, 3, 4]. Predictive models can effectively
inform policy to coordinate social policy responses to infectious
pandemics[5]. The CDC has embraced the use of mitigation strategies in schools
to allow communities to keep preschools and K-12 schools open[6, 7]. Here we
focus on simulating infection responses to the novel coronavirus (SARS-CoV-2,
which is responsible for COVID-19) based on observed interactions in a
classroom setting. To our knowledge, this is the first study to apply
Susceptible, Exposure, Infected, Recovered (SEIR) to physically-defined
interactions in these contexts. Results may inform policy decisions related to
vaccination, school attendance, and school closure and reopening.

SEIR models are used to quantify the spread of epidemics at a societal
level[8, 9]. SEIR models of COVID-19 quantify growth rate with respect to
sample and population characteristics that may change with time[10]. Given
their population focus, these models do not characterize individuals as
transmission agents, and do not model their activity over time in a physical
space inhabited by other agents. As SEIR models do not model granular
transmission processes, they may be limited in the degree to which they speak
to policy initiatives involving the repeated activity of individuals in
physical space, such as schools.

While movement data collected by GPS [11, 12], mobile phone records [13, 14,
15, 16], and public transportation [17] has been used to investigate the
transmission of SARS-CoV-2 and other viruses, the spatial resolution of those
sensors is largely limited to ranges of meters to kilometers, which are
unsuitable to the study of classroom and other enclosed space transmission.
Two recent reports demonstrate the potential of more granular analysis. One
analysis modeled reported infection outcomes based on estimates of exposure
duration and the physical density of individuals derived from published
outbreak data[18]. This report contributed a model of infection saturation in
small, relatively defined samples, but did not quantify transmission dynamics
between individuals in physical space. Likewise, a recent model of
transmission at a conference event quantified the importance of contact
duration during which a specific pair of individuals were within an a priori
(badge-based) physical range[19]. Similar techniques have been used to study
classroom transmissions [20]. In contrast to the current effort, physical
distance was not measured continuously nor was the relative physical
orientation of individuals modeled.

Data from 2020 indicate that young children may acquire COVID-19 (including
confirmed, asymptomatic cases) in childcare settings and subsequently infect
their household members[21]. Moreover, contact tracing revealed high levels of
transmission between individuals of similar ages, including among children
younger than four years of age[22]. Thus there is evidence that young children
may serve as vectors for COVID-19, Badge-based technology has been used to
explore the temporal dynamics of transmission among high-school students[20],
but it is not clear how physical motion facilitates infection dynamics in
classrooms. The current study capitalizes on automated observations of pre-
COVID-19 classroom physical interaction to model transmission between children
and their teachers. It quantifies the dynamics of potential classroom
outbreaks to help support efforts to mitigate outbreaks. To do so, we simulate
the consequences of both half-capacity classrooms and teacher vaccinations on
infection outcomes. These outcomes include infection saturation (percent and
number of infected individuals) over time, time to the emergence of the first,
second, and third symptomatic individual, and infection saturation levels when
these symptomatic individuals emerge (a potential trigger of school closure).

## Methods

### Data Sources and Calibration

Figure 1: (a) Ubisense tag worn by the teachers and children in the
classrooms; Layout of Classrooms for (b) Class A; (c) Class B and (d) Classes
C & D. (Classes C & D occupied the same classroom in successive years).

Simulations were based on continuous real-time location tracking of
individuals (students and teachers). A total of 34 observation periods
occurred for four preschool classes housed in three classrooms in a large
urban center in the United States. Observations were conducted in inclusion
classrooms for typically developing children and children with developmental
disabilities contained in both a standard public county-funded school (Classes
C and D) and a county-funded university preschool (Classes A and B). The range
of represented class sizes (10-17) and teacher-child ratio are typical of
preschool classrooms, as state regulations require preschool classrooms to
have fewer than 20 children and teacher:child ratios of at least 10:1 [23].
Ultra-Wide Radio Frequency Identification (RFID) technology such as that
embedded in child-worn Ubisense tags Fig. 1(a), allowed for efficient capture
of individuals’ location and movement in the classroom indicating when
children and teachers are in proximity. The Ubisense Dimension system was used
to track individuals’ location at $2-4$ Hz to an accuracy of $21$ cm in the
classrooms[24]. The system consists of one sensor in each corner of the
classroom, a dedicated server, and active tags worn by children and teachers.
All children and teachers in the classrooms were tracked with two tags worn on
their left and right hips (children wore a specially constructed vest and
teachers wore fanny packs). The tags’ ultra-wide-band RFID signals were used
to calculate the locations and orientations of both children and teachers in a
three-dimensional space by means of triangulation and time differences in
arrival. Specifically, each individuals’ locations were estimated as the
centroids of their two tags while their orientations were calculated from the
relative positions of the left and right tags.

Fifty 2-5-year-old children enrolled in four preschool classes and their 11
teachers participated in a study in which continuous RFID measurements of each
individual’s location were collected. The classes were observed in academic
years 2018-2019 and 2019-2020. The last observation took place on 2/28/2020,
prior to the impact of COVID-19. Recordings were collected during observations
of the entire school-day ($2.5-3.5$ hours), spaced 2-4 weeks apart in each
class. Figure 1(b)–(d) plots the layout of the classrooms. Class
characteristics are reported in Table 1.

Class |  Class A |  Class B |  Class C1 |  Class D
---|---|---|---|---
Number of observations |  13 |  12 |  5 |  4
Average data collection length |  3.1 hours |  3.2 hours |  2.2 hours |  2.3 hours
Total number of children (girls) |  17 (8) |  10 (4) |  12 (5) |  11 (7)
Total number of teachers |  3 |  2 |  3 |  3
Table 1: Sample characteristics by class.

### Ethics Declarations

The study was conducted in accordance with APA ethical standards and was
approved by the University of Miami Institutional Review Board (#20160509).
All teachers gave their informed consent and received $100 for their
participation. Parents gave their informed consent on behalf of their children
and received $75 for their participation. Teachers also received a classroom
gift of their choice (e.g., a new book shelf). All research was performed in
accordance with relevant guidelines and regulations.

### Transmission Model

We modelled the spread of SARS-CoV-2 in the classrooms using a novel
modification of the SEIR model. The typical SEIR model consists of
susceptible, exposed, infectious, and recovered individuals where the
susceptible individual becomes infectious after close contact with an
infectious individual, with infection rate $\beta$, while infectious
individuals recover at rate $\gamma$. For each individual i, the real-time
classroom movement data records the position $(x_{i},y_{i})$ and the
orientation $\theta_{i}$ for every second. To integrate the simplest SEIR
model with the movement data, we develop a model to account for the dependence
of the infection rate $\beta$ on distance
$r_{ij}=\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}}$ and the time difference
$t=t_{1}-t_{2}$ between pairs of individuals ($t_{1}>t_{2}$) with
interpersonal orientations $(\theta_{i}\&\theta_{j})$, satisfying

$\beta(r_{ij},\theta_{i},\theta_{j},t)=\beta_{max}\exp\left(\frac{-r_{ij}^{2}}{2\sigma_{r}^{2}}-\frac{\theta_{i}^{2}+\theta_{j}^{2}}{2\sigma_{\theta}^{2}}\right)e^{-\lambda
t},$ (1)

with $\sigma_{r}\approx 2m$, and $\sigma_{\theta}\approx 45^{0}$ where
$\beta_{max}$ is the maximum pairwise infection rate. That is, the rate of
infection falls directly as an exponential function of growth in the square of
the radius between the two individuals and their angular distance from face-
to-face contact[25]. A visual illustration of the distance $r_{ij}$ and
orientations $\theta_{i}$ and $\theta_{j}$ is shown in Fig. 2(a), and the
heatmap and the 3D plot of $(r_{ij},\theta_{i},\theta_{j})$ with
$x=r_{ij}cos\theta_{i}$, $y=r_{ij}sin\theta_{j}$ and $\theta_{i}=0$ are
presented in Fig. 2(b) and (c), respectively. A detailed discussion of the
infection function is found in SM Section 4.C and 4.D. The $e^{-\lambda t}$
term accounts for temporal decay, modeling declining transmission over time
$t$ which characterizes the decay rate. For airborne transmission,
$\lambda=0.34/\mathrm{hours}$ [26], whereas droplet transmission decay is
functionally instantaneous over a timescale of seconds (see SM Section 4.A for
more details).

To calibrate $\beta_{max}$ of the agent-based SEIR model, we calculate the
average infection rate

$\bar{\beta}=\frac{N\iiiint\
\beta_{max}\exp\left(-\frac{x^{2}+y^{2}}{2\sigma_{r}^{2}}-\frac{\theta_{i}^{2}+\theta_{j}^{2}}{2\sigma_{\theta}^{2}}\right)\,dx\,dy\,d\theta_{i}\,d\theta_{j}}{A(2\pi)^{2}}=\beta_{max}\sigma_{r}^{2}\sigma_{\theta}^{2}\rho.$
(2)

Note that the population infection rate $\rho_{0}$ is proportional to the
density $\rho\equiv\frac{N}{A}$ whereas $\beta_{max}$ is an intrinsic
parameter, independent of the social environment. Researchers have estimated
the average daily infection rate at the meta-population level daily
$\bar{\beta}_{daily}\approx\beta_{0}$, where $\beta_{0}=R_{0}\gamma$ with
$R_{0}$ and $\gamma$ being the reproduction number and the recovery rate ,
respectively. Combining Eq. (2), we calibrate our modeling parameter

$\beta_{max}=\frac{\bar{\beta}_{daily}}{\sigma_{r}^{2}\sigma_{\theta}^{2}\rho_{daily}}.$
(3)

SARS-CoV-2 transmission has been modelled and analyzed in many literatures
while we posit that the reported $R_{0}^{\prime}s$ variations in these works
are largely due to virus evolution/mutation, as well as social, political, and
environmental between each cohort[27, 28, 29, 30, 31, 32]. We chose a
conservative lower bound of $R_{0}=2.0$ for our numeric simulations based on
reported $R_{0}$ values, which typically fall within the $2.0-3.0$ range[28,
33, 29, 30, 31, 32]. CDC guidelines have suggested defining close contact as
being within $r=6$ feet of an infected person for a cumulative total of $T=15$
minutes or more over a 24-hour period[34]. We thus estimate the daily
population density $\rho_{daily}=\frac{N_{c}}{\pi
r^{2}}\times(T/24\mathrm{hours})$ in Eq. (3) with average number of contacts
$N_{c}\approx 10$[35].

Figure 2: (a) The definition of orientation for each pair of children,
$r_{ij}$ is the distance between two children, and the orientation for a child
$i$ is defined by angle $\theta_{i}$, which is the angle between the direction
being faced (red arrow) and the direct line between persons (blue line). A
similar concept is applied to child $j$ to obtain $\theta_{j}$. (b) The
heatmap and (c) the 3D plot of $\beta(r_{ij},\theta_{i},\theta_{j})$ with
$x=r_{ij}\cos\theta_{j}$, $y=r_{ij}\sin\theta_{j}$ and $\theta_{i}=0$ where
child $i$ is positioned at the origin and facing graph right. (d) The
transition from infection to becoming infectious was $1$ day; the incubation
period (from infection to symptomaticity) had a mean of $4$ days; the mean
recovery time (from infection to infectious to recovered, i.e., non-
infectious) had a mean of $10$ days.

### Numeric Simulations

Based on the agent-based SEIR model, the spread of SARS-CoV-2 in classrooms
can be simulated on a continuous-time basis. Here, infection only occurs
during real time classroom observations that are separated by $24$ hours for
each weekday (Monday through Friday) and $72$ hours between Friday and Monday.
For each simulation, we start with a randomly chosen patient zero who is
infectious at Day 0. For each time step $\Delta t=1$ second, we use the
empirically measured location and orientation of each individual during the
class time to compute the relative distance $r$ and orientations
$\theta_{i}\&\theta_{j}$ for every pair of individuals $i$ and $j$ (including
both teachers and children). Each susceptible individual is infected by
infectious peers with a probability $\beta\Delta t$ calculated from all
infected neighbours, where $\beta(r_{ij},\theta_{ij})$ is calculated based on
Eq. (1). The status of all individuals is updated every second based on the
real-time location and orientation data. We illustrate the time dependent
parameters related to infectious processes in Fig. 2 (d). Individuals become
infectious $24$ hours after being infected[36]. Once infected, individuals
have a 75% probability of becoming symptomatic[37]; the transition to
symptomaticity followed a poisson process with a mean of four days[38, 39,
28]. Infected individuals become non-infectious with a probability
$\gamma\Delta t$ where $\gamma=1/10$ ${days}^{-1}$ at each time step,
reflecting a mean $10$-day duration[40]. We consider non-infectious
individuals to be recovered, and assume no repeated infections. Each
simulation continued until all infected individuals are recovered. All time
periods include time in and out of the classroom.

We applied the agent-based SEIR model to real-time location data from four
preschool classes A-D. For each observation, we began with a single patient
zero who was infected at $t=0$ and simulated the observed classroom
interaction over $60$ iterations. This process was repeated for each
individual in the classroom (each was used as patient zero). For each
simulation we output the (S, E, I, R) status of each individual and a mean
individual hourly status is calculated over all simulations.

### Public Health Scenarios

We explored four public health scenarios pertaining to classroom density and
teacher vaccination. To explore the role of classroom density, SARS-CoV-2
infection was simulated for both full and half sized classroom scenarios. In
the full class scenario, every student and teacher was simulated. In the half-
class scenario, half the students in the class (rounding up for odd numbers of
students) and one teacher were simulated. For example, if an empirical
observation of Class A contained 18 children and 3 teachers, the half sized
class simulations would involve 9 children and 1 teacher. The half-class
simulations simultaneously model in which class sizes are administratively
reduced and scenarios in which fewer students opt to attend class.

At the time of writing, vaccination against COVID-19 was available for
individuals above the age of 12, but not for younger individuals[41]. Thus we
conducted a parallel set of simulations assuming that classroom teachers were
vaccinated. When vaccinated, we assumed a $85.8\%$ probability of vaccine
effectiveness ($85.8\%$ probability of no infectivity) for each simulation
run. The probability of effectiveness was based on the average efficacy of
FDA-approved vaccines at the time of writing (Pfizer $97\%,$Moderna $94.1\%$,
Johnson & Johnson $66.3\%$)[42, 43, 44]. For each of the four scenarios formed
by the crosstabulation of classroom density and teacher vaccination levels,
numeric simulations for each patient zero were repeated $60$ times. There were
no infections in simulations involved a vaccinated teacher (whose vaccination
was effective) as the patient zero. For both full and half size class
scenarios, the simulations are conditional on one infected person attending
class. However, if we consider a school with all half-sized classes and double
the number of classes, the chance of each class having a patient zero is
reduced by a factor of two, which reduces the infections further.

### Policy-Focused Outcomes

Numerical simulations were intended to inform decisions and policy designed to
curb the impact of COVID-19 and other infectious agents in classrooms. We
first calculated the time course of infection saturation in full-class and
half-class scenarios with and without teacher vaccination over 28 days. Next,
we estimated the proportion of simulations yielding a first, second, and third
symptomatic individual, and the number of days required for their emergence.
We reasoned that the emergence of a symptomatic child or teacher would be the
first overt sign of COVID-19 infection in a classroom. Thus the number of days
until a first symptomatic case emerged would index ‘damage done’ from a policy
perspective. Likewise, the time to subsequent symptomatic cases would indicate
the cost of ignoring a first infection.

Figure 3: SARS-CoV-2 Infection over Time. The proportion (right) and total
number (left) of infections over time in days . Trajectories are simulated for
four classes for the full versus half-sized scenarios and not vaccinated
versus teacher-vaccinated scenarios. The red curves and grey areas represent
the mean and standard deviation of the proportion/total number of infections
for all observations, respectively. $\rho$’s in each legend represent the
classroom densities in the unit of people per $m^{2}$.

## Results

### Infection Over Time

To investigate the spread of SARS-CoV-2, we plot the proportion (and number)
of infected individuals over one month in Fig. 3 for both full/half classes
and non-vaccinated and teacher-vaccinated scenarios. The mean level (red line)
and standard deviation of infection (grey area) over simulation runs is
presented. The classroom population density, $\rho$, the number of individuals
per square meter of classroom space is calculated for each observation (see SM
Section 2 for more details). Figure 3 shows lower infection levels over time
in half-sized than full size simulations. Likewise, the teacher vaccinated
scenarios yield lower infection levels over time than the not vaccinated
scenarios. In addition, Figure 3 suggests an association between the classroom
density $\rho$ and infection levels over time. Specifically, scenarios with
higher densities produced higher proportions of infected individuals. These
findings suggest that classroom density plays an important role in controlling
SARS-CoV-2 spread.

### Transmission Likelihood and Saturation

To quantitatively investigate the impact of classroom density, we characterize
the infection patterns of the numerical simulations using the methodology of
Tupper et al[18]. In the current model, the time-dependent infection number is
determined directly by the agent-based simulation (see Fig. 3), providing
exposure duration and similar transmission measures[18]. Here we focus on the
association between two transmission measures, saturation and average
transmission rate. Saturation is the cumulative number of infected individuals
(normalized by the total population) at the end of the simulation (final
proportion infected). The transmission rate $\hat{\beta}$ is computed by
averaging from Eq. 1 at one second intervals across individuals. The product
of the transmission rate $\hat{\beta}$ and the duration time $T$ reflects the
likelihood of an individual being infected (see SM Section 2 for more
details).

Figure 4 plots the transmission likelihood $\hat{\beta}T$ against the
saturation for full-class vs. half-class and for not vaccinated vs. teacher
vaccinated scenarios. Saturation parameters appeared to grow with transmission
likelihood and, at the highest transmission likelihood, converged around
$75\%$ and $55\%$ for the not vaccinated and teacher vaccinated scenarios,
respectively. The area of each circle in Fig. 4 represents the classroom
density of the corresponding observation, suggesting a positive correlation
between density and saturation. The low density classes C and D (Full Class)
show similar infection patterns to the half-sized high density classes A and B
(Half Class). It nevertheless may be the case that simulated reductions in
classroom density slightly underestimate the effect of actual reduced density,
perhaps because simulations do not account for changes in actual social
interaction (see SM Section 4.F).

Figure 4: Scatter plot between the transmission likelihood $\hat{\beta}T$ and
saturation, where each of the 68 circles represents a full-class or a half-
class simulation of one of 34 observations in four classes for not vaccinated
and teacher vaccinated scenarios (See Table 1). The area of each circle is
proportional to the specific classroom density for a specific observation day,
and the color represents Classes A–D for the full and half-sized simulations,
respectively.

A mixed effects regression model, with observations nested in classes,
predicted saturation and transmission likelihood from class size and
vaccination (see Fig. 4). Half-sized scenarios were associated with a 18.2%
(Range: 16.0% - 20.3%) and 62.9% (Range: 49.72%. - 70.87%) mean reduction in
saturation and transmission likelihood compared to full-sized scenarios,
respectively. These full-size to half-size reductions in transmission
likelihood, $B=-.001,se=.00006,t=-17.73,p<.00001$, and in mean saturation
levels, $B=-.11,se=.01,t=-10.21,p<.00001$, were statistically significant.
Likewise, the teacher-vaccinated scenario was associated with a 25.3% (Range:
23.4% - 27.2%) and 29.4% (Range: 6.72% - 45.37%) mean reduction in saturation
and transmission likelihood compared to the not vaccinated scenario. Effects
of vaccination on lower transmission likelihood,
$B=-.0004,se=.00006,t=-6.79,p<.00001$, and mean saturation,
$B=-.13,se=.01,t=-12.59,p<0.00001$, were also significant. Of note, the
overall 25.3% reduction for the teacher-vaccinated scenario is a combination
of sources: 1) the 13.1% direct effect of reduced infection of vaccinated
teachers, and 2) the 12.2% indirect effect of vaccinated teachers infecting
fewer individuals (children and other teachers). All statistical analyses were
performed at the level of observations ($N=34$), which were nested in
classrooms. Consequently, the number of simulations performed did not
artificially reduce p-values. See SM Section 3 for additional details.

While the simulation uses the same infection parameters for each individual,
we observe transmission heterogeneity due to behavioral differences. For
example, different patient zeros lead to different infection patterns based on
individual differences in contact with others (i.e., variation in
$r_{1}\&\theta_{1}$) (see SM Section 4.G and Fig. S7 for more details).

### Policy Relevant Outcomes

To explore the criteria under which in-person schooling might be terminated
(school shut-down) after discovery of COVID-19 cases, we investigated the
timing of the emergence of infected, symptomatic individuals (children or
teachers). We first focused on the probability of a given set of simulations
yielding a first, second, or third symptomatic individuals (Fig. 5 and Table
S3). The $75\%$ symptomaticity rate used in simulations implies that the first
symptomatic case is patient zero in three of four simulations. The emergence
of the $1^{st}$, $2^{nd}$, and $3^{rd}$ symptomatic case was significantly
reduced in the half class scenarios, reflecting sensitivity to classroom
density (see Tables S3 and S4). The probability of not detecting a second
symptomatic individual, for example, was $35.2\%$ and $53.0\%$ for full and
half-sized classrooms, respectively. The corresponding probability was
$38.7\%$ and $49.5\%$ for not vaccinated and teacher vaccinated scenarios,
respectively. The reduction associated with teacher vaccination significantly
impacted the emergence of a $2^{nd}$, and $3^{rd}$ (but not a $1^{st}$)
symptomatic individual (see SM Section 3 and Table S4 for details).

Figure 5: The proportion of simulations in which the first, second, and third
symptomatic cases do not emerge in classes A-D. Teacher vaccinated and not
vaccinated as well as full and half-sized class scenarios are presented. The
horizontal axis indicates the nth symptomatic case. The vertical axes indicate
the proportion of simulations, and $\rho$’s in each legend represent the
classroom densities in the unit of people per $m^{2}$.

Next we plotted the time in days until the emergence of the first, second, and
third symptomatic individual in Fig. 6. The $75\%$ symptomaticity rate again
implied that the mean time to the onset of the first symptomatic case
typically corresponded to the $4$-day incubation period). Over all not
vaccinated classrooms, the median time to the emergence of the second
symptomatic individual was $8.51$ days in full class and $23.00$ days in half
class scenarios, reflecting sensitivity to classroom density (Fig. 6). The
median time to the emergence of the second symptomatic individual was $11.12$
days in the teacher vaccinated full-sized class scenarios and was not observed
in the teacher vaccinated half-sized class scenarios. The not-observed value
indicates that in over half the teacher vaccinated scenarios a second
symptomatic individual never emerged (see Table S5).

To explore the impact of classroom behavior, we categorized classroom time as
unstructured (free-play and transitions between activities) or as structured.
Structured activities were teacher-led and primarily occurred when children
were seated at tables (such as circle-time, shared book reading, meal-time,
and organized play). Simulations suggest that unstructured time was associated
with a higher trajectory of infection than unstructured time, presumably
because individuals were in closer proximity and were more mutually oriented
during these periods (see SM Section 4.B for more details). We also find that
transmission heterogeneity is naturally encoded in our model since individuals
behave differently and their social interactions are inhomogeneous (SM Section
4.G).

Figure 6: Time in days to the emergence of first, second, and third
symptomatic individuals in classes A–D for teacher vaccinated and not
vaccinated and for full and half-sized class scenarios. The horizontal axis
labels $n^{th}$ symptomatic case. The vertical axis indicates time in days.
Horizontal lines within each box indicate medians, and $\rho$’s in each legend
represent the classroom densities in the unit of people per $m^{2}$.

## Discussion

In the face of pandemic threats such as that posed by COVID-19, federal and
local governments must adopt policies that weigh the importance of educational
opportunities against the risks of infection within classroom settings. School
closure reduces transmission risk[45], but also limits the employment capacity
of health-care and other essential employees who are tasked with child
care[46], which reduces economic productivity[47]. Here we explore the impact
of less drastic prophylactic measures such as reducing classroom enrollment
and teacher vaccination.

Our goal was to describe the development of an agent-based SEIR model to
inform decisions and policy on curbing the impact of COVID-19 classroom
spread. Existing research examines infection as a product of expected rates of
interaction between infected and susceptible individuals, but does not examine
the role of individual (agents) in transmission. With few exceptions[19],
existing models do not describe the actual physical interactions through which
individuals infect others. Transmission, however, occurs in physical space
over time.

Since their original formulation[48], SEIR models have been used to model
infectious transmission at a population level[8, 9]. The current study’s RFID
system provided subsecond ($4Hz$) measures of the physical distance and
relative orientation of all individuals in a classroom[19]. Using the
resulting data we employed an agent-based model to examine infection in real-
time in a real-life system, the preschool classroom. The model uses actual
student and teacher interaction behavior as input and estimates infection
probabilities based on interpersonal distance and orientation.

The current policy-focused modeling scenarios indicated that classroom density
was a key parameter limiting SARS-CoV-2 transmission. Both simulations that
reduced classroom enrollment and simulations of teacher vaccination receipt
effectively reduced classroom density by decreasing the number of susceptible
individuals (see Fig. 4). The half class scenario reduced overall infection
proportion by $18.2\%$ while teacher vaccination was associated with a
$25.3\%$ reduction. In half class scenarios, the proportion of simulations in
which a second symptomatic individual never emerged rose from approximately
one third to one half of cases ($35.2\%$ to $53.0\%$) while teacher
vaccination was associated with a more modest rise ($38.7\%$ to $49.5\%$). The
timing of the emergence of these first, second, and third symptomatic cases –
canaries in the coal mine – was also delayed in both the half class and the
teacher vaccination scenarios.

A limitation of the current modeling approach is the need to make assumptions
about the physical and temporal parameters characterizing transmission.
Modeling parameters assumed that the time lag to infectiousness was one day
and the lag to symptomaticity was four days[34]. Likewise, the precise
physical parameters involved in SARS-CoV-2 transmission are not known[49]. The
current models, however, allow for flexible modeling of the physical and
temporal parameters associated with pathogen spread in enclosed spaces over
time. Note, for example, that CDC changes in safe distance recommendation from
$6$ to $3$ feet do not change key model parameters[27, 50]. Differential
infectiousness of SARS-CoV-2 variants[51] (e.g., the Delta variant,
$B.1.617.2$), changes in vaccine effectiveness and availability, and mask
mandates will all affect infection levels, as will community positivity. We
have explored the impact of parameters to the transmission patterns in SM
Section 4. That key contribution of the models is their flexibility. Indeed,
as our model has been purposely designed to be modular and therefore used as a
general framework, it allows inputs to change to reflect a multitude of future
scenarios that include different pathogens such as COVID variants and new
viral pathogens, and different transmission vectors including airborne and
aerosol particles, as well as variation in class room size/density, and
vaccination status.

## Data Availability

The datasets generated during and/or analysed during the current study are
available in the OSF repository,
https://osf.io/h7ks8/?view_only=5b03a51e79ff4c57b59c2a814d60dbd3.

## References

  * Zhu  _et al._ [2020] N. Zhu, D. Zhang, W. Wang, X. Li, B. Yang, J. Song, X. Zhao, B. Huang, W. Shi, R. Lu, _et al._ , A novel coronavirus from patients with pneumonia in china, 2019, New England Journal of Medicine 382, 727 (2020).
  * Layne  _et al._ [2020] S. P. Layne, J. M. Hyman, D. M. Morens, and J. K. Taubenberger, Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in wuhan, china: a modelling study, The Lancet 395, 689 (2020).
  * Dong  _et al._ [2020] E. Dong, H. Du, and L. Gardner, An interactive web-based dashboard to track COVID-19 in real time, The Lancet Infectious Diseases 20, 533 (2020).
  * wor [2021] COVID-19 weekly epidemiological update. World Health Organization, https://www.who.int/publications/m/item/weekly-epidemiological-update-on-covid-19---25-may-2021 (2021).
  * Giordano  _et al._ [2020] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A. Di Matteo, and M. Colaneri, Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy, Nature Medicine 26, 855 (2020).
  * cdc [2020a] Guidance for operating child care programs during covid-19. Centers for Disease Control and Prevention, https://www.cdc.gov/coronavirus/2019-ncov/science/science-briefs/scientific-brief-options-to-reduce-quarantine.html (2020a).
  * cdc [2021a] Operational strategy for k-12 schools through phased prevention. Centers for Disease Control and Prevention, https://www.cdc.gov/coronavirus/2019-ncov/community/schools-childcare/operation-strategy.html (2021a).
  * He  _et al._ [2020] S. He, Y. Peng, and K. Sun, SEIR modeling of the COVID-19 and its dynamics, Nonlinear Dynamics 101, 1667 (2020).
  * Yang  _et al._ [2020] Z. Yang, Z. Zeng, K. Wang, S.-S. Wong, W. Liang, M. Zanin, P. Liu, X. Cao, Z. Gao, Z. Mai, _et al._ , Modified SEIR and AI prediction of the epidemics trend of COVID-19 in china under public health interventions, Journal of Thoracic Disease 12, 165 (2020).
  * war [2020] A compendium of models that predict the spread of COVID-19. American Hospital Association, https://www.aha.org/guidesreports/2020-04-09-compendium-models-predict-spread-covid-19/ (2020).
  * Firth  _et al._ [2020] J. A. Firth, J. Hellewell, P. Klepac, S. Kissler, A. J. Kucharski, and L. G. Spurgin, Using a real-world network to model localized covid-19 control strategies, Nature medicine 26, 1616 (2020).
  * Nouvellet  _et al._ [2021] P. Nouvellet, S. Bhatia, A. Cori, K. E. Ainslie, M. Baguelin, S. Bhatt, A. Boonyasiri, N. F. Brazeau, L. Cattarino, L. V. Cooper, _et al._ , Reduction in mobility and covid-19 transmission, Nature communications 12, 1 (2021).
  * Gozzi  _et al._ [2021] N. Gozzi, M. Tizzoni, M. Chinazzi, L. Ferres, A. Vespignani, and N. Perra, Estimating the effect of social inequalities on the mitigation of covid-19 across communities in santiago de chile, Nature communications 12, 1 (2021).
  * Müller  _et al._ [2021] S. A. Müller, M. Balmer, W. Charlton, R. Ewert, A. Neumann, C. Rakow, T. Schlenther, and K. Nagel, Predicting the effects of covid-19 related interventions in urban settings by combining activity-based modelling, agent-based simulation, and mobile phone data, medRxiv (2021).
  * Chang  _et al._ [2021] S. Chang, E. Pierson, P. W. Koh, J. Gerardin, B. Redbird, D. Grusky, and J. Leskovec, Mobility network models of covid-19 explain inequities and inform reopening, Nature 589, 82 (2021).
  * Aleta  _et al._ [2020] A. Aleta, D. Martin-Corral, A. P. y Piontti, M. Ajelli, M. Litvinova, M. Chinazzi, N. E. Dean, M. E. Halloran, I. M. Longini Jr, S. Merler, _et al._ , Modelling the impact of testing, contact tracing and household quarantine on second waves of covid-19, Nature Human Behaviour 4, 964 (2020).
  * Kraemer  _et al._ [2020] M. U. Kraemer, C.-H. Yang, B. Gutierrez, C.-H. Wu, B. Klein, D. M. Pigott, O. C.-. D. W. Group†, L. du Plessis, N. R. Faria, R. Li, _et al._ , The effect of human mobility and control measures on the covid-19 epidemic in china, Science 368, 493 (2020).
  * Tupper  _et al._ [2020] P. Tupper, H. Boury, M. Yerlanov, and C. Colijn, Event-specific interventions to minimize covid-19 transmission, Proceedings of the National Academy of Sciences 117, 32038 (2020).
  * Stehlé  _et al._ [2011] J. Stehlé, N. Voirin, A. Barrat, C. Cattuto, V. Colizza, L. Isella, C. Régis, J.-F. Pinton, N. Khanafer, W. Van den Broeck, _et al._ , Simulation of an SEIR infectious disease model on the dynamic contact network of conference attendees, BMC Medicine 9, 1 (2011).
  * Smieszek  _et al._ [2019] T. Smieszek, G. Lazzari, and M. Salathé, Assessing the dynamics and control of droplet-and aerosol-transmitted influenza using an indoor positioning system, Scientific reports 9, 1 (2019).
  * Lopez  _et al._ [2020] A. S. Lopez, M. Hill, J. Antezano, D. Vilven, T. Rutner, L. Bogdanow, C. Claflin, I. T. Kracalik, V. L. Fields, A. Dunn, _et al._ , Transmission dynamics of COVID-19 outbreaks associated with child care facilities–Salt Lake City, Utah, April–July 2020, Morbidity and Mortality Weekly Report 69, 1319 (2020).
  * Laxminarayan  _et al._ [2020] R. Laxminarayan, B. Wahl, S. R. Dudala, K. Gopal, S. Neelima, K. J. Reddy, J. Radhakrishnan, J. A. Lewnard, _et al._ , Epidemiology and transmission dynamics of COVID-19 in two indian states, Science 370, 691 (2020).
  * Bowne  _et al._ [2017] J. B. Bowne, K. A. Magnuson, H. S. Schindler, G. J. Duncan, and H. Yoshikawa, A meta-analysis of class sizes and ratios in early childhood education programs: Are thresholds of quality associated with greater impacts on cognitive, achievement, and socioemotional outcomes?, Educational Evaluation and Policy Analysis 39, 407 (2017).
  * Irvin  _et al._ [2018] D. W. Irvin, S. A. Crutchfield, C. R. Greenwood, W. D. Kearns, and J. Buzhardt, An automated approach to measuring child movement and location in the early childhood classroom, Behavior Research Methods 50, 890 (2018).
  * Tang and Settles [2008] J. W. Tang and G. S. Settles, Coughing and aerosols, New England Journal of Medicine 359, e19 (2008).
  * Bazant and Bush [2021] M. Z. Bazant and J. W. Bush, A guideline to limit indoor airborne transmission of covid-19, Proceedings of the National Academy of Sciences 118 (2021).
  * Bertozzi  _et al._ [2020] A. L. Bertozzi, E. Franco, G. Mohler, M. B. Short, and D. Sledge, The challenges of modeling and forecasting the spread of covid-19, Proceedings of the National Academy of Sciences 117, 16732 (2020).
  * Li  _et al._ [2020a] Q. Li, X. Guan, P. Wu, X. Wang, L. Zhou, Y. Tong, R. Ren, K. S. Leung, E. H. Lau, J. Y. Wong, _et al._ , Early transmission dynamics in Wuhan, China, of novel coronavirus–infected pneumonia, New England Journal of Medicine 382, 1199 (2020a).
  * Kucharski  _et al._ [2020] A. J. Kucharski, T. W. Russell, C. Diamond, Y. Liu, J. Edmunds, S. Funk, R. M. Eggo, F. Sun, M. Jit, J. D. Munday, _et al._ , Early dynamics of transmission and control of COVID-19: a mathematical modelling study, The Lancet Infectious Diseases 20, 553 (2020).
  * Abbott  _et al._ [2020] S. Abbott, J. Hellewell, J. Munday, S. Funk, and CMMID nCoV working group, The transmissibility of novel coronavirus in the early stages of the 2019-20 outbreak in wuhan: Exploring initial point-source exposure sizes and durations using scenario analysis, Wellcome Open Research 5 (2020).
  * Read  _et al._ [2020] J. M. Read, J. R. Bridgen, D. A. Cummings, A. Ho, and C. P. Jewell, Novel coronavirus 2019-nCoV: early estimation of epidemiological parameters and epidemic predictions, Preprint https://doi.org/10.1101/2020.01.23.20018549 (2020).
  * Liu  _et al._ [2020] T. Liu, J. Hu, M. Kang, L. Lin, H. Zhong, J. Xiao, G. He, T. Song, Q. Huang, Z. Rong, _et al._ , Transmission dynamics of 2019 novel coronavirus (2019-nCoV), Preprint http://dx.doi.org/10.2139/ssrn.3526307 (2020).
  * Ferretti  _et al._ [2020] L. Ferretti, C. Wymant, M. Kendall, L. Zhao, A. Nurtay, L. Abeler-Dörner, M. Parker, D. Bonsall, and C. Fraser, Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing, Science 368 (2020).
  * cdc [2021b] Interim clinical guidance for management of patients with confirmed coronavirus disease (covid-19). Centers for Disease Control and Prevention, https://www.cdc.gov/coronavirus/2019-ncov/hcp/clinical-guidance-management-patients.html (2021b).
  * Mossong  _et al._ [2008] J. Mossong, N. Hens, M. Jit, P. Beutels, K. Auranen, R. Mikolajczyk, M. Massari, S. Salmaso, G. S. Tomba, J. Wallinga, _et al._ , Social contacts and mixing patterns relevant to the spread of infectious diseases, PLoS Medicine 5, e74 (2008).
  * cdc [2021c] COVID-19 contact tracing. Centers for Disease Control and Prevention, https://www.cdc.gov/coronavirus/2019-ncov/daily-life-coping/contact-tracing.html (2021c).
  * Johansson  _et al._ [2021] M. A. Johansson, T. M. Quandelacy, S. Kada, P. V. Prasad, M. Steele, J. T. Brooks, R. B. Slayton, M. Biggerstaff, and J. C. Butler, Sars-cov-2 transmission from people without COVID-19 symptoms, JAMA Network Open 4, e2035057 (2021).
  * Lauer  _et al._ [2020] S. A. Lauer, K. H. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. R. Meredith, A. S. Azman, N. G. Reich, and J. Lessler, The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application, Annals of Internal Medicine 172, 577 (2020).
  * Guan  _et al._ [2020] W.-j. Guan, Z.-y. Ni, Y. Hu, W.-h. Liang, C.-q. Ou, J.-x. He, L. Liu, H. Shan, C.-l. Lei, D. S. Hui, B. Du, L.-j. Li, G. Zeng, K.-Y. Yuen, R.-c. Chen, C.-l. Tang, T. Wang, P.-y. Chen, J. Xiang, S.-y. Li, J.-l. Wang, Z.-j. Liang, Y.-x. Peng, L. Wei, Y. Liu, Y.-h. Hu, P. Peng, J.-m. Wang, J.-y. Liu, Z. Chen, G. Li, Z.-j. Zheng, S.-q. Qiu, J. Luo, C.-j. Ye, S.-y. Zhu, and N.-s. Zhong, Clinical characteristics of coronavirus disease 2019 in China, New England Journal of Medicine 382, 1708 (2020).
  * cdc [2020b] Science brief: Options to reduce quarantine for contacts of persons with SARS-CoV-2 infection using symptom monitoring and diagnostic testing. Centers for Disease Control and Prevention, https://www.cdc.gov/coronavirus/2019-ncov/science/science-briefs/scientific-brief-options-to-reduce-quarantine.html (2020b).
  * [41] Interim clinical considerations for use of COVID-19 vaccines currently authorized in the united states. Centers for Disease Control and Prevention, https://www.cdc.gov/vaccines/covid-19/clinical-considerations/covid-19-vaccines-us.html.
  * pfi [2021] Real-world evidence confirms high effectiveness of pfizer-biontech COVID-19 vaccine and profound public health impact of vaccination one year after pandemic declared. Pfizer, https://www.pfizer.com/news/press-release/press-release-detail/real-world-evidence-confirms-high-effectiveness-pfizer (2021).
  * cdc [2021d] Moderna COVID-19 vaccine overview and safety. Centers for Disease Control and Prevention, https://www.cdc.gov/coronavirus/2019-ncov/vaccines/different-vaccines/Moderna.html (2021d).
  * cdc [2021e] Johnson & johnson’s janssen COVID-19 vaccine overview and safety. Centers for Disease Control and Prevention, https://www.cdc.gov/coronavirus/2019-ncov/vaccines/different-vaccines/janssen.html (2021e).
  * Honein  _et al._ [2021] M. A. Honein, L. C. Barrios, and J. T. Brooks, Data and policy to guide opening schools safely to limit the spread of SARS-CoV-2 infection, JAMA 325, 823 (2021).
  * O’Sullivan  _et al._ [2009] T. L. O’Sullivan, C. Amaratunga, K. P. Phillips, W. Corneil, E. O’Connor, L. Lemyre, and D. Dow, If schools are closed, who will watch our kids? family caregiving and other sources of role conflict among nurses during large-scale outbreaks, Prehospital and Disaster Medicine 24, 321 (2009).
  * Li  _et al._ [2020b] X. Li, W. Xu, M. Dozier, Y. He, A. Kirolos, E. Theodoratou, _et al._ , The role of children in transmission of SARS-CoV-2: A rapid review, Journal of Global Health 10, 011101 (2020b).
  * Kermack and McKendrick [1927] W. O. Kermack and A. G. McKendrick, A contribution to the mathematical theory of epidemics, Proceedings of the Royal Society of London. 115, 700 (1927).
  * Rosti  _et al._ [2020] M. Rosti, S. Olivieri, M. Cavaiola, A. Seminara, and A. Mazzino, Fluid dynamics of covid-19 airborne infection suggests urgent data for a scientific design of social distancing, Scientific reports 10, 1 (2020).
  * van den Berg  _et al._ [2021] P. van den Berg, E. M. Schechter-Perkins, R. S. Jack, I. Epshtein, R. Nelson, E. Oster, and W. Branch-Elliman, Effectiveness of three versus six feet of physical distancing for controlling spread of COVID-19 among primary and secondary students and staff: A retrospective, state-wide cohort study, Clinical Infectious Diseases , ciab230 (2021).
  * Lam-Hine  _et al._ [2021] T. Lam-Hine, S. McCurdy, L. Santora, _et al._ , Outbreak associated with SARS-CoV-2 B.1.617.2 (Delta) variant in an elementary school — marin county, california, may–june 2021. Centers for Disease Control and Prevention (2021).

## Acknowledgements

This work was financially supported by the National Science Foundation,
IBSS-L1620294; Institute of Education Sciences, R324A180203; Microsoft AI for
Health COVID-19 Grant Program and Google Cloud COVID-19 Research Credits
Program.

## Author contributions statement

C.S. and D.M. defined the overall problem and designed the study. P.W., D.M.,
L.P., Y.T. and C.S. conducted the literature review. D.M., L.P., Y.T., Y.Z.
collected and preprocessed the empirical data. P.W., Y.Z. and C.S. developed
the transmission model. Y.Z., Y.T. and M.S. prepared and conducted numerical
simulations. L.P., Y.Z. and Y.T. performed the statistical analysis and
prepared the figures. P.W., C.S., and D.M. interpreted findings and wrote the
manuscript. All authors provided critical discussion and approved its
submission.

## Additional information

Accession codes: The codes to run the simulations based on the datasets in
this study are available at
https://anonymous.4open.science/r/Classroom_COVID_Simluation-2E38.

The authors declare no competing interests.