# Correction to the photometric colors of Gaia Early Data Release 3

Zexi Niu  National Astronomical Observatories, Chinese Academy of Sciences
20A Datun Road, Chaoyang District, Beijing, China  Haibo Yuan  Department of
Astronomy, Beijing Normal University
19th Xinjiekouwai Street, Haidian District, Beijing, China  Jifeng Liu
National Astronomical Observatories, Chinese Academy of Sciences
20A Datun Road, Chaoyang District, Beijing, China

###### Abstract

In this work, we use the spectroscopy-based stellar color regression (SCR)
method with $\sim$ 0.7 million common stars between LAMOST DR7 and Gaia EDR3
to acquire color corrections in $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm RP}$. A
sub-mmag precision is achieved. Our results demonstrate that improvements in
the calibration process of the EDR3 have removed the color term in $G_{\rm
BP}-G_{\rm RP}$ and eliminated the discontinuity caused by the changes of
instrument configurations to a great extent. However, modest systematic trends
with $G$ magnitude are still detected. The corresponding color correction
terms as a function of $G$ are provided for $9.5<G<17.5$ mag and compared with
other determinations. We conclude that the corrections given in this work are
particularly suited for cases where the color-color investigations are
required while for color-magnitude investigations other corrections may be
better due to systematic with reddening. Possible applications of our results
are discussed.

Astronomy data analysis, Fundamental parameters of stars, Stellar photometry

††journal: ApJL

##  1 Introduction

Very recently, Gaia Collaboration published the Early Data Release 3 (EDR3)
(Gaia Collaboration et al., 2020) based on the first 34 months of its nominal
mission (Gaia Collaboration et al., 2016), providing $G$ band photometry for
1.8 billion sources brighter than 21 and 1.5 billion sources with $G_{BP}$ and
$G_{RP}$ photometry, with an uniform calibration at mmag level. In the Gaia
DR2 era, several works have detected the magnitude dependent systematic errors
up to 10 mmag or higher (Maíz Apellániz & Weiler, 2018; Weiler, 2018;
Casagrande & VandenBerg, 2018; Niu et al., 2021) and the discontinuities
caused by the changes of instrument configurations (Evans et al., 2018; Niu et
al., 2021). Among them, using about 0.5 million well selected common stars
between the LAMOST DR5 (Luo et al., 2015; Zhao et al., 2012) and Gaia DR2, Niu
et al. (2021) (hereafter Paper I) applied the stellar color regression (SCR)
method (Yuan et al., 2015) to calibrate the $G-G_{\rm RP}$ and $G_{\rm
BP}-G_{\rm RP}$ colors. With an unprecedented precision of about 1 mmag,
systematic trends with $G$ magnitude are revealed in great detail for both
colors, reflecting changes in instrument configurations. Color-dependent
trends are also found for the $G_{\rm BP}-G_{\rm RP}$ for stars brighter than
$G\sim$ 11.5 mag.

From DR2 (Gaia Collaboration et al., 2018) to EDR3, a number of important
improvements have been implemented to further reduce its photometric (random
and systematic) errors, including the fitting of $G$ fluxes, the processing of
BP and RP spectra, and the calibration processes (Riello et al., 2020;
Fabricius et al., 2020). The median uncertainties of $G$ magnitudes are
reduced by almost a factor of two, reaching 0.2 mmag at $10<G<14$, 0.8 mmag at
$G\sim$ 17, and 2.6 mmag at $G\sim$ 19\. In order to take full advantage of
its exquisite photometric quality, in this work, we follow the same routine of
Paper I to validate and correct possible magnitude/color dependent systematics
in EDR3.

The paper is organized as follows. We briefly describe our data and method in
Section 2, with differences from Paper I stressed. The result is presented and
discussed in Section 3. We summarize in Section 4.

##  2 Data and Method

We use the same method and follow the same steps as in Paper I, which are
briefly summarized below. Details of criteria used to select samples and the
SCR method can be found in Paper I.

We combine the newest Gaia EDR3 with the LAMOST DR7 (Luo et al., 2015) and
apply the same constraints as in Paper I, e.g., E(B$-$V) $<$ 0.05 mag
according to the Schlegel, Finkbeiner, & Davis (1998, hereafter SFD) dust
reddening map, galactic latitude $|b|>20$ deg, and vertical distance to the
galactic disk $|Z|>0.2$ kpc. In Figure 3 of Paper I, we clarified that the
signal-to-noise ratio for the $g$ band ($S/N_{\rm g}$) of 20 is sufficient for
this work because there are no systematic effects in the LAMOST stellar
parameters with the $S/N_{\rm g}$ till at very low $S/N_{\rm g}$ of about 15.
So we adopt a lower cut on $S/N_{\rm g}$ of 20 to generate a larger sample. At
$S/N_{\rm g}$ $>$ 20, the agreement of the LAMOST parameters with APOGEE is
better than 120 K for $T_{\rm eff}$, 0.15 dex for $Log$ g, and 0.1 dex for
[Fe/H], as shown in Figure 3 of Paper I.

We finally select a sample containing 779,691 main-sequence (MS) stars and
71,952 Red Giant Branch (RGB) stars, covering $9.5<G<17.5$ magnitude range.
Note that as proposed by the Gaia Collaboration et al. (2020) and Riello et
al. (2020), the corrected $G$ magnitudes for sources with 6-parameter
astrometric solutions are used in this paper.

Due to the different passbands between DR2 and EDR3, their reddening
coefficients are slightly different. Following Sun et al. (to be submitted),
we have empirically determined the temperature and reddening dependent
reddening coefficients for the EDR3 colors, as given by the Equations 1 and 2.
Note that the reddening corrections in this work do not take into account the
distance of sources and assume that all sources to be beyond than the source
of reddening. All colors referred to hereafter are dereddened using the SFD
map and the empirical coefficients.

$\displaystyle R(G_{\rm BP}-G)=1.428-0.539\times E(B-V)_{\rm SFD}+0.406$ (1)
$\displaystyle\times{E(B-V)}^{2}_{\rm SFD}-1.976\times 10^{-4}\times T_{\rm
eff}+2.004\times 10^{-5}$ $\displaystyle\times T_{\rm eff}\times E(B-V)_{\rm
SFD}-1.187\times 10^{-8}\times{T_{\rm eff}}^{2}$

$\displaystyle R(G_{\rm BP}-G_{\rm RP})=1.684-0.839\times E(B-V)_{\rm
SFD}+0.555$ (2) $\displaystyle\times{E(B-V)}^{2}_{\rm SFD}-1.390\times
10^{-4}\times T_{\rm eff}-3.142\times 10^{-6}$ $\displaystyle\times T_{\rm
eff}\times E(B-V)_{\rm SFD}+1.520\times 10^{-8}\times{T_{\rm eff}}^{2}$

As in Paper I, a control sample of $13.3<G<13.7$ are selected to define the
empirical relations between the Gaia intrinsic colors and the LAMOST DR7
stellar parameters. Then we divide the selected samples into six subsamples of
different evolution stages and colors. MS stars are divided into four
subsamples with the median $G_{\rm BP}-G_{\rm RP}$ colors of 1.12, 0.92, 0.76,
and 0.64 mag, respectively. RGB stars are divided into two subsamples, with
the median $G_{\rm BP}-G_{\rm RP}$ colors of 1.11 and 0.99 mag, respectively.
Applying these relations to a given (sub-)sample, the median values of the
color residuals between the SCR-derived colors and the corresponding Gaia EDR3
ones as a function of $G$ are obtained as the correction terms, as shown in
Figures LABEL:fig1 and LABEL:fig2. A 3$\sigma$-clipping is performed in the
process. Specifically, corrections can be expressed as:

$C^{{}^{\prime}}=C+\Delta C$ (3)

where $C^{{}^{\prime}}$ is the corrected color, $C$ is the Gaia EDR3 color and
$\Delta C$ is the color correction term.

##  3 Result and Discussion

Figure 3 plots color correction curves yielded by different subsamples. As
already demonstrated in Paper I: (1) the correction terms are independent of
the stellar evolution stages; (2) the inconsistency between the red and blue
subsamples when $G>14$ mag is due to the selection function of the LAMOST data
and the spatially dependent systematics of the SFD reddening map. Therefore,
we compute the recommended correction curves using the three blue MS
subsamples (MS 0.92, MS 0.76, and MS 0.64). The results are over-plotted in
black lines in Figure 3 and listed in Table 1. The standard deviations of the
difference between the recommended curve and the three blue MS subsamples for
$G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm RP}$ are 0.3 and 0.7 mmag, respectively.
This suggests that the calibration curve of each subsample has a typical
random error smaller than 1.0 mmag, and the random errors of the recommended
curves are even smaller. In the top panel of Figure 3, the systematic trend of
$G-G_{\rm RP}$ is smaller than the one found for Gaia DR2, with a maximum
range of less than 10 mmag in total. Feature coming from the observation modes
changing at $G\sim 16$ mag is gone and those at $G\sim$ 11 and 13 mag are
visible but at a level of about 1 mmag from the recommended curve. This
confirms the big improvement in the $G$ band photometry. The bottom panel
shows that the $G_{\rm BP}-G_{\rm RP}$ correction curves are no longer color-
dependent at the bright end, at least for F/G/K stars. However, the magnitude-
dependent trend is still significant, as in DR2.

Figure 3: Color correction curves yielded by different subsamples and the
recommended ones for $G-G_{\rm RP}$ (top) and $G_{\rm BP}-G_{\rm RP}$
(bottom).

To test our color correction curves, we select a MS sample within a narrow
[Fe/H] range of $-0.5<{\rm Fe/H}<-0.25$, and compare their distributions in
the color-color diagram before and after corrections using the recommended
curves in Table 1. The results are shown in the left panels of Figure 4.
Magnitude-dependent offsets in the color-color diagram are clearly seen from
the top left panel. After corrections, there are no offsets, and the total
width becomes narrower as expected. Furthermore, like Figure 32 in Fabricius
et al. (2020) of EDR3 and Figure 31 in Arenou et al. (2018) of DR2, we plot a
2D histogram of the $G-G_{\rm RP}$ residuals in the right panels of Figure 4
with respect to the metallicity-dependent stellar color locus. The improved
yet still subsistent trend with magnitude is clearly seen in the top panel,
especially at the bright and faint ends. After corrections of this work, the
trend becomes flat and centered at zero, demonstrating the power of our
corrections in the color-color diagram.

Figure 4: From top to bottom: before and after color corrections. left: color-
color diagram of a MS sample within a narrow [Fe/H] range. right: Residuals
from a global $G-G_{\rm RP}=f(G_{\rm BP}-G_{\rm RP},{\rm[Fe/H]})$ relation.

There are also color corrections deduced from other approaches. For example,
using 10,000 well selected Landolt standard stars (Clem & Landolt, 2013), Yang
et al. (2021) has obtained magnitude and color corrections of the Gaia EDR3 by
training the observed $UBVRI$ magnitudes into the Gaia EDR3 magnitudes and
colors. Their results (red lines in their Figure 4) are plotted in green lines
in Figure 5. Results from synthetic colors of the CALSPEC (Bohlin 2014)
spectra are also plotted in black dots. Only 61 stars with
phot$\\_$bp$\\_$rp$\\_$excess$\\_$factor $<0.1$ after Riello et al. (2020)
correction and $G>8$ mag are shown. Note that Yang et al. (2021) adopted a
control sample of $17<G<17.5$ mag, their correction curves were shifted to
match the results of the CALSPEC spctera. In this work, we have adopted a
different control sample of $13.2<G<13.6$ mag, therefore, a constant offset
may well exist. To make a straightforward comparison, the recommended curves
in this work, which are over-plotted in red solid lines, are systematically
shifted by a few mmag to match the Yang et al. (2021) corrections at 14 mag.
Given the small number of black dots and their internal scatter (about 10 mmag
for $G-G_{\rm RP}$ and 15 mmag for $G_{\rm BP}-G_{\rm RP}$), both results from
this work and Yang et al. (2021) are consistent with the black dots. However,
discrepancies between this work and Yang et al. (2021) are up to about $\pm$5
mmag for $G-G_{\rm RP}$ and $\pm$10 mmag for $G_{\rm BP}-G_{\rm RP}$, and are
correlated with the $G$ magnitude.

To investigate the origins and effects of the discrepancies between this work
and Yang et al. (2021), their differences of the color corrections without
shifting in $G_{\rm BP}-G_{\rm RP}$ are plotted against those in $G-G_{\rm
RP}$ in the bottom panel of Figure 5. A strong linear correlation is found.
The slope is in good agreement with the median value of $\frac{R(G-G_{\rm
RP})}{R(G_{\rm BP}-G_{RP})}$. Considering that reddening corrections are not
involved in Yang et al. (2021), the result suggests that the discrepancies are
mainly caused by imperfect reddening corrections in this work. We use the 2D
SFD reddening map. Despite the systematic errors that depend on spatial
position and dust temperature (e.g., Peek & Graves 2010; Sun et al. to be
submitted) in the SFD reddening map as discussed in Paper I, the map also
tends to overestimate reddening for stars that are within the Galactic dust
layer. Although in this work (and Paper I) we require stars of $|Z|>0.2$ kpc,
their reddening corrections may be overestimated to some extent around 0.01
mag in E(B-V), depending on their distances/magnitudes. It is not surprising,
as there are increasing evidences supporting the co-existence of a thin dust
disk and a thick dust disk in the Galaxy (e.g., Yuan et al. in prep.; Guo et
al. 2021; Zhang et al., to be submitted). The scale height of the thick dust
disk is about 200 - 400 pc in the solar neighborhood. Dust clouds in the
Galactic halo (Yuan et al. in prep.) and halos of other galaxies (e.g., M 31
and M 33, Zhang & Yuan 2020) are also detected. We redo the procedure with an
excessive sample of $|Z|>1$ kpc and a control sample of $15.3<G<15.5$ mag. Its
results are plotted in red dashed lines, which match the green lines much
better, consistent with the above scenario. It suggests that our correction
curves of individual colors are subject to systematic errors from reddening
correction using the SFD map. We further test other 2D reddening maps of Lenz
et al. (2017), Peek & Graves (2010), and Planck Collaboration et al. (2014),
and find no big differences with the SFD map.

Figure 5: Top and middle: Comparisons of different color correction curves.
Red solid line: $|Z|>0.2$ kpc using the SCR method. Red dashed line: $|Z|>1$
kpc using the SCR method. Black dots: from the CALSPEC spectra. Green line:
from Yang et al. (2021). Both red solid and dashed lines are systematically
shifted by a few mmag to match the green lines at 14 mag for comparison.

bottom: Differences of the color correction terms in $G_{\rm BP}-G_{\rm RP}$
versus those in $G-G_{\rm RP}$. The two dashed lines have the same slope,
determined by the median value of $\frac{R(G-G_{\rm RP})}{R(G_{\rm
BP}-G_{RP})}$.

Fortunately, even though both correction curves of the two Gaia colors
suffering systematic errors in reddening correction, their systematic errors
are largely canceled in the color-color diagram, as the reddening vector is
almost parallel to the stellar locus in the color-color diagram. We estimate
the standard deviation of the differences between the purple points and their
corresponding gray line, which is 0.3 mmag and contributed by the random
errors of corrections from Yang et al. (2021) and this work. Given that the
typical errors of color correction curves of Yang et al. (2021) are about 0.2
– 0.4 mmag, the small standard deviation of 0.3 mmag suggests that the typical
uncertainties in this work are much smaller than 0.3 mmag. Moreover, with 0.7
million stars, our sample enables a much higher resolution in $G$ magnitude,
yielding improved corrections in the color-color diagram compared to Yang et
al. (2021) (See their Figure 5).

It is worth clarifying that our color corrections in this work and in Paper I
are precise to sub-mmag only in the cases where a color-color diagram is used.
In other cases, for example where a color-magnitude diagram is used, the
corrections from Yang et al. (2021) are preferred. Another thing should be
mentioned is that we compute the curves within approximately $9.5<G<17.5$ and
$0.6<G_{\rm BP}-G_{\rm RP}<1.3$, so for stars outside the above ranges, e.g.,
very blue stars brighter than G $\sim$ 13 (Riello et al., 2020), the curves
should be used with caution. Nevertheless, this corrected color-color diagram
could be helpful in a number of studies, for example, determining reliable
photometric metallicities for an enormous and magnitude-limited sample of
stars from Gaia (Xu et al. to be submitted) and estimating binary fractions of
a volume-limited sample of stars (Niu et al. to be submitted).

Table 1: Color correction curves. The first column is $G$ magnitude. The second and third columns are respectively recommended $G-G_{\rm RP}$ and $G_{\rm BP}-G_{\rm RP}$ calibration terms in units of mmag.  G | recom | G | recom | G | recom | G | recom | G | recom | G | recom
---|---|---|---|---|---|---|---|---|---|---|---
9.50 |  $-$ 2.54 | 2.95 | 10.85 | 0.17 | 6.00 | 12.20 | 1.68 | 3.84 | 13.55 |  $-$ 0.04 |  $-$ 0.02 | 14.90 |  $-$ 4.65 |  $-$ 7.07 | 16.25 |  $-$ 5.41 |  $-$ 10.68
9.53 |  $-$ 2.48 | 3.01 | 10.88 | 0.30 | 6.19 | 12.23 | 1.68 | 3.68 | 13.58 |  $-$ 0.04 |  $-$ 0.10 | 14.93 |  $-$ 4.74 |  $-$ 7.19 | 16.28 |  $-$ 5.38 |  $-$ 10.75
9.56 |  $-$ 2.39 | 3.10 | 10.91 | 0.43 | 6.30 | 12.26 | 1.69 | 3.55 | 13.61 |  $-$ 0.06 |  $-$ 0.21 | 14.96 |  $-$ 4.81 |  $-$ 7.26 | 16.31 |  $-$ 5.33 |  $-$ 10.90
9.59 |  $-$ 2.33 | 3.16 | 10.94 | 0.60 | 6.36 | 12.29 | 1.68 | 3.39 | 13.64 |  $-$ 0.07 |  $-$ 0.32 | 14.99 |  $-$ 4.90 |  $-$ 7.28 | 16.34 |  $-$ 5.28 |  $-$ 11.09
9.62 |  $-$ 2.26 | 3.24 | 10.97 | 0.74 | 6.38 | 12.32 | 1.66 | 3.29 | 13.67 |  $-$ 0.09 |  $-$ 0.41 | 15.02 |  $-$ 4.96 |  $-$ 7.34 | 16.37 |  $-$ 5.24 |  $-$ 11.28
9.65 |  $-$ 2.18 | 3.32 | 11.00 | 0.87 | 6.42 | 12.35 | 1.64 | 3.23 | 13.70 |  $-$ 0.11 |  $-$ 0.47 | 15.05 |  $-$ 5.01 |  $-$ 7.46 | 16.40 |  $-$ 5.17 |  $-$ 11.53
9.68 |  $-$ 2.12 | 3.38 | 11.03 | 1.00 | 6.79 | 12.38 | 1.62 | 3.23 | 13.73 |  $-$ 0.14 |  $-$ 0.54 | 15.08 |  $-$ 5.07 |  $-$ 7.64 | 16.43 |  $-$ 5.12 |  $-$ 11.66
9.71 |  $-$ 2.04 | 3.46 | 11.06 | 1.15 | 6.66 | 12.41 | 1.59 | 3.20 | 13.76 |  $-$ 0.18 |  $-$ 0.66 | 15.11 |  $-$ 5.12 |  $-$ 7.83 | 16.46 |  $-$ 5.07 |  $-$ 11.71
9.74 |  $-$ 1.98 | 3.53 | 11.09 | 1.25 | 6.50 | 12.44 | 1.56 | 3.03 | 13.79 |  $-$ 0.21 |  $-$ 0.78 | 15.14 |  $-$ 5.17 |  $-$ 7.97 | 16.49 |  $-$ 5.00 |  $-$ 11.75
9.77 |  $-$ 1.92 | 3.60 | 11.12 | 1.34 | 6.32 | 12.47 | 1.53 | 2.86 | 13.82 |  $-$ 0.26 |  $-$ 0.94 | 15.17 |  $-$ 5.22 |  $-$ 8.13 | 16.52 |  $-$ 4.94 |  $-$ 11.77
9.80 |  $-$ 1.84 | 3.68 | 11.15 | 1.41 | 6.20 | 12.50 | 1.50 | 2.71 | 13.85 |  $-$ 0.32 |  $-$ 1.09 | 15.20 |  $-$ 5.26 |  $-$ 8.21 | 16.55 |  $-$ 4.89 |  $-$ 11.80
9.83 |  $-$ 1.77 | 3.75 | 11.18 | 1.46 | 6.15 | 12.53 | 1.47 | 2.55 | 13.88 |  $-$ 0.38 |  $-$ 1.20 | 15.23 |  $-$ 5.30 |  $-$ 8.28 | 16.58 |  $-$ 4.82 |  $-$ 11.85
9.86 |  $-$ 1.69 | 3.83 | 11.21 | 1.50 | 6.21 | 12.56 | 1.45 | 2.46 | 13.91 |  $-$ 0.46 |  $-$ 1.31 | 15.26 |  $-$ 5.35 |  $-$ 8.41 | 16.61 |  $-$ 4.77 |  $-$ 11.87
9.89 |  $-$ 1.60 | 3.90 | 11.24 | 1.51 | 6.27 | 12.59 | 1.44 | 2.41 | 13.94 |  $-$ 0.53 |  $-$ 1.38 | 15.29 |  $-$ 5.39 |  $-$ 8.55 | 16.64 |  $-$ 4.69 |  $-$ 11.85
9.92 |  $-$ 1.52 | 3.93 | 11.27 | 1.52 | 6.33 | 12.62 | 1.42 | 2.33 | 13.97 |  $-$ 0.60 |  $-$ 1.43 | 15.32 |  $-$ 5.43 |  $-$ 8.73 | 16.67 |  $-$ 4.63 |  $-$ 11.85
9.95 |  $-$ 1.43 | 4.00 | 11.30 | 1.52 | 6.33 | 12.65 | 1.40 | 2.25 | 14.00 |  $-$ 0.72 |  $-$ 1.48 | 15.35 |  $-$ 5.48 |  $-$ 8.92 | 16.70 |  $-$ 4.57 |  $-$ 11.92
9.98 |  $-$ 1.38 | 4.06 | 11.33 | 1.52 | 6.20 | 12.68 | 1.39 | 2.15 | 14.03 |  $-$ 0.81 |  $-$ 1.52 | 15.38 |  $-$ 5.50 |  $-$ 9.04 | 16.73 |  $-$ 4.50 |  $-$ 12.05
10.01 |  $-$ 1.35 | 4.13 | 11.36 | 1.50 | 5.97 | 12.71 | 1.38 | 2.06 | 14.06 |  $-$ 0.90 |  $-$ 1.56 | 15.41 |  $-$ 5.54 |  $-$ 9.10 | 16.76 |  $-$ 4.44 |  $-$ 12.24
10.04 |  $-$ 1.31 | 4.23 | 11.39 | 1.48 | 5.76 | 12.74 | 1.36 | 1.96 | 14.09 |  $-$ 1.02 |  $-$ 1.68 | 15.44 |  $-$ 5.56 |  $-$ 9.17 | 16.79 |  $-$ 4.38 |  $-$ 12.51
10.07 |  $-$ 1.26 | 4.28 | 11.42 | 1.47 | 5.61 | 12.77 | 1.33 | 1.89 | 14.12 |  $-$ 1.14 |  $-$ 1.86 | 15.47 |  $-$ 5.57 |  $-$ 9.30 | 16.82 |  $-$ 4.32 |  $-$ 12.79
10.10 |  $-$ 1.21 | 4.29 | 11.45 | 1.43 | 5.49 | 12.80 | 1.30 | 1.86 | 14.15 |  $-$ 1.24 |  $-$ 2.03 | 15.50 |  $-$ 5.59 |  $-$ 9.45 | 16.85 |  $-$ 4.26 |  $-$ 12.97
10.13 |  $-$ 1.15 | 4.25 | 11.48 | 1.41 | 5.50 | 12.83 | 1.25 | 1.81 | 14.18 |  $-$ 1.39 |  $-$ 2.22 | 15.53 |  $-$ 5.60 |  $-$ 9.55 | 16.88 |  $-$ 4.18 |  $-$ 13.11
10.16 |  $-$ 1.09 | 4.18 | 11.51 | 1.39 | 5.61 | 12.86 | 1.20 | 1.75 | 14.21 |  $-$ 1.51 |  $-$ 2.39 | 15.56 |  $-$ 5.61 |  $-$ 9.70 | 16.91 |  $-$ 4.12 |  $-$ 13.24
10.19 |  $-$ 1.04 | 4.15 | 11.54 | 1.38 | 5.65 | 12.89 | 1.14 | 1.69 | 14.24 |  $-$ 1.63 |  $-$ 2.54 | 15.59 |  $-$ 5.62 |  $-$ 9.79 | 16.94 |  $-$ 4.05 |  $-$ 13.40
10.22 |  $-$ 0.99 | 4.16 | 11.57 | 1.37 | 5.57 | 12.92 | 1.05 | 1.54 | 14.27 |  $-$ 1.79 |  $-$ 2.67 | 15.62 |  $-$ 5.63 |  $-$ 9.73 | 16.97 |  $-$ 3.95 |  $-$ 13.60
10.25 |  $-$ 0.96 | 4.23 | 11.60 | 1.36 | 5.49 | 12.95 | 0.98 | 1.41 | 14.30 |  $-$ 1.91 |  $-$ 2.77 | 15.65 |  $-$ 5.63 |  $-$ 9.63 | 17.00 |  $-$ 3.88 |  $-$ 13.77
10.28 |  $-$ 0.93 | 4.27 | 11.63 | 1.34 | 5.40 | 12.98 | 0.91 | 1.30 | 14.33 |  $-$ 2.06 |  $-$ 2.95 | 15.68 |  $-$ 5.62 |  $-$ 9.61 | 17.03 |  $-$ 3.80 |  $-$ 13.89
10.31 |  $-$ 0.89 | 4.17 | 11.66 | 1.34 | 5.33 | 13.01 | 0.83 | 1.17 | 14.36 |  $-$ 2.22 |  $-$ 3.18 | 15.71 |  $-$ 5.61 |  $-$ 9.63 | 17.06 |  $-$ 3.71 |  $-$ 13.97
10.34 |  $-$ 0.87 | 4.08 | 11.69 | 1.34 | 5.25 | 13.04 | 0.72 | 1.03 | 14.39 |  $-$ 2.35 |  $-$ 3.39 | 15.74 |  $-$ 5.61 |  $-$ 9.72 | 17.09 |  $-$ 3.62 |  $-$ 13.99
10.37 |  $-$ 0.86 | 3.92 | 11.72 | 1.34 | 5.21 | 13.07 | 0.64 | 0.99 | 14.42 |  $-$ 2.52 |  $-$ 3.71 | 15.77 |  $-$ 5.60 |  $-$ 9.86 | 17.12 |  $-$ 3.53 |  $-$ 14.03
10.40 |  $-$ 0.84 | 3.69 | 11.75 | 1.34 | 5.12 | 13.10 | 0.55 | 1.01 | 14.45 |  $-$ 2.65 |  $-$ 3.94 | 15.80 |  $-$ 5.59 |  $-$ 10.02 | 17.15 |  $-$ 3.39 |  $-$ 14.18
10.43 |  $-$ 0.84 | 3.50 | 11.78 | 1.34 | 4.97 | 13.13 | 0.45 | 0.99 | 14.48 |  $-$ 2.78 |  $-$ 4.14 | 15.83 |  $-$ 5.57 |  $-$ 10.13 | 17.18 |  $-$ 3.27 |  $-$ 14.39
10.46 |  $-$ 0.83 | 3.39 | 11.81 | 1.35 | 4.83 | 13.16 | 0.36 | 0.94 | 14.51 |  $-$ 2.95 |  $-$ 4.40 | 15.86 |  $-$ 5.56 |  $-$ 10.19 | 17.21 |  $-$ 3.15 |  $-$ 14.53
10.49 |  $-$ 0.83 | 3.41 | 11.84 | 1.37 | 4.70 | 13.19 | 0.25 | 0.83 | 14.54 |  $-$ 3.08 |  $-$ 4.58 | 15.89 |  $-$ 5.55 |  $-$ 10.19 | 17.24 |  $-$ 2.99 |  $-$ 14.73
10.52 |  $-$ 0.81 | 3.48 | 11.87 | 1.39 | 4.62 | 13.22 | 0.18 | 0.73 | 14.57 |  $-$ 3.23 |  $-$ 4.74 | 15.92 |  $-$ 5.55 |  $-$ 10.15 | 17.27 |  $-$ 2.86 |  $-$ 15.01
10.55 |  $-$ 0.78 | 3.63 | 11.90 | 1.42 | 4.54 | 13.25 | 0.12 | 0.63 | 14.60 |  $-$ 3.39 |  $-$ 4.91 | 15.95 |  $-$ 5.55 |  $-$ 10.15 | 17.30 |  $-$ 2.72 |  $-$ 15.37
10.58 |  $-$ 0.76 | 3.73 | 11.93 | 1.47 | 4.39 | 13.28 | 0.06 | 0.52 | 14.63 |  $-$ 3.52 |  $-$ 5.09 | 15.98 |  $-$ 5.55 |  $-$ 10.21 | 17.33 |  $-$ 2.57 |  $-$ 15.84
10.61 |  $-$ 0.71 | 3.80 | 11.96 | 1.50 | 4.23 | 13.31 | 0.02 | 0.44 | 14.66 |  $-$ 3.69 |  $-$ 5.36 | 16.01 |  $-$ 5.55 |  $-$ 10.25 | 17.36 |  $-$ 2.45 |  $-$ 16.56
10.64 |  $-$ 0.64 | 3.85 | 11.99 | 1.53 | 4.08 | 13.34 |  $-$ 0.01 | 0.34 | 14.69 |  $-$ 3.82 |  $-$ 5.57 | 16.04 |  $-$ 5.54 |  $-$ 10.31 | 17.39 |  $-$ 2.28 |  $-$ 17.13
10.67 |  $-$ 0.56 | 3.94 | 12.02 | 1.56 | 3.95 | 13.37 |  $-$ 0.03 | 0.22 | 14.72 |  $-$ 3.94 |  $-$ 5.78 | 16.07 |  $-$ 5.54 |  $-$ 10.44 | 17.42 |  $-$ 2.15 |  $-$ 16.85
10.70 |  $-$ 0.43 | 4.17 | 12.05 | 1.59 | 3.87 | 13.40 |  $-$ 0.03 | 0.11 | 14.75 |  $-$ 4.10 |  $-$ 6.04 | 16.10 |  $-$ 5.52 |  $-$ 10.55 | 17.45 |  $-$ 2.03 |  $-$ 16.86
10.73 |  $-$ 0.33 | 4.41 | 12.08 | 1.61 | 3.86 | 13.43 |  $-$ 0.04 | 0.03 | 14.78 |  $-$ 4.21 |  $-$ 6.21 | 16.13 |  $-$ 5.50 |  $-$ 10.59 | 17.48 |  $-$ 1.86 |  $-$ 16.85
10.76 |  $-$ 0.21 | 4.74 | 12.11 | 1.62 | 3.90 | 13.46 |  $-$ 0.04 | 0.01 | 14.81 |  $-$ 4.33 |  $-$ 6.45 | 16.16 |  $-$ 5.48 |  $-$ 10.64 | 17.51 |  $-$ 1.73 |  $-$ 16.83
10.79 |  $-$ 0.06 | 5.32 | 12.14 | 1.64 | 3.96 | 13.49 |  $-$ 0.04 | 0.01 | 14.84 |  $-$ 4.45 |  $-$ 6.71 | 16.19 |  $-$ 5.47 |  $-$ 10.68 |  |  |
10.82 | 0.05 | 5.71 | 12.17 | 1.66 | 3.97 | 13.52 |  $-$ 0.04 | 0.01 | 14.87 |  $-$ 4.54 |  $-$ 6.89 | 16.22 |  $-$ 5.45 |  $-$ 10.68 |  |  |

##  4 Summary

Following Paper I, by combining $\sim$ 0.7 million high-quality common stars
in LAMOST DR7 with the SCR method, we obtain $G-G_{\rm RP}$ and $G_{\rm
BP}-G_{\rm RP}$ color corrections as a function of $G$ magnitude for
$9.5<G<17.5$ to sub-mmag precision for Gaia EDR3. Our results confirm the
improvements in the calibration process of the EDR3. The color term of the
$G_{\rm BP}-G_{\rm RP}$ for bright stars is removed. The discontinuity caused
by the changes of instrument configurations is significantly reduced. Yet
modest systematic trends with $G$ magnitude are detected.

By comparing with the work of Yang et al. (2021), we find that our color
corrections of individual colors are subject to systematic errors in reddening
correction with the SFD map. In the cases of color-color diagram, our
corrections still achieve an unprecedented sub-mmag precision. Our work could
be beneficial to studies where a high-precision color-color diagram is
required, including estimates of Gaia photometric metallicities and
discrimination between binaries and single stars.

We acknowledge the anonymous referee for his/her valuable comments that
improve the quality of this paper significantly. This work is supported by
National Key Research and Development Program of China (NKRDPC) under grant
numbers 2019YFA0405503, 2019YFA0405504, and 2016YFA0400804, National Science
Foundation of China (NSFC) under grant numbers 11603002, 11988101, and
113300034, and Beijing Normal University grant No. 310232102. This work has
made use of data products from the Guoshoujing Telescope (the Large Sky Area
Multi-Object Fiber Spectroscopic Telescope, LAMOST). LAMOST is a National
Major Scientific Project built by the Chinese Academy of Sciences. Funding for
the project has been provided by the National Development and Reform
Commission. LAMOST is operated and managed by the National Astronomical
Observatories, Chinese Academy of Sciences. This work has made use of data
from the European Space Agency (ESA) mission Gaia
(https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and
Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

## References

  * Arenou et al. (2018) Arenou, F., Luri, X., Babusiaux, C., et al. 2018, A&A, 616, A17, doi: 10.1051/0004-6361/201833234
  * Bohlin (2014) Bohlin, R. C. 2014, AJ, 147, 127, doi: 10.1088/0004-6256/147/6/127
  * Casagrande & VandenBerg (2018) Casagrande, L., & VandenBerg, D. A. 2018, MNRAS, 479, L102, doi: 10.1093/mnrasl/sly104
  * Clem & Landolt (2013) Clem, J. L., & Landolt, A. U. 2013, AJ, 146, 88, doi: 10.1088/0004-6256/146/4/88
  * Evans et al. (2018) Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4, doi: 10.1051/0004-6361/201832756
  * Fabricius et al. (2020) Fabricius, C., Luri, X., Arenou, F., et al. 2020, arXiv e-prints, arXiv:2012.06242.  https://arxiv.org/abs/2012.06242
  * Gaia Collaboration et al. (2020) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2020, arXiv e-prints, arXiv:2012.01533.  https://arxiv.org/abs/2012.01533
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Guo et al. (2021) Guo, H. L., Chen, B. Q., Yuan, H. B., et al. 2021, ApJ, 906, 47, doi: 10.3847/1538-4357/abc68a
  * Lenz et al. (2017) Lenz, D., Hensley, B. S., & Doré, O. 2017, ApJ, 846, 38, doi: 10.3847/1538-4357/aa84af
  * Luo et al. (2015) Luo, A. L., Zhao, Y.-H., Zhao, G., et al. 2015, Research in Astronomy and Astrophysics, 15, 1095, doi: 10.1088/1674-4527/15/8/002
  * Maíz Apellániz & Weiler (2018) Maíz Apellániz, J., & Weiler, M. 2018, A&A, 619, A180, doi: 10.1051/0004-6361/201834051
  * Niu et al. (2021) Niu, Z., Yuan, H., & Liu, J. 2021, arXiv e-prints, ApJ, accepted, arXiv:2101.04290.  https://arxiv.org/abs/2101.04290
  * Peek & Graves (2010) Peek, J. E. G., & Graves, G. J. 2010, ApJ, 719, 415, doi: 10.1088/0004-637X/719/1/415
  * Planck Collaboration et al. (2014) Planck Collaboration, Abergel, A., Ade, P. A. R., et al. 2014, A&A, 571, A11, doi: 10.1051/0004-6361/201323195
  * Riello et al. (2020) Riello, M., De Angeli, F., Evans, D. W., et al. 2020, arXiv e-prints, arXiv:2012.01916.  https://arxiv.org/abs/2012.01916
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Weiler (2018) Weiler, M. 2018, A&A, 617, A138, doi: 10.1051/0004-6361/201833462
  * Yang et al. (2021) Yang, L., Yuan, H., Zhang, R., et al. 2021, arXiv e-prints, ApJ, accepted, arXiv:2101.00750.  https://arxiv.org/abs/2101.00750
  * Yuan et al. (2015) Yuan, H., Liu, X., Xiang, M., et al. 2015, ApJ, 799, 133, doi: 10.1088/0004-637X/799/2/133
  * Zhang & Yuan (2020) Zhang, R., & Yuan, H. 2020, ApJ, 905, L20, doi: 10.3847/2041-8213/abccc4
  * Zhao et al. (2012) Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Research in Astronomy and Astrophysics, 12, 723, doi: 10.1088/1674-4527/12/7/002


$\tilde{\boldsymbol{\mu}}$ will evolve as a moving average of incoming sensory
data $\tilde{\mathbf{y}}$ subject to a dynamical bias or "drag", which is a
consequence of the latent belief that hidden states $\mathbf{x}$ continuously
relax towards a fixed point at $\boldsymbol{\eta}$. Specifically, the beliefs
are constantly pulled closer to the data in order to minimize sensory
prediction errors $\tilde{\boldsymbol{\xi}}_{z}$; however, this process itself
incurs state prediction errors $\tilde{\boldsymbol{\xi}}_{\omega}$ that will
pull the beliefs back towards the fixed point. This constant tug of war
between sensory and process prediction errors can be shifted
disproportionately in one direction by adjusting the relative precisions of
the likelihood vs. dynamical models, respectively. If the process precision
$\tilde{\Pi}^{\boldsymbol{\omega}}$ is high relative to the observation
precision $\tilde{\Pi}^{\mathbf{z}}$, then the beliefs will tend to their
expected fixed point of $\boldsymbol{\eta}$. A similar enhancement of prior
bias can be achieved by increasing the drift rate $\alpha$ of the dynamics
model, which increases the force driving $\boldsymbol{\mu}$ towards
$\boldsymbol{\eta}$ — this was the approach taken in [96], for example.

Note that when numerically integrating the differential equation in (A.37)
with a forwards Euler scheme, one uses a finite number of iterations to update
the variational means $\tilde{\boldsymbol{\mu}}$, which we term
$n_{\text{InferIter}}$, and a step-size $\kappa_{\mu}$ which scales the size
of the increment to $\tilde{\boldsymbol{\mu}}$ [88]. In all simulations shown
here, we set $n_{\text{InferIter}}=1,\kappa_{\mu}=0.1$ (see Table E.1 for
details).

####  A.4.2 Closing the loop with observations and action

In order to interpret the random variables of the generative model as
representing behaviorally-relevant features of an agent’s world, we now turn
to specfiying the generative process, i.e., the actual physics of the world
that our self-propelled particle agents will inhabit. In this section we
detail how the observations $\tilde{\mathbf{y}}$ for a single agent are
generated from the positions and velocities of other active inference agents,
and how actions can be generated through active inference, which in this
contexts means changing continuous control variables using a gradient descent
on the same free energy used to derive the belief update equations of the
previous section.

We now shift our perspective to that of a single agent, hereafter referred to
as the focal individual or focal agent, and specify how its sensory data
$\tilde{\mathbf{y}}$ are generated. We start by describing univariate hidden
states and corresponding observations, where the true hidden variable is an
average nearest-neighbor distance $x_{h}$. We add the $h$ subscript to
distinguish these ‘real’ variables (hidden states, observations, noise terms)
from their representations in the generative model (e.g., $\tilde{x}$,
$\tilde{y}$).

We indicate the focal individual with index $i$; so the agent $i$-relative
hidden state $x_{h,i}$ denotes the average nearest-neighbor distance from the
perspective of agent $i$. This average distance $x_{h,i}$ is calculated from
the $K$ neighbors that form the interaction set $N_{in}$ of the $i$th focal
individual. How to define the interaction set $N_{in}$ is a choice to make in
each simulation, but for the case of recapitulating classical, distance-
dependent social forces models, we define $N_{in}$ as those neighbors that are
within a fixed distance $R_{0}$ of the focal individual’s position:

$\displaystyle x_{h,i}$ $\displaystyle\triangleq\frac{1}{K}\sum_{j\in
N_{in}}\|\Delta\mathbf{r}_{ij}\|$ $\displaystyle\ \textrm{where}\hskip
5.69054ptN_{in}$ $\displaystyle\triangleq\left\\{j\neq
i\,:\,\|\Delta\mathbf{r}_{ij}\|\leq R_{0}\right\\}$ (A.38) $\displaystyle\ K$
$\displaystyle\triangleq|N_{in}|$ $\displaystyle\ \Delta\mathbf{r}_{ij}$
$\displaystyle\triangleq\mathbf{r}_{j}-\mathbf{r}_{i}$ (A.39)

An additional filter on $N_{in}$ that is common to self-propelled particle
models, is to only include neighbors that subtend some angular extent (also
known as a ‘vision cone’ or ‘visual field’) relative to the focal agent’s
velocity vector $\mathbf{v}_{i}$. This is the approach taken in [97], for
instance, and in the simulations examined in the main text we do the same.

The vector $\mathbf{r}_{i}$ denotes the 2-D coordinate of the focal agent, and
$\mathbf{r}_{j}$ is that of neighbor $j$. $\mathbf{r}_{ij}$ thus represents
the relative displacement vector of neighbour $j$, from the perspective of the
focal agent $i$.

We also define the first temporal derivative of the local average distance
$x^{\prime}_{h,i}$:

$\displaystyle\tilde{x}_{h,i}$
$\displaystyle\triangleq(x_{h,i},x^{\prime}_{h,i})$ $\displaystyle
x^{\prime}_{h,i}$
$\displaystyle\triangleq\frac{dx_{h,i}}{dt}=\nabla_{\mathbf{r}_{i}}x_{h,i}\cdot\mathbf{v}_{i}+\sum_{j\in
N_{in}}\left(\nabla_{\mathbf{r}_{j}}x_{h,i}\cdot\mathbf{v}_{j}\right)$ (A.40)

where $\mathbf{v}_{j}$ is the velocity or heading vector of neighbour $j$. The
expression in (A.40) means that we can compute the first derivative or
velocity of the distance $x^{\prime}_{h,i}$ as a function of the positions and
velocities of all agents, as opposed to some discrete-time approximation,
e.g., $x^{\prime}_{h,i}\approx\frac{x_{h,i}(t+\Delta t)-x_{h,i}(t)}{\Delta t}$
for some small $\Delta t$. Note that this expression for $x^{\prime}_{h,i}$
assumes a local linearization of $x_{h,i}$ at the radius defined by $R_{0}$,
i.e., this linearization will be a poor predictor of the actual change in the
state $x_{h,i}(t+\Delta t)-x_{h,i}(t)$ when neighbors are instantaneously
leaving or entering the interaction set $N_{in}$. Observations
$\tilde{y}_{h,i}$ are perturbed versions of the hidden states with additive
generalised fluctuations $\tilde{z}_{h,i}$:

$\displaystyle y_{h,i}$ $\displaystyle=x_{h,i}+z_{h,i}$ $\displaystyle
y^{\prime}_{h,i}$ $\displaystyle=x^{\prime}_{h,i}+z^{\prime}_{h,i}$
$\displaystyle\textrm{where}\hskip 14.22636ptp(\tilde{z}_{h,i})$
$\displaystyle=N(\tilde{z}_{h,i};\mathbf{0},\tilde{\Sigma}_{z,h})$ (A.41)

In all simulations we parameterize the $\tilde{z}_{h,i}$ as independent
Gaussian variables, i.e.,

$\displaystyle\tilde{\Sigma}_{z,h}$
$\displaystyle=\begin{bmatrix}\sigma^{2}_{z,h}&0\\\
0&\sigma^{2}_{z^{\prime},h}\end{bmatrix}$ (A.42)

where the two variances $\sigma^{2}_{z,h}$ and $\sigma^{2}_{z^{\prime},h}$ can
be set independently. The ‘perception’ step of our active inference process
proceeds by providing these observations to the filtering equations in (A.37).
The result is that posterior means $\tilde{\mu}$ appears to track
$\tilde{x}_{h,i}$ over time, while additionally estimating its higher-order
motion (acceleration) via $\mu^{\prime\prime\prime}$.

Finally, we now furnish a scheme for updating actions by mapping the control
variables $a$ and sensorimotor contingency terms of (A.22) to the case of our
distance-tracking self-propelled agent.

We let actions be identifiable with the heading vector $\mathbf{v}_{i}$ of the
focal individual, i.e., $a=\mathbf{v}_{i}$. For the simulations presented in
the current paper, we always asserted that this heading have unit magnitude,
but in general this constraint is not necessary.

Given this definition of actions, we can unpack the sensorimotor contingency
term $\frac{\partial y(a)}{\partial a}$ that appeared in the active inference
control equation of (A.22), now letting $a=\mathbf{v}$ and turning partial
derivatives into Jacobians to account for vectorial nature of actions (being a
velocity in 2-D) and observations (being comprised of two generalized
coordinates):

$\displaystyle\frac{d\mathbf{v}_{i}}{dt}$
$\displaystyle=-\nabla_{\mathbf{v}_{i}}\tilde{y}_{h,i}(\mathbf{v}_{i})^{\top}\nabla_{\tilde{y}_{h,i}(\mathbf{v}_{i})}F_{L}$
(A.43)

Note here that observations $\tilde{y}_{h,i}$ are a function of actions; this
is because observations are a linear function of hidden states, which
themselves are linear in the velocity vector of the focal individual
$\mathbf{v}_{i}$ via the relation in (A.40). Importantly, however, the
distance observation $y_{h,i}$ does not directly depend on the
$\mathbf{v}_{i}$ — only the distance velocity $y^{\prime}_{h,i}$ does. This
means the sensorimotor contingency in (A.43) is comprised of non-zero partial
derivatives only for $y^{\prime}_{h,i}$:

$\displaystyle\nabla_{\mathbf{v}_{i}}\tilde{y}_{h,i}(\mathbf{v}_{i})$
$\displaystyle=\begin{bmatrix}\nabla_{\mathbf{v}_{i}}y_{h,i}(\mathbf{v}_{i})\\\
\nabla_{\mathbf{v}_{i}}y^{\prime}_{h,i}(\mathbf{v}_{i})\end{bmatrix}=\begin{bmatrix}\mathbf{0}\\\
\nabla_{\mathbf{r}_{i}}x_{h,i}\end{bmatrix}$ (A.44)

This has an important consequence for action, when we consider the form of the
second part of the action update in (A.43), the free energy gradient term
$\nabla_{\tilde{y}_{h,i}}F_{L}$:

$\displaystyle\nabla_{\tilde{y}_{h,i}}F_{L}$
$\displaystyle=\tilde{\xi}_{z}=\tilde{\Pi}_{z}\tilde{\varepsilon}_{z}=\begin{bmatrix}\Gamma_{z}(y_{h,i}-\mu)\\\
2\Gamma_{z}\lambda_{z}^{2}(y^{\prime}_{h,i}-\mu^{\prime})\end{bmatrix}$ (A.45)

The free energy gradient with respect to observations is simply the
generalized (precision-weighted) sensory error $\tilde{\xi}_{z}$, which we
have written in terms of the observations $\tilde{y}_{h,i}$, posterior beliefs
$\tilde{\mu}$ and precision parameters $\Gamma_{z},\lambda_{z}$. The sparse
form of the sensorimotor contingency in (A.44) means that the 0th-order
prediction error $\xi_{z}$ will have no effect on behavior and only the
velocity prediction errors $\xi^{\prime}_{z}$ will be relevant for the update
to $\mathbf{v}_{i}$, i.e.,

$\displaystyle\frac{d\mathbf{v}_{i}}{dt}$
$\displaystyle=-\left(\xi_{z}\underbrace{\nabla_{\mathbf{v}_{i}}y_{h,i}(\mathbf{v}_{i})}_{=\mathbf{0}}+\xi^{\prime}_{z}\nabla_{\mathbf{v}_{i}}y^{\prime}_{h,i}(\mathbf{v}_{i})\right)$
$\displaystyle=-\xi^{\prime}_{z}\nabla_{\mathbf{r}_{i}}x_{h,i}$
$\displaystyle=2\Gamma_{z}\lambda_{z}^{2}(y^{\prime}_{h,i}-\mu^{\prime})\Delta\hat{\mathbf{r}}$
$\displaystyle\textrm{where}\hskip 5.69054pt\Delta\hat{\mathbf{r}}$
$\displaystyle=\frac{1}{K}\sum_{j\in
N_{in}}\frac{\Delta\mathbf{r}_{ij}}{\|\Delta\mathbf{r}_{ij}\|}$ (A.46)

Note that, as for the inference update in (A.37), we update $\mathbf{v}_{i}$
using a fixed number of action iterations $n_{\text{ActionIter}}$ and step-
size $\kappa_{a}$, where here we set $n_{\text{ActionIter}}=1,\kappa_{a}=0.1$.
This action update equation has a few key implications for the behavior of
active inference agents equipped with this type of generative model, and its
relationship to ‘classical’ self-propelled particle models like the Couzin-
Aoki model and the Reynolds or BOIDS model [98, 97, 99]. The first is the fact
that the sensorimotor contingency is identical to the ‘social force’ vector
used to drive interactions in self-propelled particle models
$\Delta\hat{\mathbf{r}}$; this the average of the vectors pointing from each
neighbor in the interacting set to the focal agent’s position
$\mathbf{r}_{i}$. The sign of the precision-weighted prediction error
$\epsilon^{\prime}_{z}$ determines whether the social force is attractive
(pointing towards other agents) or repulsive (pointing away from other
agents). Secondly, the fact that actions only depend on velocity observations,
rather than state observations, means that agents will adjust their heading
according to how the (sensed) distance is instantaneously changing (its
velocity), rather than its value. This lends action a predictive, anticipatory
power and accounts for why we observe robust polarized motion in the absence
of an explicit alignment term like in classic self-propelled particle models
[100, 97]. The alignment-like forces emerges from the fact that the velocity
vectors of other agents $v_{j},j\in N_{in}$ are integrated into the
computation of $y^{\prime}_{h,i}$ via the relation in the second line of
(A.40).

One of the defining features of other self-propelled particle models like the
Couzin-Aoki model [98, 97] is the presence and priorization of interaction
zones. The two main zones used in these models, and which on their own are
sufficient for group cohesion, are a narrow repulsion zone defined by some
radius $r_{r}$ and a wider attraction zone with radius $r_{a}$, where
$r_{a}>r_{r}$. Neighboring agents within the repulsive radius exert repulsive
forces on the focal agent, while those beyond the repulsion radius but within
the attraction zone exert attractive forces, where the difference between
attraction and repulsion is given by the sign of the force vector
$\Delta\hat{\mathbf{r}}$. The active inference model leads to an effective
notion of zones, but rather than being explicitly encoded, these zones emerge
through the fixed-point attractor $\eta$ parameterizing the generative model’s
dynamics model $\mathbf{f}$. This is made clear when we examine the precision-
weighted prediction error $\xi^{\prime}_{z}$, which itself is a function of
velocity observations $y^{\prime}_{h,i}$ and velocity beliefs $\mu^{\prime}$.
Consider the limiting case of when inference is strongly biased by the
dynamics model $\mathbf{f}$ (i.e., in the case that
$\Gamma_{\omega}>\Gamma_{z}$ or large $\alpha$); the generalised beliefs
$\tilde{\mu}$ will be strongly drawn to the setpoint $\eta$ of the dynamics
prior, i.e.,

$\displaystyle\tilde{\mu}=\begin{bmatrix}\mu\\\ \mu^{\prime}\\\
\mu^{\prime\prime}\end{bmatrix}$ $\displaystyle\approx\begin{bmatrix}\eta\\\
0\\\ 0\end{bmatrix}$

Under this assumption, the precision-weighted prediction error
$\xi^{\prime}_{z}$ approximates $2\Gamma_{z}\lambda_{z}^{2}y^{\prime}_{h,i}$,
and thus signals whether neighbors are instantaneously approaching or moving
away from the focal agent, where $\xi^{\prime}_{z}<0$ indicates they are
approaching and $\xi^{\prime}_{z}>0$ indicates they are moving away. This in
turn determines whether the update to the focal agent’s action
$\mathbf{v}_{i}$ is repulsive or attractive, as its sign determines the
direction of the social force vector $\Delta\hat{\mathbf{r}}$. Although the
first order distance $y_{h,i}$ does not directly drive action, it does so
indirectly through its effect on inference of $\mu^{\prime}$. If we consider
the case when the sensed distance $y_{h,i}$ drops below the setpoint $\eta$,
then one can reason through the cascade of prediction errors that ultimately
lead to a repulsive force. As a direct consequence of a drop in $y_{h,i}$
below $\mu$, sensory prediction errors $\xi_{z}$ will become negative, whose
minimization will require $\mu$ to move below $\eta$. This process in turn
incurs slower-moving (negative) model prediction errors $\xi_{\omega}$, whose
minimization drives $\mu$ back to its fixed point of $\eta$, given the dynamic
constraint for the beliefs to relax to their fixed point. In order to
accomplish this upward movement of $\mu$, either the rate of change of $\mu$
or the sensed distance itself must be positive, i.e., $\dot{\mu}>0$ or
$y_{h,i}>0$. In the absence of positive $y_{h,i}$, model prediction errors
will drive $\mu^{\prime}$ (and hence $\dot{\mu}$) above $0$. This temporarily
sets a larger radius of repulsion, i.e., a larger range of $y^{\prime}_{h,i}$
for which $\xi^{\prime}_{z}$ is negative and for which repulsive forces impact
the focal agent’s velocity. This causes the agent to move away from its
neighbors and thus further increase $y_{h,i}$, under the assumption that the
agent’s prediction of the distance dynamics are correlated with the true
change in $x_{h,i}$. Belief updating and action thus work together to
accelerate the return of $\mu$ towards $\eta$ and
$\tilde{\xi}_{z},\tilde{\xi}_{\omega}$ towards $0$; for this reason active
inference is often described as an account of action and perception driven by
‘self-fulfilling prophecy’ [88].

In order to imbue action with a more direct coupling to the neighbors‘
distances as is done in the classical self-propelled particle models, rather
than the velocity of the distance, one could hand-craft the sensorimotor
contingency term $\nabla_{\mathbf{v}_{i}}\tilde{y}_{h,i}$ to enforce a
coupling between $y_{h,i}$ and $\mathbf{v}_{i}$. This would render the action
rule equivalent to a ‘soft’-form of PD control [96], where errors on both the
first order state $(y_{h,i}-\mu)\approx(y_{h,i}-\eta)$ and its derivative
$(y^{\prime}_{h,i}-\mu^{\prime})\approx y^{\prime}_{h,i}$ would drive changes
to the velocity.

###  A.5 Extending to multiple sensory sectors

The results of the previous sections can be straightforwardly extended to the
multivariate case as explored in the main text. The focal agent now senses the
local distance computed across a set of distinct sensory sectors. For the
model explored in the current work, we split up the computation of the local
distance variable into a set of $L$ sensory adjacent sectors that comprise an
arc of a given angle, relative to the agent’s heading vector $\mathbf{v}_{i}$.
We define the multivariate distance hidden state as follows (dropping the
focal agent index $i$ from the sector-specific hidden states to avoid
subscript overload):

$\displaystyle\mathbf{x}_{h,i}$ $\displaystyle=\begin{bmatrix}x_{h,1}\\\
x_{h,2}\\\ \vdots\\\ x_{h,L}\end{bmatrix}$ (A.47)
$\displaystyle\textrm{where}\hskip 5.69054ptx_{h,l}$
$\displaystyle\triangleq\frac{1}{K_{l}}\sum_{j\in
N_{l}}\|\Delta\mathbf{r}_{ij}\|$

where $N_{l}$ is the set of neighbors in the $l$th sensory sector, and
$K_{l}=|N_{l}|$ (c.f., (A.39)). As with the scalar hidden state defined above,
we also equip the vector of sector distances $\mathbf{x}_{h,i}$ with
corresponding sector-specific, generalized observations
$\tilde{\mathbf{y}}_{h,i}$, i.e.

$\displaystyle\mathbf{y}_{h,i}$
$\displaystyle=\mathbf{x}_{h,i}+\mathbf{z}_{h,i}$ (A.48)
$\displaystyle\mathbf{y}^{\prime}_{h,i}$
$\displaystyle=\mathbf{x}^{\prime}_{h,i}+\mathbf{z}^{\prime}_{h,i}$
$\displaystyle\textrm{where}\hskip 14.22636ptp(\tilde{\mathbf{z}}_{h,i})$
$\displaystyle=\mathcal{N}(\tilde{\mathbf{z}}_{h,i};\mathbf{0},\tilde{\Sigma}_{\mathbf{z},h})$
(A.49)

such that the focal individual now observes a vector of local (noise-
perturbed) distances and their first orders of motion. Note that the
generalized covariance matrix here $\tilde{\Sigma}_{\mathbf{z},h}$ is now a
$2L\times 2L$ size matrix, that encodes the covariance structure between
sector-specific noise and their generalized orders. For all simulations we
generated uncorrelated noise across the different sectors, although spatially-
smooth noise could be modelled by introducing off diagonal elements in
$\tilde{\Sigma}_{\mathbf{z},h}$, i.e., $\mathbb{E}[z_{h,l}z_{h,k}]\neq 0$.

The agent’s generative model is also extended to the multivariate state-space
formulation we began with, using a vector of generalised hidden states
$\tilde{\mathbf{x}}=(\tilde{x}_{1},\tilde{x}_{2},...,\tilde{x}_{L})$ to
estimate the local distance within each sensory sector. Belief-updating
consists in updating a vector of generalised means $\tilde{\boldsymbol{\mu}}$
through integration of (A.37).

The action update has an identical form as before, except now the sensorimotor
contingency term
$\nabla_{\mathbf{v}_{i}}\tilde{\mathbf{y}}_{h,i}(\mathbf{v_{i}})$ is a
collection of partial derivative vectors, one for each sensory sector:

$\displaystyle\frac{d\mathbf{v}_{i}}{dt}$
$\displaystyle=-\nabla_{\mathbf{v}_{i}}\tilde{\mathbf{y}}_{h,i}(\mathbf{v}_{i})^{\top}\nabla_{\tilde{\mathbf{y}}_{h,i}}F_{L}$
$\displaystyle\nabla_{\mathbf{v}_{i}}\tilde{\mathbf{y}}_{h,i}(\mathbf{v_{i}})$
$\displaystyle=\begin{bmatrix}\nabla_{\mathbf{v}_{i}}\mathbf{y}_{h,i}(\mathbf{v}_{i})\\\
\nabla_{\mathbf{v}_{i}}\mathbf{y}^{\prime}_{h,i}(\mathbf{v}_{i})\end{bmatrix}=\begin{bmatrix}\mathbf{0}\\\
\vdots\\\ \mathbf{0}\\\ \nabla_{\mathbf{r}_{i}}x_{h,1}\\\
\nabla_{\mathbf{r}_{i}}x_{h,2}\\\ \vdots\\\
\nabla_{\mathbf{r}_{i}}x_{h,L}\end{bmatrix}$ (A.50)

The last $L$ rows of this Jacobian matrix encode the gradients of the sector-
specific distance velocities $y^{\prime}_{h,l}$ with respect to the focal
agent’s action; these partial derivatives are vectors pointing from the
average position of the neighbors in sector $l$ towards the focal individual.
When we combine the Jacobian matrix in (A.50) with the sensory prediction
error term $\tilde{\mathbf{y}}_{h,i}$ (i.e., the free energy gradients
$\nabla_{\tilde{\mathbf{y}}_{h,i}}F_{L}$), we are left with the following
update for the velocity:

$\displaystyle\frac{d\mathbf{v}_{i}}{dt}$
$\displaystyle=\boldsymbol{\xi}^{\prime}_{z}\cdot\Delta\hat{\mathbf{R}}=\begin{bmatrix}\xi^{\prime}_{z,1}&\xi^{\prime}_{z,2}&\ldots&\xi^{\prime}_{z,L}\end{bmatrix}\cdot\begin{bmatrix}\Delta\hat{\mathbf{r}}_{1}\\\
\Delta\hat{\mathbf{r}}_{2}\\\ \vdots\\\
\Delta\hat{\mathbf{r}}_{L}\end{bmatrix}$
$\displaystyle=\sum_{l=1}^{L}\xi^{\prime}_{z,l}\Delta\hat{\mathbf{r}}_{l}=2\Gamma_{z}\lambda_{z}^{2}\sum_{l=1}^{L}(y^{\prime}_{h,l}-\mu^{\prime}_{l})\Delta\hat{\mathbf{r}}_{l}$
(A.51) $\displaystyle\textrm{where}\hskip 5.69054pt\Delta\hat{\mathbf{r}}_{l}$
$\displaystyle=\frac{1}{K_{l}}\sum_{j\in
N_{l}}\frac{\Delta\mathbf{r}_{ij}}{||\Delta\mathbf{r}_{ij}||}$

The action thus becomes a weighted sum of ‘sector-vectors’
$\Delta\hat{\mathbf{r}}_{l}$, which are vectors pointing from the focal
agent’s position $\mathbf{r}_{i}$ towards the average position of the
neighbors in $N_{l}$. The weights that scale each $\Delta\hat{\mathbf{r}}_{l}$
are the precision-weighted prediction errors associated with velocity
observations emanating from the appropriate sector
$\xi^{\prime}_{z,l}\propto(y^{\prime}_{h,l}-\mu^{\prime}_{l})$. The fact we
can pull the spatiotemporal precision terms $2\Gamma_{z}\lambda_{z}^{2}$
outside the sum over sector-vectors, inherits from a between-sector
independence assumption, built into the agent’s sensory likelihood model
$P(\tilde{\mathbf{y}}|\tilde{\mathbf{x}})$ (see (A.31)). If the generative
model allowed for between-sector correlations (i.e. $\Sigma^{\mathbf{z}}$ was
not diagonal), then the action update would include cross-terms that couple
prediction errors from one sector to the sector-vector from another sector.

An active inference agent equipped with such a multivariate representation of
the local neighbor-distances thus engages in a sort of ‘predictive balancing-
act’, differentially responding more or less to each part of its sensory field
in accordance with how much sensations deviate from their posterior
expectations $\mu^{\prime}_{l}$, where the sign and degree of this deviation
is scored by $\xi^{\prime}_{z,l}$.

### Gaussianity of the generative model

One may question the use of a Gaussian form for the generalized noises
$\tilde{\mathbf{z}},\tilde{\boldsymbol{\omega}}$ and a linear form for
$\tilde{\mathbf{f}}$ and $\tilde{\mathbf{g}}$. These assumptions guarantee a
simple generative model whose free energy gradients (with respect to both the
state belief $\tilde{\boldsymbol{\mu}}$ and action $\mathbf{v}$) are linear in
the generalized distance observations $\tilde{\mathbf{y}}$. The reason we use
a Gaussian form here, is to achieve those simple, linear free energy
gradients, i.e., those which align with the classical social forces seen in
SPP models (more specifically, the selective attraction-and-repulsion forces
first described in [101]); in other words, it is exactly an active inference
model equipped with this sort of model (and observations with two generalized
coordinates, aka position and velocity) that acts in a way equivalent to being
driven by vectorial social forces. Note that this also means the Laplace
approximation to the posterior over $\tilde{\mathbf{x}}$ is not an
approximation at all, but leads to exact inference.

Existing social force models where the forces are linear in the state, can
thus be read as special cases of active inference agents equipped with
Gaussian generative models – this is a close cousin to the relationship
between linear PID controllers and Gaussian active inference models derived in
[96]. However, if we were to use a different, non-Gaussian of the generative
model (either through breaking Gaussianity of the noise model or by
introducing nonlinear forms of $\tilde{\mathbf{f}},\tilde{\mathbf{g}}$), then
the updates to the beliefs and actions would no longer be linear functions of
the generalized posterior means $\tilde{\boldsymbol{\mu}}$ and observations
$\tilde{\mathbf{y}}$, but could be arbitrary, nonlinear functions of these
variables – and the resulting social forces would have no guaranteed
interpretation in terms of attraction and repulsion. Such extensions are an
interesting avenue for future work; one interesting possibility would be to
explore new extensions of social forces, where the noises
$\tilde{\mathbf{z}},\tilde{\boldsymbol{\omega}}$ are assumed to belong to the
exponential family and that the dynamics model $\tilde{\mathbf{f}}$ is smooth
and contains attracting fixed points (i.e., regions of $\mathbf{\tilde{x}}$
where the derivatives of $\tilde{\mathbf{f}}$ vanish). A obvious consequence
of such assumptions, is that when the posterior belief
$\tilde{\boldsymbol{\mu}}$ is near these fixed points, i.e.,
$\tilde{\boldsymbol{\mu}}\approx\arg\min_{\tilde{\mathbf{x}}}\nabla_{\tilde{\mathbf{x}}}\tilde{\mathbf{f}}(\tilde{\mathbf{x}})\implies\tilde{\boldsymbol{\mu}}^{\prime}=0$,
we would still recover attractive and repulsive social forces whose magnitude
would be approximately linear in the distance observations
$\tilde{\mathbf{y}}$. This is equivalent to making a locally-quadratic (i.e.,
Laplace) approximation to the free energy landscape around its minima – this
corresponds to just those points in belief-space where the posterior is well-
approximated by a Gaussian.

##  Appendix B Alignment forces from active inference on angles

In previous sections we have shown how repulsive and attractive forces emerge
from active inference models in which the agent entertains a latent
representation of the average local distance between itself and its neighbors,
and how its heading direction couples to (the derivative of) that variable. In
this section we derive alignment-based social forces, like those that appear
in the Reynolds, Couzin, and Vicsek models [99, 97, 100], as a special case of
active inference, where an agent infers the (cosine) angle between its own
heading and that of its neighbors, and acts under the prior belief that this
angle tends to $0$.

As before, we start with a generative model that represents a generalised
latent variable $\tilde{x}_{\phi}$ that evolves in time with Gaussian additive
fluctuations $\tilde{\omega}_{\phi}$. We use the $\phi$ subscript to
distinguish this angle-tracking latent variable from the distance-tracking
variable of the previous section. We truncate the generalized representation
of this state at second order, i.e.
$\tilde{x}_{\phi}=\\{x_{\phi},x^{\prime}_{\phi}\\}$, leading to a dynamical
equation and corresponding likelihood of the following form:

$\displaystyle\dot{x}_{\phi}$
$\displaystyle=-\alpha_{\phi}(x_{\phi}-1)+\omega_{\phi}$
$\displaystyle\dot{x}^{\prime}_{\phi}$
$\displaystyle=-\alpha_{\phi}x_{\phi}^{\prime}+\omega^{\prime}_{\phi}$
$\displaystyle\implies p(D\tilde{x}_{\phi}|\tilde{x}_{\phi})$
$\displaystyle=\mathcal{N}(D\tilde{x}_{\phi};\tilde{f}_{\phi},\tilde{\Sigma}_{\omega_{\phi}})$
(B.52) $\displaystyle\textrm{where}\hskip 5.69054pt\tilde{f}_{\phi}$
$\displaystyle=\begin{bmatrix}-\alpha_{\phi}(x_{\phi}-1)\\\
-\alpha_{\phi}x^{\prime}_{\phi}\end{bmatrix},\hskip
5.69054pt\tilde{\Sigma}_{\omega_{\phi}}=\begin{bmatrix}\sigma^{2}_{\omega_{\phi}}&0\\\
0&\sigma^{2}_{\omega^{\prime}_{\phi}}\end{bmatrix}$ (B.53)

The observation model describes a mapping from the $0$th-order state to a
corresponding observation thereof, perturbed again by Gaussian innovations:

$\displaystyle y_{\phi}$ $\displaystyle=x_{\phi}+z_{\phi}$
$\displaystyle\implies p(y_{\phi}|x_{\phi})$
$\displaystyle=\mathcal{N}(y_{\phi};x_{\phi},\sigma^{2}_{z_{\phi}})$ (B.54)

Following the same steps as we did previously for the multivariate, distance-
tracking generative model, we can write down the Laplace-approximated
variational free energy of this model as a quadratic function of the
observations and generalized means $\tilde{\mu}_{\phi}$:

$\displaystyle F_{L}$
$\displaystyle\propto\varepsilon_{z_{\phi}}^{\top}\Pi_{z_{\phi}}\varepsilon_{z_{\phi}}+\tilde{\varepsilon}_{\omega_{\phi}}^{\top}\tilde{\Pi}_{\omega_{\phi}}\tilde{\varepsilon}_{\omega_{\phi}}$
$\displaystyle\textrm{where}\hskip 5.69054pt\varepsilon_{z_{\phi}}$
$\displaystyle\triangleq y_{\phi}-\mu_{\phi}$
$\displaystyle\tilde{\varepsilon}_{\omega_{\phi}}$ $\displaystyle\triangleq
D\tilde{\mu}_{\phi}-\tilde{f}_{\phi}$

The agent performs a gradient descent on $F_{L}$ to infer the value of
$\tilde{x}_{\phi}$ in light of sensory observations. This inference is encoded
by a Gaussian variational posterior with mean $\tilde{\mu}_{\phi}$. As before,
we can tune model parameters such that inference is strongly biased by the
dynamics model $\tilde{f}_{\phi}$, where the zeroth-order of motion
$\mu_{\phi}\approx 1$. The reason we set the set-point at $1$ becomes evident
when we consider the generation of sensory data and actions.

Assume that the focal agent with index $i$ observes the local average cosine
angle between its own heading vector $\mathbf{v}_{i}$ and those of its
neighbors $\mathbf{v}_{j},j\in N_{in}$, where neighbors are once again defined
by membership in some interaction zone333For notational convenience and
because it doesn’t change the derivations, we omit observation noise on
$y_{\phi}$.:

$\displaystyle y_{\phi}$ $\displaystyle=\frac{1}{K}\sum_{j\in
N_{in}}\mathbf{v}_{i}^{\top}\mathbf{v}_{j}=\langle\cos(\theta_{ij}\rangle_{N_{in}}$
(B.55)

where the equivalence between the dot products and the cosine angle is assured
when we assume all $\mathbf{v}_{k},k\in\\{i\\}\cup N_{in}$ have unit
magnitude. Recall that if two unit-magnitude vectors
$\mathbf{v}_{i},\mathbf{v}_{j}$ are parallel, their dot product (cosine angle)
is 1. When we once again assume that agents act by adjusting their heading
direction, then the action update given the continuous active inference rule
in (A.22) has the following form:

$\displaystyle\frac{d\mathbf{v}_{i}}{dt}$
$\displaystyle=-\frac{1}{\sigma^{2}_{z_{\phi}}}(y_{\phi}-\mu_{\phi})\hat{\mathbf{v}}\approx(1-y_{\phi})\hat{\mathbf{v}}$
(B.56) $\displaystyle\textrm{where}\hskip 5.69054pt\hat{\mathbf{v}}$
$\displaystyle=\frac{1}{K}\sum_{j\in N_{in}}\mathbf{v}_{j}$ (B.57)

The approximation in the first line holds when we assume the sensory variance
$\sigma^{2}_{z_{\phi}}$ is $1$ and the dynamics prior (either via increasing
$\alpha$ or decreasing $\sigma^{2}_{\omega_{\phi}}$) dominates inference such
that $\mu_{\phi}\approx 1$. In this case, the focal agent $i$ then updates its
velocity using the average neighbor velocity. This is proportional to the
alignment force in e.g. [100, 97], except that it is also scaled by how
unaligned the focal individual is with its neighbourhood, scored by
$1-y_{\phi}$.

##  Appendix C Online parameter estimation

In this section we derive update rules for the generative model parameters
using a simple gradient descent scheme on the Laplace-approximated variational
free energy. In the active inference literature this process of updating
parameters, as opposed to beliefs about states, is often analogized to online
learning or neural plasticity [102, 103].

###  C.1 Updating sensory smoothness

In this section we derive an update equation for the sensory smoothness
parameter $\lambda_{z}$, which captures the generative model’s assumptions
about the temporal autocorrelation structure of sensory noise $z$.

Recall the formulation of state space models in generalized coordinates of
motion in Section A.1. In addition to providing a concise description of local
paths of the state $\vec{x}_{t}$ in terms of its higher derivatives
$x^{\prime},x^{\prime\prime},...,x^{[n]}$, stochastic differential equations
in generalized coordinates also allow one to express serial correlations in
the noises at the first order $z$, by assuming that it can be differentiated
(has non-zero, smooth autocovariance) and represented in terms of hierarchical
or generalized noises
$z^{\prime},z^{\prime\prime},z^{\prime\prime\prime},...,z^{[n]}$.

Recall the parameterization of the generalized sensory precision
$\tilde{\Pi}^{z}$ as a factorization into two precision matrices, that
respectively represent agent’s beliefs about the ‘spatial’ and ‘temporal’
covariance structure. We parameterize these with the two precision parameters
$\Gamma_{z}$ and $\lambda_{z}$. $\Gamma_{z}$ encodes the agent’s belief about
the overall magnitude of the fluctuations, and $\lambda_{z}$ encodes beliefs
about their their serial correlations in time, assuming a Gaussian form for
their autocorrelation:

$\displaystyle\tilde{\Pi}^{z}$
$\displaystyle=S(\lambda_{z})\otimes\Pi(\Gamma_{z})$
$\displaystyle\Pi(\Gamma_{z})$ $\displaystyle=\begin{bmatrix}\Gamma_{11}&\\\
&\Gamma_{22}&&\\\ &&\ddots&\\\ &&&\Gamma_{LL}\end{bmatrix}$ $\displaystyle
S(\lambda_{z})$
$\displaystyle=\begin{bmatrix}1&0&-\frac{1}{2\lambda^{2}_{z}}&\ldots\\\
0&\frac{1}{2\lambda^{2}_{z}}&0&\\\
-\frac{1}{2\lambda^{2}_{z}}&0&\frac{3}{4\lambda^{4}_{z}}&\\\
\vdots&&&\ddots\end{bmatrix}^{-1}$ (C.58)

We implement a form of behavioral plasticity by allowing agents to update
$\lambda_{z}$ using observations. We accomplish this using a gradient descent
on variational free energy:

$\displaystyle\frac{d\lambda_{z}}{dt}$
$\displaystyle=-\kappa_{\theta}\frac{\partial F}{\partial\lambda_{z}}$ (C.59)

where the ‘learning rate’ $\kappa_{\theta}$ is typically set to be at least an
order of magnitude lower than the update rate of inference $\kappa_{\mu}$; in
all simulations we use $\kappa_{\theta}=0.001$ and $n_{\text{LearnIter}}=1$
iteration. This enforces a separation of timescales that is typical in
generalized filtering and state-space models that perform simultaneous state-
and parameter-estimation [87, 85, 96].

To compute the gradients of the variational free energy with respect to
$\lambda_{z}$, we can start by expressing those components of the (Laplace-
approximated) variational free energy that depend on $\lambda_{z}$:

$\displaystyle F(\lambda_{z})$
$\displaystyle=\tilde{\varepsilon}_{\mathbf{z}}^{\top}\tilde{\Pi}^{\mathbf{z}}\tilde{\varepsilon}_{\mathbf{z}}-\ln\left(\det{\tilde{\Pi}^{\mathbf{z}}}\right)$
(C.60)

where we only have included the terms that depend on the sensory precision
$\tilde{\Pi}^{z}$ due to its dependence on $\lambda_{z}$. The full gradient is
then simply:

$\displaystyle\frac{\partial
F}{\partial\lambda_{z}}=\frac{\tilde{\varepsilon}_{\mathbf{z}}^{\top}\tilde{\Pi}^{\mathbf{z}}\tilde{\varepsilon}_{\mathbf{z}}}{\partial\lambda_{z}}-\frac{\partial\ln\left(\det\tilde{\Pi}^{\mathbf{z}}\right)}{\partial\lambda_{z}}$
(C.61)

Starting with the case of a single sensory sector $L=1$, then the generalized
prediction error $\tilde{\varepsilon}_{\mathbf{z}}$ is a vector of prediction
errors, one for each order of motion:
$\tilde{\varepsilon}_{z}=\\{\varepsilon_{z},\varepsilon^{\prime}_{z},\varepsilon^{\prime\prime}_{z},...\\}$
where a sensory prediction error at a given order of motion is simply:
$\varepsilon^{[n]}_{z}=y^{[n]}-\tilde{g}^{[n]}$, where the $n$ subscript
refers to an order of differentiation. In the case of $3$ generalized
coordinates for the simple scalar case:

$\displaystyle\frac{\partial F}{\partial\lambda_{z}}$
$\displaystyle=4\Gamma_{z}\lambda_{z}(\varepsilon^{\prime}_{z})^{2}+\varepsilon^{\prime\prime}_{z}(8\Gamma_{z}\varepsilon^{\prime\prime}_{z}\lambda_{z}^{3}+2\Gamma_{z}\varepsilon_{z}\lambda_{z})+2\Gamma_{z}\lambda_{z}\varepsilon_{z}\varepsilon^{\prime\prime}_{z}-\frac{6}{\lambda_{z}}$
$\displaystyle=4\Gamma_{z}\lambda_{z}(\varepsilon_{z}\varepsilon^{\prime\prime}_{z}+(\varepsilon^{\prime}_{z})^{2}+2\lambda_{z}^{2}(\varepsilon^{\prime\prime}_{z})^{2})-\frac{6}{\lambda_{z}}$
(C.62)

Meaning that the update for the $\lambda_{z}$ parameter can be simplified to
(omitting the learning rate $\kappa_{\theta}$):

$\displaystyle\frac{d\lambda_{z}}{dt}$
$\displaystyle=-4\Gamma_{z}\lambda_{z}(\varepsilon_{z}\varepsilon^{\prime\prime}_{z}+(\varepsilon^{\prime}_{z})^{2}+2\lambda_{z}^{2}(\varepsilon^{\prime\prime}_{z})^{2})+\frac{6}{\lambda_{z}}$
(C.63)

In the case of the distance-tracking generative model we explore in the main
text, we assume that the agents can only observe the $0$th (position, $y$) and
1st (velocity, $y^{\prime}$) orders of motion of the hidden states
$\tilde{x}$. This means there are no longer 2nd-order prediction errors
$\varepsilon^{\prime\prime}_{z}$ and the update becomes even simpler:

$\displaystyle\frac{d\lambda_{z}}{dt}$
$\displaystyle=-4\Gamma_{z}\lambda_{z}(\varepsilon^{\prime}_{z})^{2}+\frac{6}{\lambda_{z}}$
$\displaystyle\approx-4\Gamma_{z}\lambda_{z}(y^{\prime}_{h,i})^{2}+\frac{6}{\lambda_{z}}$
(C.64)

where approximation in the second line results in the case of ‘biased’
inference, i.e., $\mu\approx\eta\implies\mu^{\prime}\approx 0$, allowing us to
replace the velocity prediction error $y^{\prime}_{h,i}-\mu^{\prime}$ with
$y^{\prime}_{h,i}$.

Given that spatial and temporal precisions are independent from each other due
to the factorization of the generalized precision matrix, and further given
the diagonal structure of the spatial precision $\Pi_{\mathbf{z}}$ (i.e.,
independence in random fluctuations across sensory sectors), we can write an
update for $\lambda_{z}$ that is a sum of squared prediction errors across
sensory sectors:

$\displaystyle\frac{d\lambda_{z}}{dt}$
$\displaystyle\approx-4\Gamma_{z}\lambda_{z}\sum_{l=1}^{L}(y^{\prime}_{h,l})^{2}-\frac{6L}{\lambda_{z}}$
(C.65)

The quadratic form of this update means that the update to the smoothness
parameter decays in proportion with the overall magnitude of the velocity
prediction errors, regardless of its sign. This means that if the distance is
fluctuating quickly in any direction, then the agent will infer that
fluctuations are slightly-less serially-correlated at the $0$th order,
reflected by a decrease in $\lambda_{z}$.

##  Appendix D Adding a target representation into the generative model

As described in the main text, it is straightforward to add an additional
observation model and dynamics model to an agent’s generative model to
represent the distance between itself and some abstract spatial target, which
in the context of the collective information transfer experiments, we
represent with T:

$\displaystyle\dot{x}_{\text{target}}$
$\displaystyle=-\alpha_{\text{target}}x_{\text{target}}+\omega_{\text{target}}$
$\displaystyle y_{\text{target}}$
$\displaystyle=x_{\text{target}}+z_{\text{target}}$
$\displaystyle\dot{x}^{\prime}_{\text{target}}$
$\displaystyle=-\alpha_{\text{target}}x^{\prime}_{\text{target}}+\omega^{\prime}_{\text{target}}$
$\displaystyle y^{\prime}_{\text{target}}$
$\displaystyle=x^{\prime}_{\text{target}}+z^{\prime}_{\text{target}}$ (D.66)

We truncate the generalized hidden states at third order
$\tilde{x}_{\text{target}}=(x_{\text{target}},x^{\prime}_{\text{target}},x^{\prime\prime}_{\text{target}})$
and the observations at second order
$\tilde{y}_{\text{target}}=(y_{\text{target}},y^{\prime}_{\text{target}})$.
When the agent assumes the generalized noises $\tilde{\omega}_{\text{target}}$
and $\tilde{z}_{\text{target}}$ are zero-mean and normally-distributed with
covariances $\tilde{\Sigma}^{\omega_{\text{target}}}$ and
$\tilde{\Sigma}^{z_{\text{target}}}$ and leverage the Laplace approximation
exactly as we did in the previous section, then we can supplement the Laplace-
approximated free energy in (A.4.1) with additional terms corresponding to
target-related prediction errors:

$\displaystyle F_{L}$
$\displaystyle\propto\frac{1}{2}\left[\tilde{\boldsymbol{\varepsilon}}_{z\text{-Soc}}^{\top}\tilde{\Pi}^{\boldsymbol{z}\text{-Soc}}\tilde{\boldsymbol{\varepsilon}}_{z\text{-Soc}}+\tilde{\boldsymbol{\varepsilon}}_{\omega\text{-Soc}}^{\top}\tilde{\Pi}^{\boldsymbol{\omega}\text{-Soc}}\tilde{\boldsymbol{\varepsilon}}_{\omega\text{-Soc}}+\tilde{\varepsilon}_{z\text{-Tar}}^{\top}\tilde{\Pi}^{z\text{-Tar}}\tilde{\varepsilon}_{z\text{-Tar}}+\tilde{\varepsilon}_{\omega\text{-Tar}}^{\top}\tilde{\Pi}^{\omega\text{-Tar}}\tilde{\varepsilon}_{\omega\text{-Tar}}\right]+C$
(D.67)

Here we use the suffixes "-Soc" or "-Tar" to indicate ‘social’ relevant
information (related to the average neighbor distance) and the ‘target’
prediction errors. $C$ captures all the additional terms (log determinants of
precision matrices, etc.) that are constant with respect to the posterior
means
$\tilde{\boldsymbol{\mu}}=(\tilde{\boldsymbol{\mu}}_{\text{Social}},\tilde{\mu}_{\text{Target}})$.
Following the same reasoning as used to derive the inference and action rules
for the case of the social distance hidden states and
observations($\tilde{\mathbf{x}}_{\text{Social}},\tilde{\mathbf{y}}_{\text{Social}})$,
we can do the same to derive active inference rules for the target-relevant
hidden states and observations
$\tilde{x}_{\text{target}},\tilde{y}_{\text{target}}$:

$\displaystyle\frac{d\tilde{\boldsymbol{\mu}}_{\text{Social}}}{dt}$
$\displaystyle=D\tilde{\boldsymbol{\mu}}_{\text{Social}}-\nabla_{\tilde{\boldsymbol{\mu}}_{\text{Social}}}F_{L}(\tilde{\boldsymbol{\mu}}_{\text{Social}},\tilde{\mathbf{y}}_{\text{Social}})$
$\displaystyle\frac{d\mathbf{v}}{dt}$
$\displaystyle=-\left(\nabla_{\mathbf{v}}F_{L}(\tilde{\boldsymbol{\mu}}_{\text{Social}},\tilde{\mathbf{y}}_{\text{Social}})+\nabla_{\mathbf{v}}F_{L}(\tilde{\mu}_{\text{Target}},\tilde{y}_{\text{Target}})\right)$
$\displaystyle\frac{d\tilde{\mu}_{\text{Target}}}{dt}$
$\displaystyle=D\tilde{\mu}_{\text{Target}}-\nabla_{\tilde{\mu}_{\text{Target}}}F_{L}(\tilde{\mu}_{\text{Target}},\tilde{y}_{\text{Target}})$
(D.68)

Where expanding the free energy gradients on the right equation leads to an
expression for the action update in terms of a precision-weighted sum of
vectors, appearing in (18) in the main text.

##  Appendix E Numerical methods

We used a forwards Euler-Maruyama scheme to the integrate a (Itô-style)
stochastic differential equation for the positions of all agents over time:

$\displaystyle d\mathbf{r}_{t}$
$\displaystyle=\mathbf{v}_{t}dt+\sigma_{a}dW_{t}$ (E.69)

where the variance of ‘action noise’ $\sigma^{2}_{a}$ was set to $0.01$ for
all experiments unless explicitly stated otherwise. We used a step size of
$\Delta t=0.01s$ in the integration. For the current timestep $\tau$ in
‘simulation time’, we used a simple forwards Euler scheme to integrate the
differential equations used for belief updating (see (A.37)) and action (see
(A.43)) for each agent in parallel. We use the positions and heading vectors
of all agents from the previous integration timestep ($\tau-\Delta t$) to
generate the observations for the current timestep.

The collective information transfer experiments were performed using custom
Julia code, and all other simulations were implemented in JAX using custom
code. To accelerate the parameter scans over $p_{inf}$,
$\Gamma_{z\text{-Social}}$, and $\Gamma_{z\text{-Target}}$ to create the
results in Figure 3 in the main text, we used the high-performance computing
clusters (Cobra and Draco) provided by the Max Planck Computing and Data
Facility.

Table E.1: Default parameter configuration used in numerical simulations (unless otherwise stated). First column denotes the name of the parameter, second column denotes its default value and third column indicates whether the parameter concerns the generative process (the physics of the simulation environment), the generative model used for active inference, or a hyperparameter (e.g., used in the active inference algorithm). Parameter | Value | Type
---|---|---
$\Delta t$ (Euler integration step, in seconds) | $0.01$ | generative process
Number of sensory sectors | $4$ | generative process
Sector angle (in degrees ∘) | $60$ | generative process
$R_{0}$ (interaction radius, in arbitrary units) | $5$ | generative process
$\sigma^{2}_{a}$ | $0.01$ | generative process
$\sigma^{2}_{z,h}$ | $0.01$ | generative process
$\sigma^{2}_{z^{\prime},h}$ | $0.01$ | generative process
$\Gamma_{z}$ | $1.0$ | generative model
$\Gamma_{\omega}$ | $1.0$ | generative model
$\lambda_{z}$ | $1.0$ | generative model
$\lambda_{\omega}$ | $1.0$ | generative model
$\alpha$ | $0.5$ | generative model
$\alpha_{\text{target}}$ | $0.5$ | generative model
$\eta$ | $1.0$ | generative model
Number of generalised coordinates ($x$) | $3$ | hyperparameter
Number of generalised coordinates ($y$) | $2$ | hyperparameter
$\kappa_{\mu}$ | $0.1$ | hyperparameter
$n_{\text{InferIter}}$ | $1$ | hyperparameter
$\kappa_{a}$ | $0.1$ | hyperparameter
$n_{\text{ActionIter}}$ | $1$ | hyperparameter
$\kappa_{\theta}$ | $0.001$ | hyperparameter
$n_{\text{LearnIter}}$ | $1$ | hyperparameter
Figure E.1: Fragmentation probability as a function of the two precision
parameters $\log\Gamma_{z}$ and $\lambda_{z}$. Fragmentation probability was
quantified as the proportion of trials (out of 500 independent trials per
condition) where the group fragmented. A trial was considered fragmented if
least one individual was further than 2.0 dimensionless units from all other
individuals for at least 3 of the last 10 seconds of the 15-second trial. All
other parameters are identical to those used in Figure 2B in the main text.

## Supplemental Movie Legends

##### Movie 1

Example of a simulation of $N=96$ agents that includes a dynamic transition
from polarized to milling regime. Parameters:
$\sigma^{2}_{z^{\prime},h}=0.05$; Sector angle $=80^{\circ}$; $R_{0}=10$
units; $\kappa_{a}=0.2$; $\lambda_{\omega}=0.5$; $\lambda_{z}=2.0$. Unless
specified, all remaining parameters are as listed in Table E.1.

##### Movie 2

Example of a polarized group of $N=64$ agents. Parameters: Sector angle
$=80^{\circ}$;; $\kappa_{a}=0.2$; $\lambda_{\omega}=0.5$; $\lambda_{z}=1.5$.

##### Movie 3

Example of a milling regime observed in $N=64$ agents. Parameters:
$\sigma^{2}_{z^{\prime},h}=0.04$; Sector angle $=80^{\circ}$; $\alpha=1.0$;
$\kappa_{a}=0.2$; $\lambda_{\omega}=0.8$; $\lambda_{z}=1.8$.

##### Movie 4

Example of a disordered regime observed in $N=96$ agents. Parameters: Number
of sensory sectors $=2$; Sector angle $=160^{\circ}$; $R_{0}=10.0$ units;
$\alpha=0.2$; $\eta=0.5$, $\kappa_{a}=0.2$; $\lambda_{\omega}=0.1$;
$\lambda_{z}=1.787$.

##### Movie 5

Metastable ‘snaking’ configuration observed in $N=64$ agents. Parameters:
$\sigma^{2}_{z^{\prime},h}=0.04$; Sector angle $=80^{\circ}$; $\alpha=0.1$;
$\kappa_{a}=0.2$; $\lambda_{\omega}=0.5$; $\lambda_{z}=2.2$.