* Kowdle et al (2012) Kowdle A, Batra D, Chen WC, Chen T (2012) imodel: Interactive co-segmentation for object of interest 3D modeling. In: Proceedings of the 11th European Conference on Trends and Topics in Computer Vision - Volume Part II, Springer-Verlag, ECCV’10, pp 211–224, DOI 10.1007/978-3-642-35740-4˙17, URL http://dx.doi.org/10.1007/978-3-642-35740-4_17
  * Kuang and Åström (2013) Kuang Y, Åström K (2013) Pose estimation with unknown focal length using points, directions and lines. In: International Conference on Computer Vision, IEEE, pp 529–536
  * Kuang et al (2014) Kuang Y, Oskarsson M, Åström K (2014) Revisiting trifocal tensor estimation using lines. In: Pattern Recognition (ICPR), 2014 22nd International Conference on, IEEE, pp 2419–2423
  * Kunsberg and Zucker (2014) Kunsberg B, Zucker SW (2014) Why shading matters along contours. In: Neuromathematics of Vision, Lecture Notes in Morphogenesis, Springer, pp 107–129, DOI 10.1007/978-3-642-34444-2˙3, URL http://dx.doi.org/10.1007/978-3-642-34444-2_3
  * Lebeda et al (2014) Lebeda K, Hadfield S, Bowden R (2014) 2D or not 2D: Bridging the gap between tracking and structure from motion. In: Proc. of ACCV
  * Li and Zucker (2003) Li G, Zucker SW (2003) A differential geometrical model for contour-based stereo correspondence. In: Proc. IEEE Workshop on Variational, Geometric, and Level Set Methods in Computer Vision, Nice, France
  * Li and Zucker (2006) Li G, Zucker SW (2006) Contextual inference in contour-based stereo correspondence. International Journal of Computer Vision 69(1):59–75
  * Lin et al (2009) Lin WY, Tan GC, Cheong LF (2009) When discrete meets differential. International Journal of Computer Vision 86(1):87–110, DOI 10.1007/s11263-009-0260-y, URL http://dx.doi.org/10.1007/s11263-009-0260-y
  * Litvinov et al (2012) Litvinov V, Yu S, Lhuillier M (2012) 2-manifold reconstruction from sparse visual features. In: 3D Imaging (IC3D), 2012 International Conference on, IEEE, pp 1–8
  * Longuet-Higgins (1981) Longuet-Higgins HC (1981) A computer algorithm for reconstructing a scene from two projections. Nature 293:133–135
  * Longuet-Higgins and Prazdny (1980) Longuet-Higgins HC, Prazdny K (1980) The interpretation of a moving retinal image. Proc R Society London 208(1173):385–397
  * Lowe (2004) Lowe DG (2004) Distinctive image features from scale-invariant keypoints. Int Journal of Computer Vision 60(2):91–110
  * Ma et al (2004) Ma Y, Soatto S, Kosecka J, Sastry SS (2004) An invitation to 3D vision. Springer
  * Mattingly et al (2015) Mattingly WA, Chariker JH, Paris R, jen Chang D, Pani JR (2015) 3D modeling of branching structures for anatomical instruction. Journal of Visual Languages & Computing 29(0):54 – 62, DOI http://dx.doi.org/10.1016/j.jvlc.2015.02.006, URL http://www.sciencedirect.com/science/article/pii/S1045926X15000191
  * Maybank (1992) Maybank S (1992) Theory of Reconstruction from Image Motion. Springer-Verlag New York, Inc., Secaucus, NJ, USA
  * Mendonça et al (2001) Mendonça PRS, Wong KYK, Cipolla R (2001) Epipolar geometry from profiles under circular motion. IEEE Trans Pattern Anal Mach Intell 23(6):604–616, DOI http://dx.doi.org/10.1109/34.927461
  * Mikolajczyk and Schmid (2004) Mikolajczyk K, Schmid C (2004) Scale and affine invariant interest point detectors. IJCV 60(1):63–86
  * Mikolajczyk and Schmid (2005) Mikolajczyk K, Schmid C (2005) A performance evaluation of local descriptors. IEEE Trans Pattern Anal Mach Intell 27(10):1615–1630
  * Moravec (1977) Moravec HP (1977) Towards automatic visual obstacle avoidance. In: Proc. of the 5th International Joint Conference on Artificial Intelligence, p 584
  * Moreels and Perona (2007) Moreels P, Perona P (2007) Evaluation of features detectors and descriptors based on 3d objects. Int J Comput Vision 73(3):263–284, DOI http://dx.doi.org/10.1007/s11263-006-9967-1
  * Ohta and Kanade (1985) Ohta Y, Kanade T (1985) Stereo by intra- and inter-scanline search using dynamic programming. IEEE Trans Pattern Analysis and Machine Intelligence 7(2):139–154
  * Papadopoulo (1996) Papadopoulo T (1996) Motion analysis of 3D rigid curves from monocular image sequences. Tech. Rep. RR–2779 (PhD Thesis), INRIA
  * Papadopoulo and Faugeras (1996) Papadopoulo T, Faugeras OD (1996) Computing structure and motion of general 3D curves from monocular sequences of perspective images. In: Proceedings of the $4^{th}$ European Conference on Computer Vision, Springer-Verlag, London, UK, pp 696–708
  * Pollefeys et al (2004) Pollefeys M, Van Gool L, Vergauwen M, Verbiest F, Cornelis K, Tops J, Koch R (2004) Visual modeling with a hand-held camera. International Journal of Computer Vision 59(3):207–232, DOI http://dx.doi.org/10.1023/B:VISI.0000025798.50602.3a
  * Ponce and Genc (1998) Ponce J, Genc Y (1998) Epipolar geometry and linear subspace methods: a new approach to weak calibration. International Journal of Computer Vision 28(3):223–243
  * Porrill and Pollard (1991) Porrill J, Pollard S (1991) Curve matching and stereo calibration. Image and Vision Computing 9(1):45–50, DOI http://dx.doi.org/10.1016/0262-8856(91)90048-T
  * Pötsch and Pinz (2011) Pötsch K, Pinz A (2011) 3D geometric shape modeling by ‘3D contour cloud’ reconstruction from stereo videos. In: 16th Computer Vision Winter Workshop, Citeseer, p 99
  * Rao et al (2012) Rao D, Chung SJ, Hutchinson S (2012) CurveSLAM: An approach for vision-based navigation without point features. In: IEEE/RSJ Intelligent Robots and Systems (IROS), pp 4198–4204, DOI 10.1109/IROS.2012.6385764
  * Restrepo et al (2014) Restrepo MI, Ulusoy AO, Mundy JL (2014) Evaluation of feature-based 3-D registration of probabilistic volumetric scenes. ISPRS Journal of Photogrammetry and Remote Sensing 98:1–18
  * Reyes and Bayro Corrochano (2005) Reyes L, Bayro Corrochano E (2005) The projective reconstruction of points, lines, quadrics, plane conics and degenerate quadrics using uncalibrated cameras. Image and Vision Computing 23(8):693–706
  * Robert and Faugeras (1991) Robert L, Faugeras OD (1991) Curve-based stereo: figural continuity and curvature. In: Proceedings of Computer Vision and Pattern Recognition, pp 57–62
  * Schmid and Zisserman (2000) Schmid C, Zisserman A (2000) The geometry and matching of lines and curves over multiple views. International Journal of Computer Vision 40(3):199–233, DOI http://dx.doi.org/10.1023/A:1008135310502
  * Schneevoigt et al (2014) Schneevoigt T, Schroers C, Weickert J (2014) A dense pipeline for 3d reconstruction from image sequences. In: Pattern Recognition, Springer, pp 629–640
  * Seitz et al (2006) Seitz S, Curless B, Diebel J, Scharstein D, Szeliski R (2006) A comparison and evaluation of multi-view stereo reconstruction algorithms. In: Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, IEEE Computer Society, pp 519–528
  * Shashua (1994) Shashua A (1994) Trilinearity in visual recognition by alignment. In: Proceedings of the third European conference on Computer vision, Springer-Verlag, Secaucus, NJ, USA, pp 479–484
  * Sherman and Peleg (1990) Sherman D, Peleg S (1990) Stereo by incremental matching of contours. IEEE Trans on Pattern Analylis and Machine Intelligence 12(11):1102–1106, DOI http://dx.doi.org/10.1109/34.61711
  * Shinozuka and Saito (2014) Shinozuka Y, Saito H (2014) Sharing 3D object with multiple clients via networks using vision-based 3D object tracking. In: Proceedings of the 2014 Virtual Reality International Conference, ACM, New York, NY, USA, VRIC ’14, pp 34:1–34:4, DOI 10.1145/2617841.2620723, URL http://doi.acm.org/10.1145/2617841.2620723
  * Simoes et al (2012) Simoes F, Almeida M, Pinheiro M, dos Anjos R (2012) Challenges in 3D reconstruction from images for difficult large-scale objects. 2014 XVI Symposium on Virtual and Augmented Reality 0:74–83, DOI http://doi.ieeecomputersociety.org/10.1109/SVR.2012.5
  * Sinha et al (2004) Sinha SN, Pollefeys M, McMillan L (2004) Camera network calibration from dynamic silhouettes. In: Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition, IEEE Computer Society, pp 195–202
  * Spetsakis and Aloimonos (1991) Spetsakis M, Aloimonos JY (1991) A multi-frame approach to visual motion perception. Int J Comput Vision 6(3):245–255, DOI http://dx.doi.org/10.1007/BF00115698
  * Stewénius et al (2007) Stewénius H, Engels C, Nistér D (2007) An efficient minimal solution for infinitesimal camera motion. In: Computer Vision and Pattern Recognition, 2007\. CVPR’07. IEEE Conference on, IEEE, pp 1–8
  * Tamrakar (2008) Tamrakar A (2008) Image contour extraction using geometric consistency. Ph.D. dissertation, Division Of Engineering, Brown University, Providence, RI, 02912
  * Tamrakar and Kimia (2007) Tamrakar A, Kimia BB (2007) No grouping left behind: From edges to curve fragments. In: ICCV2007 (2007)
  * Teney and Piater (2012) Teney D, Piater J (2012) Sampling-based multiview reconstruction without correspondences for 3D edges. In: Proceedings, IEEE, Washington, DC, USA, 3DIMPVT ’12, pp 160–167, DOI 10.1109/3DIMPVT.2012.28, URL http://dx.doi.org/10.1109/3DIMPVT.2012.28
  * Tian et al (1996) Tian TY, Tomasi C, Heeger DJ (1996) Comparison of approaches to egomotion computation. In: Computer Vision and Pattern Recognition, 1996. Proceedings CVPR’96, 1996 IEEE Computer Society Conference on, IEEE, pp 315–320
  * Triggs (1999) Triggs B (1999) Differential Matching Constraints. In: Proceedings of the IEEE International Conference on Computer Vision, IEEE Computer Society Press, KerKyra, Greece
  * Usumezbas et al (2016) Usumezbas A, Fabbri R, Kimia BB (2016) From multiview image curves to 3d drawings. In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, IEEE Computer Society Press, Las Vegas, Nevada, USA, submitted
  * Valgaerts et al (2012) Valgaerts L, Bruhn A, Mainberger M, Weickert J (2012) Dense versus sparse approaches for estimating the fundamental matrix. International Journal of Computer Vision 96(2):212–234, DOI 10.1007/s11263-011-0466-7, URL http://dx.doi.org/10.1007/s11263-011-0466-7
  * Vieville and Faugeras (1996) Vieville T, Faugeras O (1996) The First Order Expansion of Motion Equations in the Uncalibrated Case. Computer Vision and Image Understanding 64(1):128–146
  * Viéville and Faugeras (1995) Viéville T, Faugeras OD (1995) Motion analysis with a camera with unknown, and possibly varying intrinsic parameters. In: Computer Vision, 1995. Proceedings., Fifth International Conference on, IEEE, pp 750–756
  * Wang et al (2014) Wang R, Choi J, Medioni G (2014) 3D modeling from wide baseline range scans using contour coherence. In: Computer Vision and Pattern Recognition (CVPR), 2014 IEEE Conference on, IEEE, pp 4018–4025
  * Waxman and Ullman (1985) Waxman AM, Ullman S (1985) Surface structure and three-dimensional motion from image flow kinematics. International Journal of Robotics Research 4(3):72–94
  * Wong and Cipolla (2004) Wong KY, Cipolla R (2004) Reconstruction of sculpture from its profiles with unknown camera positions. IEEE Transactions on Image Processing 13(3):381–9
  * Wong et al (2001) Wong KYK, Mendonça PRS, Cipolla R (2001) Head model acquisition from silhouettes. In: Arcelli C, Cordella L, di Baja GS (eds) IWVF, Springer, Capri, Italy, pp 787–796
  * Yi Ma (1998) Yi Ma SS Jana Koseck (1998) Motion recovery from image sequences: Discrete viewpoint vs. differential viewpoint. In: Proceedings of European Conference on Computer Vision, Springer, Lecture Notes in Computer Science, vol 1407, p 337
  * Zhang (2013) Zhang L (2013) Line primitives and their applications in geometric computer vision. PhD thesis, Department of Computer Science, Univ.
  * Zhuang and Haralick (1984) Zhuang X, Haralick RM (1984) Rigid body motion and the optic flow image. In: Proceedings of the First Conference on Artificial Intelligence Applications, IEEE Computer Society, pp 366–375
  * Zucker (2014) Zucker S (2014) Stereo, shading, and surfaces: Curvature constraints couple neural computations. Proceedings of the IEEE 102(5):812–829


# Self-consistent stochastic dynamics for finite-size networks of spiking
neurons

Gianni V. Vinci  Natl. Center for Radiation Protection and Computational
Physics, Istituto Superiore di Sanità, 00169 Roma, Italy  PhD Program in
Physics, “Tor Vergata” University of Rome, 00133 Roma, Italy  Roberto Benzi
Dept. of Physics and INFN, “Tor Vergata” University of Rome, 00133 Roma, Italy
Maurizio Mattia<EMAIL_ADDRESS>Natl. Center for Radiation Protection
and Computational Physics, Istituto Superiore di Sanità, 00169 Roma, Italy

(December 29, 2021)

###### Abstract

Despite the huge number of neurons composing a brain network, ongoing activity
of local cell assemblies composing cortical columns is intrinsically
stochastic. Fluctuations in their instantaneous rate of spike firing $\nu(t)$
scale with the size of the assembly and persist in isolated network, i.e., in
absence of external source of noise. Although deterministic chaos due to the
quenched disorder of the synaptic couplings likely underlies this seemingly
stochastic dynamics, an effective theory for the network dynamics of a finite
ensemble of spiking neurons is lacking. Here, we fill this gap by extending
the so-called population density approach including an activity- and size-
dependent stochastic source in the Fokker-Planck equation for the membrane
potential density. The finite-size noise embedded in this stochastic partial
derivative equation is analytically characterized leading to a self-consistent
and non-perturbative description of $\nu(t)$ valid for a wide class of spiking
neuron networks. Its power spectra of $\nu(t)$ are found in excellent
agreement with those from detailed simulations both in the linear regime and
across a synchronization phase transition, when a size-dependent smearing of
the critical dynamics emerges.

## Finite-size physical systems. —

Complex systems in statistical physics usually deal with a huge number $N$ of
interacting bodies like molecules leading to have wide applicability of
effective mean-field theories. However, the finite size (e.g., volume) of the
system can have a significant impact in the emergent collective dynamics [39].
This can be due to the fact that finite systems have specific boundaries
(finite surfaces and volumes) affecting the system behavior when phase
transitions are approached [7, 20]. Besides, the presence of a finite number
of elements can be incorporated as a stochastic field whose fluctuation size
depend on $N$ [21, 10]. In this way, the continuum formalism valid in the
thermodynamic limit may result to be effective in the description of small-
scale phenomena across phase transitions [12]. Biological and ecological
systems can be even more challenging as they may incorporate peculiar boundary
conditions and heterogeneity in space and time posing them continuously
outside equilibrium [31, 19]. This is the case of biological networks of
neurons for which the probability current (see below) does not vanish even
under stationary condition when it matches the frequency of spikes emitted per
neuron, i.e. the firing rate $\nu$ [22, 1, 4, 14]. In this Letter we will show
that finite-size fluctuations can be effectively taken into account in this
challenging system via a self-consistent definition of the noise to be
embedded into the mean-field population dynamics.

## Population density and mean-field approximation. —

In the thermodynamic limit ($N\to\infty$) networks of single-compartment
spiking neurons have collective dynamics statistically described by the
probability density $p(v,t)$ of realizations with membrane potential $v$ at
time $t$ [1, 23, 4, 14], obeying the Fokker-Planck equation

$\partial_{t}p=-\partial_{v}S_{p}=-\partial_{v}[(F+\mu)p]+\frac{1}{2}\partial^{2}_{v}(\sigma^{2}p)\,.$
(1)

In this continuity equation the density $p$ changes according to the
divergence of the probability current
$S_{p}(v,t)=(F+\mu)p+\frac{1}{2}\partial_{v}(\sigma^{2}p)$. Here, the membrane
potential $V(t)$ of a neuron follows the nonlinear Langevin equation

$dV=\left[F(V)+\mu\right]dt+\sigma dW\,,$ (2)

where $F(V)$ is the drifting current determining the model-specific relaxation
dynamics [17]. In what follows, as workbench we adopt the ‘leaky’ integrate-
and-fire neuron with $F(V)=-V/{\tau}$ and decay time ${\tau}$. The total
synaptic and ionotropic input current is a Gaussian white noise $W(t)$
[$\langle W(t)W(t^{\prime})\rangle=\delta(t-t^{\prime})$] inhomogeneously
modulated in time to have infinitesimal mean $\mu(V,t)$ and variance
$\sigma^{2}(V,t)$. Such diffusion approximation holds in the limit of large
number $K$ of presynaptic contacts and small average synaptic efficacy $J$
[38], well describing cortical networks [11, 25]. Following the Grigelionis
central limit theorem [18], this is due to having pooled together presynaptic
spike trains with arbitrary inter-spike interval statistics provided that the
$K$ input neurons are independent and fire at low rates $\nu$, eventually
giving rise to a inhomogeneous Poisson process. The independence hypothesis
guaranteed by having small $J$, underlies also the mean-field approximation.
In this framework all neurons are independent realizations of the same
stochastic process (same $F$, $\mu$ and $\sigma$) [2]. Synaptic interactions
are incorporated in the $\nu$-dependent moments of the current, which in
current-based models like LIF neurons are

$\begin{array}[]{rcl}\mu(\nu)&=&KJ\nu(t)+\mu_{\mathrm{ext}}\\\
\sigma^{2}(\nu)&=&KJ^{2}\nu(t)+\sigma^{2}_{\mathrm{ext}}\end{array}\,.$ (3)

Here, the firing rate $\nu(t)=S_{p}({v_{\mathrm{thr}}},t)$ is the flux of
neurons crossing the emission threshold ${v_{\mathrm{thr}}}$, and
$\mu_{\mathrm{ext}}$ and $\sigma^{2}_{\mathrm{ext}}$ are the moments of the
synaptic current due to the spikes incoming from neurons external to the
network.

## Why are finite-size fluctuations important? —

Having a finite number $N$ of neurons emitting spike trains
$S_{i}(t)=\sum_{k}{\delta(t-t_{i,k})}$, leads to a fluctuating rate
$\nu_{N}(t)=\sum_{i=1}^{N}{S_{i}(t)}/N=\mathcal{N}(t)$ of action potentials
fired per unit time and per neuron. According to the aforementioned central
limit, $\mathcal{N}(t)$ is an inhomogeneous Poission process with mean
$N\,\nu(t)$, such that $\lim_{N\to\infty}\nu_{N}(t)=\nu(t)$. Infinitesimal
mean and variance of $\nu_{N}$ are $\nu$ and $\nu/N$, respectively, and for
large enough $N$ ($\sim 100$) like those observed in cortical minicolumns
[30], holds the Gaussian approximation

$\nu_{N}(t)=\nu(t)+\sqrt{\frac{\nu(t)}{N}}W(t)\equiv\nu(t)+\eta(t)\,.$ (4)

Here, the finite-size noise $\eta(t)$ results from a modulation of a white
noise $W(t)$ independent from the one in Eq. (2). The correlation structure of
$\eta$ will be self-consistently derived in the following. Taking into account
the fluctuating $\nu_{N}$ into Eq. (3), the moments of the input current
(i.e., the mean-field) are no longer deterministic [4]. Fluctuations of the
mean $\mu$ have relative size
$\mathrm{Var}[\mu(\nu_{N})]^{1/2}/\mathrm{E}[\mu(\nu_{N})]=1/\sqrt{N\nu}$ of
about 10% in cortical networks of interest where $\nu\sim 1$ Hz [42].
Remarkably, such variability is of the same order of the changes in $\mu$
associated with the coding of sensorial stimuli [32] or other relevant
information [33], which usually involves only a sparse set of tuned neurons
with variations of few Hz in their firing rates. Thus, finite-size
fluctuations may have a not negligible impact disturbing or nonlinearly
amplifying the encoding dynamics of cortical networks.

## How can finite-size fluctuations be incorporated? —

To understand the impact of such fluctuations, according to [4] we incorporate
the stochastic moments $\mu_{N}\equiv\mu(\nu_{N})$ and
$\sigma_{N}\equiv\sigma(\nu_{N})$ directly into Eq. (1). The resulting
stochastic Fokker-Planck equation describes now an infinite set of independent
neurons all driven by the same fluctuating mean-field. Differently from other
stochastic Smoluchowski equations [21, 10] and their coarse-grained versions
[5, 6], stochasticity here appears as an additional probability current with
drift and diffusion coefficients differently affected by $\eta(t)$. Not only,
an additional source of noise must be incorporated as a fluctuating source of
realizations in $v={v_{\mathrm{res}}}$ [27]. This is due to the fact that the
finite flux of neurons crossing the threshold ${v_{\mathrm{thr}}}$ reenters at
the reset potential, eventually leading to

$\partial_{t}p=-\partial_{v}[(F+\mu_{N})p]+\frac{1}{2}\partial^{2}_{v}(\sigma_{N}^{2}p)+\delta(v-{v_{\mathrm{res}}})\,\eta(t)\,.$
(5)

We remark that this stochastic Fokker-Planck (SFP) equation is nonlinear as
both $\mu_{N}$ and $\sigma_{N}$ depend on the density $p$ via the firing rate
$\nu$.

## Self-consistent derivation of $\eta(t)$. —

To determine the statistical features of $\eta(t)$ we refer to the specific
case of a set of $N$ uncoupled neurons ($J=0$) driven by a stationary external
input ($\dot{\mu}_{\mathrm{ext}}=\dot{\sigma}_{\mathrm{ext}}=0$). In this case
neurons are renewal processes, and the probability density $\rho(t)$ of their
inter-spike intervals (ISI) fully characterize the statistics of the spike
trains they emit [8]. Pooling together these spike trains gives $\nu_{N}(t)$
(see above) such that its power spectral density is [24, 17]

$P_{\nu}^{(\mathrm{RT})}(\omega)=|\nu_{N}(\omega)|^{2}=\frac{\nu_{0}}{N}\,\mathrm{Re}\left[\frac{1+\rho(\omega)}{1-\rho(\omega)}\right]\,.$
(6)

Here $f(\omega)=\int_{-\infty}^{\infty}{f(t)e^{-i\omega t}dt}$ is the Fourier
transform of any function $f(t)$, and $\nu_{0}=1/\langle T\rangle$ is the mean
firing rate, i.e., the inverse of the mean ISI.

This exact result must be also obtained from Eq. (5). To carry out the power
spectral density $|\nu_{N}(\omega)|^{2}$ in this case we resort to the
spectral expansion approach introduced in [27] giving

$P_{\nu}^{(\mathrm{SE})}(\omega)=\,\left|1+\vec{f}\cdot(i\omega\mathbf{I}-\mathbf{\Lambda})^{-1}{\vec{\psi}_{\mathrm{res}}}\right||\eta(\omega)|^{2}$
(7)

If the finite-size noise is assumed to be white,
$|\eta(\omega)|^{2}=\nu_{0}/N$ leading to an overestimate of the power
$P_{\nu}(\omega)$ at relatively low-$\omega$ at least under ‘suprathreshold’
regime, i.e., when neurons emit spikes even if $\sigma=0$ [27].

In Eq. (7), we make use of the eigenfunctions $\phi_{n}(v)$ of the non-
Hermitian Fokker-Planck operator $\mathcal{L}$ defined from Eq. (1) as
$\partial_{t}p\equiv\mathcal{L}\,p$. This operator has an infinite spectrum of
discrete eigenvalues $\lambda_{n}$ ($n\in\mathbb{Z}$) such that
$\mathcal{L}\,\phi_{n}=\lambda_{n}\,\phi_{n}$ [1, 23, 27]. The flux
$f_{n}=1/2\partial_{v}(\sigma\phi_{n})|_{v={v_{\mathrm{thr}}}}$ of
nonstationary ($n\neq 0$) eigenfunctions are the infinite elements of
$\vec{f}$, while the matrix $\mathbf{\Lambda}$ is diagonal with
$\\{\mathbf{\Lambda}\\}_{nn}=\lambda_{n}$. The elements of
${\vec{\psi}_{\mathrm{res}}}$ are instead the eigenfunction $\psi_{n}(v)$ of
the adjoint operator $\mathcal{L}^{\dagger}$ with same $\lambda_{n}$, computed
in $v={v_{\mathrm{res}}}$ [23, 27]. Note that all these coefficients are
state-depend being functions of the current moments $\mu$ and $\sigma$.

Figure 1: Normalized power spectra of the finite-size noise $\eta$ for $N$
uncoupled LIF neurons with firing rate $\nu_{0}$. (a)
$|\eta(\omega)|^{2}N/\nu_{0}$ as a function of the mean synaptic current
$\mu$. $\sigma$ is chosen to keep the mean firing rate unchanged at
$\nu_{0}=20$ Hz. (b) In noise-dominate regime
($\mu{\tau}<{v_{\mathrm{thr}}}=20$ mV) finite-size noise is essentially white
(bottom). In drift-dominate regime, power is low at low-$\omega$ as ISIs are
more regular ($c_{v}<1$). In both cases, a two-dimensional Markovian embedding
(black) faithfully reproduce theoretical spectra (red) from Eq. (9).

For any IF neuron model, the series in Eq. (7) can be summed as a function of
$\rho(\omega)$ [41]:

$\vec{f}\cdot(i\omega\mathbf{I}-\mathbf{\Lambda})^{-1}{\vec{\psi}_{\mathrm{res}}}=\frac{\rho(\omega)}{1-\rho(\omega)}-\frac{\nu_{0}}{i\omega}\,.$
(8)

Using it in Eq. (7) and requiring the equivalence
$P_{\nu}^{(\mathrm{SE})}(\omega)=P_{\nu}^{(\mathrm{RT})}(\omega)$, we obtain
with Eq. (6) a self-consistent expression for the power spectrum of $\eta$:

$|\eta(\omega)|^{2}=\frac{\nu_{0}}{N}\left[1-\left|\frac{(i\omega+\nu_{0})\rho(\omega)-\nu_{0}}{\nu_{0}\rho(\omega)+i\omega-\nu_{0}}\right|^{2}\right]\,.$
(9)

In the limit $\omega\to\infty$, $\rho(\omega)\to 0$ and the finite-size noise
is white-like with variance $\nu_{0}/N$. For $\omega\to 0$, the l.h.s. of Eq.
(8) reduces to
$-\vec{f}\cdot\mathbf{\Lambda}^{-1}{\vec{\psi}_{\mathrm{res}}}=(c_{v}^{2}-1)/2$
[41], function of the coefficient of variation $c_{v}$ of the ISIs, leading to

$|\eta(0)|^{2}=\frac{\nu_{0}}{N}\frac{4c_{v}^{2}}{(1+c_{v}^{2})^{2}}\,.$ (10)

These limits suggest a sigmoid-like power spectra of the finite-size noise
which is confirmed in LIF neurons [Fig. 1(a)]. As the firing regimes moves
from noise- to drift-dominated (i.e., from sub- to suprathreshold) regimes by
increasing the mean current $\mu$, the sigmoidal shape becomes increasingly
more apparent. Indeed, ISIs are more and more regular leading to a decrease of
their $c_{v}$, and hence to a lower power $|\eta(0)|^{2}$.

## Markovian embedding of $\eta(t)$. —

Having derived in Eq. (9) the correlation structure of the finite-size noise
for a set of uncoupled neurons under stationary condition a question arises:
can we generalize this result to networks of synaptically coupled neurons and
far from equilibrium?

To answer this question we remark that the spectra from Eq. (9) shown in Fig.
1 are constant (white noise) with a power reduction at low-$\omega$ possibly
resulting by subtracting a Lorentzian-shaped function. Following [40], such
kind of spectra are very well approximated by two-dimensional Ornstein-
Uhlenbeck processes $\vec{u}(t)$ driven by and interfering with the same white
noise $W(t)$:

$\begin{array}[]{rcl}d\vec{u}&=&\mathbf{A}\,\vec{u}\,dt+\mathbf{B}\,dW\\\
\eta&=&\displaystyle\vec{1}\cdot\vec{u}+\sqrt{\frac{\nu_{0}}{N}}\,W\end{array}\,.$
(11)

In Fig. 1(b) we show that this Markovian embedding faithfully reproduces the
correlation structure of $\eta(t)$ in the simple case of $\mathbf{B}$ with
only one non-zero element, $B_{11}=\sqrt{\nu_{0}/N}$, and of a three-parameter
$\mathbf{A}$ defined as $A_{11}=A_{22}=a$, $A_{12}=a_{\uparrow}$ and
$A_{21}=a_{\downarrow}$. Instead of fitting
$\\{a,a_{\uparrow},a_{\downarrow}\\}$, we carry them out analytically as in
[40] by matching the power of $\eta$ from Eq. (11) with the exact one in Eq.
(9) at the frequencies: $\omega=\\{0,\pi\nu_{0},2\pi\nu_{0}\\}$. The
parameters change according to $\mu$ and $\sigma$ leading to a state-dependent
$\mathbf{A}(\mu,\sigma)$.

This dynamical description of $\eta(t)$ in principle allows to overcome the
renewal hypothesis, as the memory embedded in the network activity is
reintroduced via the dependence on $\nu(t)$ of the current moments in Eq. (3).
We then conjecture that taken together Eqs. (5) and (11) provide a complete
statistical description of the out-of-equilibrium dynamics of a finite-size
network of spiking neurons.

## Interacting neurons and linear-response theory. —

Figure 2: Bifurcation analysis of a network of LIF excitatory neurons
($N\to\infty$) at varying synaptic efficacy $J$. (a) Real and imaginary parts
of the poles $\mathrm{sd}_{\pm 1}$ and $\mathrm{st}_{1}$ of $\nu(s)$ solving
Eq. (12). Weakly coupled networks (small $KJ$, blue range) have a stable focus
in $\nu_{0}$. Beyond $KJ\simeq 11$ mV (red) limit cycles arise via a
supercritical Hopf bifurcation, eventually crossing a saddle-node bifurcation
(grey). (b) Numerical integration of the Fokker-Planck equation with
$p(v,0)=\delta(v-{v_{\mathrm{res}}})$ for $KJ$ pointed out in (a) (top
arrows). The persistence of the fixed-point $\nu_{0}=20$ Hz results from
keeping constant the moments $\mu\,{\tau}=21$ mV and
$\sigma\,{\tau}^{1/2}=2.665$ mV, and changing $\mu_{\mathrm{ext}}$ and
$\sigma_{\mathrm{ext}}$ according to $J$. Other parameters: $K=10^{3}$,
${\tau}=20$ ms, ${v_{\mathrm{thr}}}=20$ mV, $\delta_{\mathrm{min}}=2$ ms and
${\tau_{\mathrm{\delta}}}=1$ ms.

We test this conjecture using as test bench a network of LIF neurons with a
fixed-point at $\nu_{0}=20$ Hz in the thermodynamic limit ($N\to\infty$). We
balance the increase in the excitatory synaptic efficacy $J$ varying the
moments $\mu_{\mathrm{ext}}$ and $\sigma_{\mathrm{ext}}$ of the external
current to keep $\nu_{0}$ unchanged. An exponential distribution of axonal
delays $\delta>\delta_{\mathrm{min}}$ in delivering emitted spikes is taken
into account and having the same impact as non-instantaneous synaptic
transmission [26]. The incoming rate of spikes $\tilde{\nu}(t)$ has dynamics
$\dot{\tilde{\nu}}=(\nu-\tilde{\nu})/{\tau_{\mathrm{\delta}}}$, and replaces
$\nu(t)$ in Eq. (3). The stability analysis of the equilibrium point $\nu_{0}$
of the network can be effectively carried out resorting to the spectral
expansion of $p(v,t)$ mentioned above. Indeed, the linear dynamics of the
perturbations $\nu(t)-\nu_{0}$ is fully captured by the time course of the
projections on the two slowest eigenmodes [29]. Stability are then determined
by the poles $s$ of the Laplace transform
$\nu(s)=\int_{0}^{\infty}{e^{-st}\nu(t)dt}$ of the firing rate solving the
equation [27]

$1-g(s)\left[\Phi^{\prime}+\frac{s}{{\tau}}\left(\frac{c_{1}}{s-\lambda_{1}}+\frac{c_{2}}{s-\lambda_{2}}\right)\right]=0\,.$
(12)

Here $g(s)=\frac{e^{-(s-\delta_{\mathrm{min}})/\tau_{\delta}}}{\tau_{\delta}}$
is the Laplace transform of the delay distribution,
$\Phi^{\prime}=\partial_{\nu}\Phi$ is the slope of the current-to-rate gain
function $\Phi(\mu,\sigma)=S_{\phi_{0}}({v_{\mathrm{thr}}})$ computed in
$\nu_{0}$, and
$c_{n}=\int_{v_{\mathrm{min}}}^{v_{\mathrm{thr}}}{\phi_{0}\partial_{\nu}\psi_{n}dv}$
with $n\in\\{1,2\\}$ are coupling coefficients proportional to $J$ [27, 29].
Real and imaginary parts of the three poles solving Eq. (12) are shown in Fig.
2(a) as a function of $KJ$ . Three phases are singled out starting from a
stable focus in $\nu_{0}$ for weak couplings which with increasing $KJ$
destabilizes via a supercritical Hopf bifurcation giving rise to a stable
limit cycle. Large enough $KJ$ eventually lead to a saddle-node bifurcation
where $\nu(t)$ is asymptotically attracted to low firing rates. The numerical
integrations [3] of the Fokker-Planck Eq. (1) for three representative
networks are shown in Fig. 2(b) confirming the reliability of the approximated
bifurcation analysis.

## Linear response regime for finite networks. —

Figure 3: Power spectra $P_{\nu}(\omega)$ (a) and firing rates $\nu_{N}(t)$
(b) of weakly coupled networks ($KJ=5$ mV) with varying size $N$. (a) Top,
$P_{\nu}(\omega)$ from the numerical integration of the stochastic Fokker-
Planck equation (SFP, pink), the equivalent spiking neuron network (SNN)
simulation (black) and linear perturbation theory (dashed blue) ($N=10^{4}$).
Dashed black, white noise with variance $\nu_{0}/$. Bottom, normalized spectra
of $\nu$ from SFP of networks with different $N$ ($=10^{3},10^{4},10^{5}$).
Theoretical power spectrum (dotted blue) is from Eq. (13). (b) $\nu(t)$ from
SNN simulations (black) and SFP varying $N$ ($=10^{3},10^{4},10^{5}$) (left).
Right, related histograms of $\nu(t)$. Other parameters as in Fig. 2.

In weakly coupled networks with relatively small $KJ$, finite-size noise
induces stochastic perturbations of $\nu_{N}(t)$ around an equilibrium point
amenable to linear response theory. This is the case of our example network
with $KJ=5$ mV for which the power spectrum of $\nu_{N}$ can be carried out as
a linear transformation of the finite-size noise [27, 28]:

$P_{\nu}(\omega)=|\nu_{N}(\omega)|^{2}=|H_{\nu}(\omega)|^{2}|\eta(\omega)|^{2}\,.$
(13)

Here $H_{\nu}(\omega)$ is the Fourier transfer function giving the linear
response to sinusoidal modulations of the input firing rate which is
analytically known for LIF neurons [4, 41]. In Fig. 3(a)-top we show a
remarkable agreement between such theoretical expression, and
$P_{\nu}(\omega)$ estimated from both the simulations of spiking neuron
networks (SNN, [13]) and the numerical integration of the SFP Eq. (5). To
computed the SFP integration we extended a standard deterministic approach [3]
by incorporating the Markovian embedding (11) of $\eta(t)$. SFP integration
for different network sizes $N$ confirms in Fig. 3(a)-bottom what expected
from the linear theory, that is, the spectrum shape does not change with $N$
once normalized by the variance $\nu_{0}/N$. The equivalence in these case
between SNN simulation and SFP integration is also apparent in the direct
comparison of $\nu_{N}(t)$ time series. Indeed, in Fig. 3(b) no differences
emerge and both time series display a variance scaling as $1/N$.

## Beyond linear response theory. —

Figure 4:  Power spectra $P_{\nu}(\omega)$ (a,c) and firing rates $\nu_{N}(t)$
(b,d) of strongly coupled (nonlinear) networks with varying size $N$. (a, b)
$KJ=10$ mV, mean-field dynamics predicts a stable focus at $\nu_{0}=20$ Hz.
(c, d) $KJ=12$ mV, networks beyond a supercritical Hopf bifurcation with a
stable limit cycle. See Figs. 2 and 3 for additional details.

By further increasing synaptic coupling $KJ$, a supercritical Hopf bifurcation
is approached (Fig. 2) and finite-size networks have $P_{\nu}(\omega)$ no
longer fully described by linear theory [i.e., Eq. (13)]. The mismatch is
apparent in Fig. 4(a)-top as an unpredicted second-harmonic peak at
$\omega/2\pi\simeq 2\nu_{0}$ arises both in SNN simulations and in the SFP
integration. Despite such footprint of nonlinear dynamics, SNN and SFP keep
displaying a remarkable overlap not only in $P_{\nu}$ but also in the
stochastic dynamics of $\nu_{N}(t)$ shown in Fig. 4(b) where coherent
oscillations become more and more apparent as $N$ decreases. The ongoing
finite-size fluctuations of $\nu$ in this case continuously stimulate the
oscillating relaxation of the stable-focus which is relatively slow due to the
nearby critical point. The $N$-dependent coherence of the fluctuation-driven
oscillation is even more apparent in the increase of power of the second-
harmonic peak in Fig. 4(a)-bottom as the size of the network is reduced.

The same remarkable match between SFP and SNN is shown in Fig. 4(c)-top beyond
the Hopf bifurcation ($KJ=12$ mV) when a stable limit cycle is expected at
$N\to\infty$. Global oscillations in this case are strongly nonlinear giving
rise to several high-order harmonic peaks and finite-size noise contribute to
dephase them limiting their coherence in time [Fig. 4(d)]. According to this,
the resonant peaks in $P_{\nu}(\omega)$ have power lowering with decreasing
$N$ [Fig. 4(c)-bottom].

## Conclusions. —

The agreement found between SFP and SNN in a synchronization phase transition
is an encouraging evidence that the stochastic population dynamics we derived
in Eqs. (5), (11) and (4) has the potential to faithfully describe out-of-
equilibrium networks composed of a finite number $N$ of spiking neurons.
Indeed, it is important to stress that the theoretical framework we introduce
is not perturbative and takes into account the point-like (spike-based) nature
of the cell-to-cell interactions. Furthermore, the correlation structure of
$\eta(t)$ here results from the spiking statistics of isolated neurons,
implying that finite-size noise is not $J$\- dependent. Not only, its
Markovian embedding allows in principle to derive single-neuron ISI statistics
also under nonstationary conditions, i.e., when renewal hypothesis is broken.

Alternative approaches dealing with finite-size networks of spiking neurons
are those relying on the refractory density method (RDM) [36]. Neurons in this
case are quasi-stationary renewal processes and the population dynamics is
fully described by the probability density of a single-cell to be into a
refractory state [15, 16]. In the RDM, the integro-differential equation
governing the population activity at the mesoscopic scale can incorporate
fluctuations to describe a finite number of neurons [37, 35]. However, the
renewal hypothesis underlying RDM may in principle limit its applicability by
making our theoretical framework preferable in dealing with out-of-equilibrium
conditions.

In conclusion, we remark that the theoretical framework we developed applies
to a wide class of spiking neuron models, and the same approach can be
extended to other physical systems with a finite size. Indeed, the inclusion
of fluctuations in population density theories has the potential to further
advance our understanding of noise-driven out-of-equilibrium dynamics [9, 12].

Besides, numerical integrations of SFP in the examples shown computationally
outperform detailed simulations of SNNs. In principle, this paves the way to
implement large-scale heterogeneous networks of neuronal population making
affordable detailed in silico experiments on macroscopic brain regions [34].

## Acknowledgments. —

Work in part funded by EU H2020 Research and Innovation Programme, Grant
945539 (HBP SGA3) to MM.

## References

  * [1] L. F. Abbott and C. van Vreeswijk.  Asynchronous states in networks of pulse-coupled oscillators.  Phys. Rev. E, 48(2):1483–1490, 1993.
  * [2] D. J. Amit and N. Brunel.  Model of global spontaneous activity and local structured activity during delay periods in the cerebral cortex.  Cereb. Cortex, 7(3):237–52, 1997.
  * [3] M. Augustin, J. Ladenbauer, F. Baumann, and K. Obermayer.  Low-dimensional spike rate models derived from networks of adaptive integrate-and-fire neurons: Comparison and implementation.  PLoS Comput. Biol., 13(6):e1005545, 2017.
  * [4] N. Brunel and V. Hakim.  Fast global oscillations in networks of integrate-and-fire neurons with low firing rates.  Neural Comput., 11(7):1621–71, 1999.
  * [5] P.-H. Chavanis.  Hamiltonian and Brownian systems with long-range interactions: V. Stochastic kinetic equations and theory of fluctuations.  Physica A, 387(23):5716–5740, 2008.
  * [6] P.-H. Chavanis.  Generalized stochastic Fokker-Planck equations.  Entropy, 17(5):3205–3252, 2015.
  * [7] C.-C. Chen, A. B. Herhold, C. S. Johnson, and A. P. Alivisatos.  Size dependence of structural metastability in semiconductor nanocrystals.  Science, 276(5311):398–401, 1997.
  * [8] D. R. Cox and H. D. Miller.  The theory of stochastic processes, volume 134.  CRC press, 1977.
  * [9] S. K. Das, S. Roy, S. Majumder, and S. Ahmad.  Finite-size effects in dynamics: Critical vs. coarsening phenomena.  Europhys. Lett., 97(6):66006, 2012.
  * [10] D. S. Dean.  Langevin equation for the density of a system of interacting Langevin processes.  J. Phys. A, 29(24), 1996.
  * [11] J. DeFelipe, L. Alonso-Nanclares, and J. I. Arellano.  Microstructure of the neocortex: comparative aspects.  J. Neurocytol., 31(3-5):299–316, 2002.
  * [12] M. A. Durán-Olivencia, P. Yatsyshin, B. D. Goddard, and S. Kalliadasis.  General framework for fluctuating dynamic density functional theory.  New J. Phys., 19(12), 2017.
  * [13] T. Fardet, S. B. Vennemo, J. Mitchell, H. Mørk, S. Graber, J. Hahne, and et al.  NEST 2.20.0.  Zenodo, 2020.  https://doi.org/10.5281/zenodo.3605514.
  * [14] S. Fusi and M. Mattia.  Collective behavior of networks with linear (VLSI) integrate-and-fire neurons.  Neural Comput., 11(3):633–52, 1999.
  * [15] W. Gerstner.  Time structure of the activity in neural network models.  Phys. Rev. E, 51(1):738–758, 1995.
  * [16] W. Gerstner.  Population dynamics of spiking neurons: fast transients, asynchronous states, and locking.  Neural Comput., 12(1):43–89, 2000.
  * [17] W. Gerstner, W. M. Kistler, R. Naud, and L. Paninski.  Neuronal dynamics: From single neurons to networks and models of cognition.  Cambridge University Press, 2014.
  * [18] B. Grigelionis.  On the convergence of sums of random step processes to a Poisson process.  Theory Probab., 8(2):177–182, 1963.
  * [19] S. Gupta, A. Campa, and S. Ruffo.  Kuramoto model of synchronization: equilibrium and nonequilibrium aspects.  J. Stat. Mech., 2014(8):R08001, 2014.
  * [20] I. S. Hwang, S. H. Chang, C. K. Fang, L. J. Chen, and T. T. Tsong.  Observation of finite-size effects on a structural phase transition of 2D nanoislands.  Phys. Rev. Lett., 93(10):1–4, 2004.
  * [21] K. Kawasaki.  Stochastic model of slow dynamics in supercooled liquids and dense colloidal suspensions.  Physica A, 208(1):35–64, 1994.
  * [22] B. W. Knight.  Dynamics of encoding in a population of neurons.  J. Gen. Physiol., 59(6):734–66, 1972.
  * [23] B. W. Knight, D. Manin, and L. Sirovich.  Dynamical models of interacting neuron populations in visual cortex.  In E.C. Gerf, editor, Symposium on Robotics and Cybernetics: Computational Engineering in Systems Applications, pages 1–5, Cite Scientifique, Lille, France, 1996. Cite Scientifique.
  * [24] B. Lindner.  Superposition of many independent spike trains is generally not a Poisson process.  Phys. Rev. E, 73(2):1–4, 2006.
  * [25] H. Markram, E. Muller, S. Ramaswamy, M. W. Reimann, M. Abdellah, and et al.  Reconstruction and simulation of neocortical microcircuitry.  Cell, 163(2):456–492, 2015.
  * [26] M. Mattia, M. Biggio, A. Galluzzi, and M. Storace.  Dimensional reduction in networks of non-Markovian spiking neurons: Equivalence of synaptic filtering and heterogeneous propagation delays.  PLoS Comput. Biol., 15(10):e1007404, 2019.
  * [27] M. Mattia and P. Del Giudice.  Population dynamics of interacting spiking neurons.  Phys. Rev. E, 66(5 Pt 1):051917, 2002.
  * [28] M. Mattia and P. Del Giudice.  Finite-size dynamics of inhibitory and excitatory interacting spiking neurons.  Phys. Rev. E, 70(5 Pt 1):052903, 2004.
  * [29] M. Mattia and G. V. Vinci.  Low dimensional dynamics of spiking neuron networks.  Zenodo, 2021.  https://doi.org/10.5281/zenodo.5518215.
  * [30] V. B. Mountcastle.  The columnar organization of the neocortex.  Brain, 120(Pt 4):701–22, 1997.
  * [31] M. Pascual and S. A. Levin.  From individuals to population densities: Searching for the intermediate scale of nontrivial determinism.  Ecology, 80(7):2225–2236, 1999.
  * [32] C. Poo and J. S. Isaacson.  Odor representations in olfactory cortex: ”sparse” coding, global inhibition, and oscillations.  Neuron, 62(6):850–861, 2009.
  * [33] M. Rigotti, O. Barak, M. R. Warden, X.-J. Wang, N. D. Daw, E. K. Miller, and S. Fusi.  The importance of mixed selectivity in complex cognitive tasks.  Nature, 497(7451):585–90, 2013.
  * [34] P. Sanz Leon, S. A. Knock, M. M. Woodman, L. Domide, J. Mersmann, A. R. Mcintosh, and V. Jirsa.  The virtual brain: A simulator of primate brain network dynamics.  Front. Neuroinform., 7:10, 2013.
  * [35] V. Schmutz, E. Löcherbach, and T. Schwalger.  On a finite-size neuronal population equation.  arXiv, 2021.  2106.14721.
  * [36] T. Schwalger and A. V. Chizhov.  Mind the last spike — firing rate models for mesoscopic populations of spiking neurons.  Curr. Opin. Neurobiol., 58:155–166, 2019.
  * [37] T. Schwalger, M. Deger, and W. Gerstner.  Towards a theory of cortical columns: From spiking neurons to interacting neural populations of finite size.  PLoS Comput. Biol., 13(4):e1005507, 2017.
  * [38] H. C. Tuckwell.  Introduction to theoretical neurobiology: volume 2, nonlinear and stochastic theories, volume 8.  Cambridge University Press, 1988.
  * [39] N. G. van Kampen.  Stochastic Processes in Physics and Chemistry.  North Holland, 3rd edition, 2007.
  * [40] S. Vellmer and B. Lindner.  Theory of spike-train power spectra for multidimensional integrate-and-fire neurons.  Phys. Rev. Res., 1(2):023024, 2019.
  * [41] G. V. Vinci and M. Mattia.  A ‘Rosetta stone’ for the population dynamics of spiking neuron networks.  Zenodo, 2021.  https://doi.org/10.5281/zenodo.5519083.
  * [42] B. O. Watson, D. Levenstein, J. P. Greene, J. N. Gelinas, and G. Buzsáki.  Network homeostasis and state dynamics of neocortical sleep.  Neuron, 90(4):839–52, 2016.


# Identifying the Gamma-ray Emission of the Nearby Galaxy M83

Yi Xing  Key Laboratory for Research in Galaxies and Cosmology, Shanghai
Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road,
Shanghai 200030, China<EMAIL_ADDRESS><EMAIL_ADDRESS>Zhongxiang
Wang  Department of Astronomy, School of Physics and Astronomy, Yunnan
University, Kunming 650091, China  Key Laboratory for Research in Galaxies and
Cosmology, Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80
Nandan Road, Shanghai 200030, China<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We report on the detection of a $\gamma$-ray source at the position of the
nearby star-forming galaxy (SFG) M83, which is found from our analysis of 14
years of the data obtained with the Large Area Telescope (LAT) on-board Fermi
Gamma-ray Space Telescope (Fermi). The source is weakly detected, with a
significance of $\sim 5\sigma$, and its emission can be described with an
exponentially cutoff power law. At a distance of 4.61 Mpc, the source’s
$\gamma$-ray luminosity is $\sim 1.4\times 10^{39}$ erg s-1, roughly along the
correlation line between the $\gamma$-ray and IR luminosities determined for
nearby SFGs. Because of the weak detection, the source spectrum can not be
used for checking its similarity with those of other SFGs. Given the
positional matches and the empirical expectation for $\gamma$-ray emission
from M83 due to the galaxy’s star-forming activity, we conclude that the
$\gamma$-ray source is the likely counterpart to m83. The detection thus adds
another member to the group of approximately a dozen SFGs, whose $\gamma$-ray
emissions mostly have a cosmic-ray origin.

Gamma-ray sources (633); Starburst galaxies (1570)

##  1 Introduction

Among more than 6000 $\gamma$-ray sources detected with the Large Area
Telescope (LAT) on board the Fermi Gamma-ray Space Telescope (Fermi) in all
sky, the dominant class is active galactic nuclei (AGN; Abdollahi et al.
2022), whose high-energy emission is mostly radiated from their jets. Non-
active galaxies thus do not show such emission. However there are
approximately a dozen of the galaxies, either within the local group or
nearby, have been detected at $\gamma$-rays (Ajello et al., 2020; Xi et al.,
2020). While there are complications in the production of the $\gamma$-ray
emissions observed from these galaxies, for example AGN possibly hiding in
some of them (Peng et al., 2019) and the $\gamma$-ray emission of the local-
group galaxy M31 being considered consisting of different components (Li et
al., 2016; Pshirkov et al., 2016; Ackermann et al., 2017; Karwin et al., 2019;
Zimmer et al., 2022; Xing et al., 2023), $\gamma$-ray luminosities of most of
them well correlate with infrared (IR) or radio 1.4 GHz luminosities (Abdo et
al., 2010b; Ackermann et al., 2012; Ajello et al., 2020; Xi et al., 2020).
This correlation is considered as an indicator for the cosmic ray (CR) origin
of the gamma-ray emissions. Supernova remnants (SNRs) produce CRs as their
shock fronts serve as particle accelerators (e.g., Bykov et al. 2018), and the
accelerated particles emit at radio frequencies through the synchrotron
process and at high energies through the proton-proton collisions and/or
leptonic processes (i.e., bremsstrahlung or inverse Compton scattering; e.g.,
Dermer 1986). On the other hand, SNRs are the results of massive stars
($M\gtrsim 8\,M_{\sun}$) evolving to the end of their lives on relatively
short, $\sim 10^{7}$ yr timescales, and their densities are thus closely
related to star-formation of a galaxy. Given that the IR luminosities are an
indicator of star-formation rates of a galaxy, a correlation between them and
corresponding $\gamma$-ray luminosities is naturally expected for star-forming
galaxies (see, e.g., Domingo-Santamaría & Torres 2005; Lacki et al. 2010;
Ackermann et al. 2012, and references therein).

Along this expected correlation, efforts have been made to detect nearby star-
forming galaxies at $\gamma$-rays. Thus far, approximately a dozen of them
have been reported with the detection (see, e.g., Ajello et al. 2020; Xi et
al. 2020, and references therein). Here we report on the likely detection of
another one, the nearby galaxy M83 (also known as NGC 5236 or the Southern
Pinwheel).

Figure 1: TS maps of a $3\arcdeg\times 3\arcdeg$ region centered at M83 in the
energy ranges of 0.1–500 GeV (left panel) and 0.5–500 GeV (middle panel). The
only 4FGL-DR3 catalog source in the region, a blazar candidate (marked by
green pluses), is removed in the maps. The green crosses mark the position of
M83, which is within the 2$\sigma$ error circle (marked by the green circles)
determined for the residual excess emission at the location. There is also
weak excess emission (whose 2$\sigma$ error circle is marked by green dashed
circles) NE to the center. Right panel: the same as the middle panel, but the
NE source is removed as a point source. For each TS map, the image scale is
0$\fdg$05 pixel-1, and the color bar indicates the TS value range.

M83 is often referred to as a grand-design spiral galaxy. It is nearly face-on
to us ($i\sim 25\arcdeg$, Sofue et al. 1999), at a distance of 4.61 Mpc (Saha
et al., 2006). There have been extensively studies of the galaxy over the
whole wavelength range. Its star formation is relatively active, having a
total star-formation rate of 5 $M_{\sun}$ yr-1 (Kennicutt, 1998). Thus it was
listed in the star-forming galaxy sample selected by Ackermann et al. (2012),
to be searched for CR-induced $\gamma$-ray emission. A flux upper limit was
reported in Ackermann et al. (2012), while only 3 years of the Fermi-LAT data
were used.

Now with 14 years of the data having been collected with Fermi-LAT, we
conducted a search for M83’s $\gamma$-ray emission. We found a likely
counterpart and report the results. Below we describe the details of our
analysis and provide the corresponding results in Section 2. The results are
discussed in Section 3.

##  2 Analysis and Results

###  2.1 Fermi-LAT Data and Source Model

We selected 0.1–500 GeV LAT events from the updated Fermi Pass 8 database in a
region of interest (RoI) that has a size of 20∘ $\times$ 20∘ and the center at
the central position of M83. Since the galaxy appears to have a size $\sim
12\arcmin$ in the sky, which is not resolvable in the LAT data due to its
large point spread function (PSF), we treat M83 as a point source through this
paper. The time period of the LAT data was from 2008-08-04 15:43:39 (UTC) to
2022-09-26 23:16:35 (UTC), slightly more than 14 yrs. The CLEAN event class
was used. We included the events with zenith angles less than 90 deg and
excluded the events with quality flags of ‘bad’. Both these are recommended by
the LAT team111http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/.

We constructed the source model by including all sources within 20 deg from
M83. The positions and the spectral parameters of these sources are provided
in the Fermi LAT 12-year source catalog (4FGL-DR3; Abdollahi et al. 2022). We
set the spectral parameters of the sources within 5 deg from M83 free, and
froze the other parameters at their catalog values. The spectral model
gll_iem_v07.fit was used for the Galactic diffuse emission, and the spectral
file iso_P8R3_CLEAN_V3_v1.txt for the extragalactic diffuse emission. The
normalizations of these two diffuse components were set free in the following
analyses.

Table 1: Likelihood analysis results with the PL, PLEC, and LP models Model | Best-fit parameters | $\log(L)$ | TS
---|---|---|---
PL | $\Gamma=2.0\pm 0.2$ | 265522.5 | 14
PLEC | $\Gamma=0.00\pm 0.03$ | 265526.9 | 22
| $E_{c}=2.8\pm 0.8$ |  |
LP | $\alpha=0.00\pm 0.06$ | 265526.1 | 20
| $\beta=0.67\pm 0.14$ |  |

###  2.2 Likelihood Analysis

Figure 2: Same as the middle panel of Figure 1, but the blazar candidate is
kept in the maps. The left panel TS map is calculated from the whole LAT data,
and the right panel one from the time period of MJD 55042–59000 (see Section
2.3 and Figure 3).

We performed the standard binned likelihood analysis to the whole data in
0.1–500 GeV and updated the parameter values for the sources within 5 deg from
M83. With the obtained results, we calculated a 0.1–500 GeV residual Test
Statistic (TS) map of a $\mathrm{3^{o}\times 3^{o}}$ region centered at M83.
The field is rather clean, with only one catalog source (4FGL J1335.3$-$2949;
see Section 2.3) in the region. This source was included in the source model
and removed in the TS map, which is shown in the left panel of Figure 1. As
can be seen, there is weak excess emission at the position of M83. The maximum
TS value is approximately 14, corresponding to a $\sim$3.7$\sigma$ detection
significance. Because in the analysis below (Section 2.4 and Table 2) we have
found that the low energy range for significant detection is $\gtrsim$0.5 GeV,
we also calculated the 0.5–500 GeV residual TS map (shown in the middle panel
of Figure 1). The maximum TS value of the excess emission at M83 is improved
to be $\sim$25, now at a $\sim$5$\sigma$ detection significance. We ran
gtfindsrc in the Fermitools to the 0.5–500 GeV data to determine the position,
and obtained R.A.=204$\fdg$23, Decl.=$-$29$\fdg$83 (equinox J2000.0), with a
1$\sigma$ nominal uncertainty of 0$\fdg$04\. M83 is 0$\fdg$04 away from this
position and thus within the 1$\sigma$ error circle.

We then added this new source, the possible counterpart to M83, in the source
model as a point source and repeated the likelihood analysis in 0.1–500 GeV.
Given the low TS value of the source, we considered three models to fit its
emission. One is a simple power law (PL), $dN/dE=N_{0}E^{-\Gamma}$, where
$\Gamma$ is the photon index, and the other two are a PL with an exponential
cutoff (PLEC), $dN/dE=N_{0}E^{-\Gamma}\exp(-E/E_{c})$, where $E_{c}$ is the
cutoff energy, and a Log-Parabola (LP) function,
$dN/dE=N_{0}(E/E_{b})^{-[\alpha+\beta\log(E/E_{b})]}$, where $E_{b}$ is a
scale parameter and was fixed at 1 GeV in our analysis. The results obtained
with the three models are given in Table 1. The PLEC model provided the
largest TS value, $\simeq$22\. By comparing the log-likelihood values, that is
$\sqrt{-2log(L_{i}/L_{j})}$, where $L_{i/j}$ are the maximum likelihood values
from model $i$ and $j$, the PLEC and LP models are found to be more favored
than the PL at $\sim$3.0$\sigma$ and $\sim$2.7$\sigma$ significances
respectively. Although $\gamma$-ray emissions of most star-forming galaxies
can be well described with a PL model (or a LP model; e.g., Ajello et al.
2020), below we adopted the PLEC model because of the largest TS value
resulting from it. The corresponding 0.1–500 GeV photon flux for the source
was $F_{0.1-500}\simeq 1.2\pm 0.4\times 10^{-10}$ photon cm-2 s-1.

Figure 3: 60-day binned light curves (upper) and TS curves (bottom) of the
blazar candidate (red) and the source at M83 (gray) in 0.1–500 GeV. For the
latter, 2-yr binned light curve and TS curve (black data points) are also
shown. Fluxes with TS$\geq$4 are kept in the light curves. Two dotted lines
mark the time period of MJD 55042–59000, during which no obvious flares are
seen from the blazar candidate.

###  2.3 Analysis for Nearby Sources

Since the emission at the position of M83 was weak, we conducted extra-checks
to ensure the detection. First as shown in the TS maps (Figure 1), another
excess emission is present, which is north-east (NE) to M83. Although it does
not appear like a point source, we ran gtfindsrc to the 0.5–500 GeV data, and
obtained a position of R.A.=205$\fdg$3, Decl.=$-$28$\fdg$9 (equinox J2000.0),
with a 1$\sigma$ nominal uncertainty of 0$\fdg$2\. This position is 1.3 deg
away from M83, nearly outside of the 68% containment angle of the LAT PSF in
$>$0.5 GeV
band222https://www.slac.stanford.edu/exp/glast/groups/canda/lat_Performance.htm.
We added this source in the source model and performed the likelihood
analysis. The source could be totally removed in the TS map (shown in the
right panel of Figure 1), and the results for the emission at M83 were nearly
the same as above. These suggest that the NE excess emission did not affect
our analysis results for the source at M83.

Second, the catalog source, 4FGL J1335.3$-$2949, is located very close to M83
(see Figure 1). It has an angular separation of 0$\fdg$35 from M83 and its
1$\sigma$ positional uncertainty is 0$\fdg$02, given in 4FGL-DR3 (Abdollahi et
al., 2022). Thus this source and M83’s source are outside of the error circle
of each other. Because it is relatively bright, it could contaminate the
detection of the source at M83. We calculated the 0.5–500 GeV TS map with the
sources (including the NE one) in the $\mathrm{3^{o}\times 3^{o}}$ region
kept. 4FGL J1335.3$-$2949 is clearly seen as the brightest one (the TS value
was $\sim$206) in the field (left panel of Figure 2).

This source is a blazar candidate but did not show significant $\gamma$-ray
variations in 12 years of the Fermi-LAT data (Abdollahi et al., 2022). We
extracted its 0.1–500 GeV light curve by setting 60-day time bins and
performing the likelihood analysis to each time-bin data. In the extraction,
only the normalization parameters of the sources within 5 deg from M83 were
set free. As a check, we also extracted a light curve of the source at M83
with the same setup. The resulting light curves and TS curves for the two
sources are shown in Figure 3, for which only flux data points with TS$\geq$4
were kept in the light curves. As can be seen, this blazar candidate had a
high TS value in the beginning of the data at $\sim$MJD 55000 and likely
showed a flaring event after MJD 59000 (note that the latter part was not
covered in the LAT 12-yr source catalog). The source at M83 did not show any
obvious variations.

To reduce any contamination possibly caused by the variations of 4FGL
J1335.3$-$2949, we selected the LAT data during the time period of MJD
55042–59000 (marked as dotted lines in Figure 3). The 0.5–500 GeV TS map of
this time period was calculated and is shown in the right panel of Figure 2.
The sources in the TS-map region were kept. The TS value for the blazar
candidate is reduced to $\sim$100, and the source at M83, appearing as an
extension of the former to the east direction, is revealed. Thus the flaring
activity of the blazar candidate could affect our analysis and weaken the
visibility of the source at M83, but when the flares were filtered out, the
detection of this source was proved to be true.

Figure 4: $\gamma$-ray spectra obtained for M83’s source (black circles) and
the blazar candidate (red diamonds), for which the model fits from the
likelihood analysis are also shown for comparison, where the black dashed and
solid curves are the PLEC and PL model fits (cf., Section 2.2), respectively,
for the former source and the red line is the PL model fit for the latter one.
The model fit to the normalized spectra of 9 star-forming galaxies determined
by Ajello et al. (2020) is shown as the blue curve, while its 1$\sigma$
uncertainty range is marked by the gray region (here this model-fit curve is
simply aligned with the first flux measurement of M83’s source).

###  2.4 Spectral Analysis

We extracted the $\gamma$-ray spectrum of the source at M83 by performing
maximum likelihood analysis to the LAT data in 10 evenly divided energy bands
in logarithm from 0.1–500 GeV. In the extraction, the spectral normalizations
of the sources within 5 deg from M83 were set as free parameters, while all
the other parameters of the sources were fixed at the values obtained from the
above maximum likelihood analysis. The emission was set to be a PL with
$\Gamma$ fixed to 2. For the result, we kept only spectral data points when
TS$\geq 4$ ($\geq$2$\sigma$ significance) and derived 95% flux upper limits
otherwise, where for the latter a Bayesian approch implemented in the Python
tool IntegralUpperLimit (provided in the Fermi Science tools) was used. The
obtained spectrum is plotted as black points in Figure 4, and the flux and TS
values of the spectral data points are provided in Table 2. In the energy band
of 0.2–0.5 GeV, the TS value is 5, but the flux has a large uncertainty,
0.42$\pm$0.62$\times$10-12 erg cm-2 s-1, likely caused by the contamination of
the nearby blazar candidate (because of large containment angles of the LAT
PSF at the low energies; see also Figure 4). Thus for this data point, we
report an upper limit instead. In addition, the spectrum of the nearby blazar
candidate was obtained with the same setup, and is also shown in Figure 4 for
comparison.

###  2.5 Variability Analysis

We checked the source at M83 for any long-term variability in 0.1–500 GeV by
calculating the variability index TSvar (Nolan et al., 2012). We set 87 time
bins with each bin consisting of 60-day data and derived the fluxes or flux
upper limits for the source (i.e., shown in Figure 3). If the emission is
constant, TSvar would be distributed as a $\chi^{2}$ distribution with 86
degrees of freedom. A variable source would be identified if TSvar is larger
than 119.4 (at a 99% confidence level). The computed TSvar for the source is
74.8, lower than the threshold value. Since this source is faint, we also
constructed its 2-yr binned light curve (see Figure 3) and checked for
variability. For 7 time bins (i.e., 6 degrees of freedom), TS${}_{var}>16.8$
is required for a variable source. We obtained TS${}_{var}\simeq 11.3$. Thus
there were no significant long-term variations found for this source.

Table 2: Flux Measurements

Band | $G_{M83}/10^{-12}$ | TS
---|---|---
(GeV) | (erg cm-2 s-1) |
0.15 (0.1–0.2) | 0.61 | 0
0.36 (0.2–0.5) | 0.74 | 5
0.84 (0.5–1.3) | 0.34$\pm$0.19 | 7
1.97 (1.3–3.0) | 0.20$\pm$0.10 | 6
4.62 (3.0–7.1) | 0.17$\pm$0.09 | 5
10.83 (7.1–16.6) | 0.31$\pm$0.14 | 13
25.37 (16.6–38.8) | 0.31 | 0
59.46 (38.8–91.0) | 0.89 | 0
139.36 (91.0–213.3) | 1.58 | 0
326.60 (213.3–500.0) | 3.69 | 0

Note: Fluxes without uncertainties are the 95$\%$ upper limits.

##  3 Discussion

By analyzing the Fermi-LAT data for the M83 region, we have found a faint
$\gamma$-ray source at the position of M83. The source’s emission is
preferably described with a curved-function model (a PLEC or a LP). Although
the detection significance for the source is low, only $\sim 4.7\sigma$, and a
nearby blazar candidate complicates the detection (particularly at the
$\lesssim$0.5 GeV low energies), our detailed analysis has ensured the
existence of the source. Given the high positional coincidence between the
$\gamma$-ray source and M83 and the expectation for $\gamma$-ray emission from
this star-forming galaxy, we conclude that this $\gamma$-ray source is the
likely counterpart to M83.

At a distance of 4.61 Mpc (Saha et al., 2006), the $\gamma$-ray luminosity (in
0.1–500 GeV) of the source is $\sim 1.4\pm 0.5\times 10^{39}$ erg s-1, higher
than that of the Milky way (Ackermann et al., 2012) and lower than that of NGC
253 (a spiral galaxy with a central starburst region; Abdo et al. 2010a;
Ajello et al. 2020). The 2–1000 $\mu$m IR luminosity of M83 is approximately
$8.7\times 10^{43}$ erg s-1 (Ackermann et al. 2012; note that a source
distance of 3.7 Mpc was used in Ackermann et al. 2012). If we use the
parameters obtained in Ajello et al. (2020) for the $\gamma$-ray–IR luminosity
correlation, the predicted $\gamma$-ray luminosity for M83 would be $\sim
7.9\times 10^{39}$ erg s-1. Further considering a dispersion of 0.3
($\gamma$-ray–luminosity residuals with respect to the correlation line in log
space; Ackermann et al. 2012; Ajello et al. 2020), the luminosity range would
be $\sim 4.0$–$16\times 10^{39}$ erg s-1. Thus the observed luminosity is
slightly below the correlation line but can be considered to be consistent
with it, given significant uncertainties such as on the distance and IR
properties (Ajello et al. 2020 and references therein).

Because of the weak detection of M83, its $\gamma$-ray spectrum only contains
4 data points (Figure 4). The model fits we obtained in Section 2.2, either
the PLEC or the LP, appear to be highly curved (for example, $\beta\approx
0.67$ in the LP model), which are different from those of the other star-
forming galaxies. The latter are mostly described with a PL without
significant curvature required (e.g., Ajello et al. 2020). The difference
could be due to the weak detection causing the emission property not well
determined, and remains to be resolved when more LAT data for M83 are
collected. Ajello et al. (2020) have normalized the $\gamma$-ray spectra of 9
star-forming galaxies by simply scaling the spectra to a common value at 1
GeV, and obtained a best-fit model that is in the form of a smoothly broken
power law. In Figure 4, we show this best fit by aligning it with the first
flux measurement of M83’s spectrum, since the energy range of the data point
is 0.5–1.3 GeV, approximately at 1 GeV. As can be seen, the spectrum of M83 is
approximately consistent with the best fit. Given this and the $\gamma$-ray
luminosity indicated above being approximately in the right range, the
$\gamma$-ray emission of M83 likely has the same origin as the other nearby
galaxies, arising from CRs and related to the star-formation activity
(Ackermann et al., 2012). Thus M83 is likely another member of this gourp of
the $\gamma$-ray–emitting, star-forming galaxies. Hopefully with more LAT data
collected in the future, both the detection significance and the quality of
the $\gamma$-ray spectrum will be improved, helping provide more information
for the high-energy properties of M83.

This research is supported by the Original Innovation Program of the Chinese
Academy of Sciences (E085021002), the Basic Research Program of Yunnan
Province No. 202201AS070005, and the National Natural Science Foundation of
China No. 12273033.

## References

  * Abdo et al. (2010a) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, ApJ, 709, L152
  * Abdo et al. (2010b) —. 2010b, A&A, 523, L2
  * Abdollahi et al. (2022) Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53
  * Ackermann et al. (2012) Ackermann, M., Ajello, M., Allafort, A., et al. 2012, ApJ, 755, 164
  * Ackermann et al. (2017) Ackermann, M., Ajello, M., Albert, A., et al. 2017, ApJ, 836, 208
  * Ajello et al. (2020) Ajello, M., Di Mauro, M., Paliya, V. S., & Garrappa, S. 2020, ApJ, 894, 88
  * Bykov et al. (2018) Bykov, A. M., Ellison, D. C., Marcowith, A., & Osipov, S. M. 2018, Space Sci. Rev., 214, 41
  * Dermer (1986) Dermer, C. D. 1986, A&A, 157, 223
  * Domingo-Santamaría & Torres (2005) Domingo-Santamaría, E., & Torres, D. F. 2005, A&A, 444, 403
  * Karwin et al. (2019) Karwin, C. M., Murgia, S., Campbell, S., & Moskalenko, I. V. 2019, ApJ, 880, 95
  * Kennicutt (1998) Kennicutt, Robert C., J. 1998, ApJ, 498, 541
  * Lacki et al. (2010) Lacki, B. C., Thompson, T. A., & Quataert, E. 2010, ApJ, 717, 1
  * Li et al. (2016) Li, Z., Huang, X., Yuan, Q., & Xu, Y. 2016, J. Cosmology Astropart. Phys, 2016, 028
  * Nolan et al. (2012) Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, ApJS, 199, 31
  * Peng et al. (2019) Peng, F.-K., Zhang, H.-M., Wang, X.-Y., Wang, J.-F., & Zhi, Q.-J. 2019, ApJ, 884, 91
  * Pshirkov et al. (2016) Pshirkov, M. S., Vasiliev, V. V., & Postnov, K. A. 2016, MNRAS, 459, L76
  * Saha et al. (2006) Saha, A., Thim, F., Tammann, G. A., Reindl, B., & Sandage, A. 2006, ApJS, 165, 108
  * Sofue et al. (1999) Sofue, Y., Tutui, Y., Honma, M., et al. 1999, ApJ, 523, 136
  * Xi et al. (2020) Xi, S.-Q., Zhang, H.-M., Liu, R.-Y., & Wang, X.-Y. 2020, ApJ, 901, 158
  * Xing et al. (2023) Xing, Y., Wang, Z., Zheng, D., & Li, J. 2023, ApJ, 945, L22
  * Zimmer et al. (2022) Zimmer, F., Macias, O., Ando, S., Crocker, R. M., & Horiuchi, S. 2022, MNRAS, 516, 4469