Besançon Galactic Model Overview

• The Besançon Galactic Model is a comprehensive population-synthesis and forward-modeling framework that simulates the Milky Way using explicit IMF, SFR, and stellar evolution tracks.
• It integrates star counts, photometric observables, kinematics, and extinction models to generate mock catalogs that validate Galactic-evolution scenarios.
• Advancements in self-consistent dynamical modeling and efficient binary-star treatment enable precise inference of the thin disc, bulge, and halo structures.

Searching arXiv for recent and foundational Besançon Galactic Model papers to ground the article. The Besançon Galactic Model (BGM) is a population-synthesis model of the Milky Way designed not only to predict star counts, but to test Galactic-evolution scenarios by specifying the stellar birth process through an initial mass function (IMF) and a star formation rate (SFR), then evolving stars along tracks and embedding them in a Galaxy whose structural parameters are constrained observationally and dynamically. In its modern forms, BGM links stellar populations, extinction, photometric observables, kinematics, and gravitational structure in a forward-modeling framework that has been used for all-sky star counts, bulge and halo structure, microlensing, Gaia-era dynamical inference, and survey mock-catalog generation (Czekaj et al., 2014).

1. Concept, scope, and historical development

BGM represents the Galaxy as a superposition of major stellar components, typically including the thin disc, thick disc, halo or stellar halo, and bulge or bar. A defining feature, emphasized especially in the renewed thin-disc work, is that BGM is not only a descriptive star-count code but a model for testing Galactic-evolution hypotheses under explicit assumptions on the IMF, SFR, stellar evolution, metallicity evolution, extinction, and kinematics. Earlier public implementations, notably the 2003 version, already combined population synthesis with dynamical constraints in the Solar neighbourhood, but they used a more rigid thin-disc generation scheme inherited from Hipparcos-era calibrations, lacked consistent binary generation, and showed systematic failures in bright-star comparisons, including overproduction of nearby A–F dwarfs and global discrepancies in Tycho-2 color distributions (Czekaj et al., 2014).

The model’s later development followed two main directions. One direction renewed the stellar-population machinery, especially for the thin disc, by turning the IMF, SFR, and stellar evolutionary tracks into explicit free ingredients of the simulation and by extending the model toward survey-scale inference. The other direction deepened the dynamical treatment by replacing empirical disc kinematics with explicit distribution functions depending on three integrals of motion in a stationary axisymmetric potential, yielding a self-consistent density–potential–kinematics loop suitable for Gaia-era analyses (Bienaymé et al., 2018).

2. Stellar population synthesis and forward modeling

At the population-synthesis level, BGM generates stars from an IMF and an SFR, evolves them along stellar tracks, and transforms intrinsic quantities into observables through atmosphere libraries and extinction models. One compact formulation used in a later application writes the star-count basis as

$N_* = \int \Phi(m)\, dm \int \Psi(t)\, dt,$

where $\Phi(m)$ is the IMF and $\Psi(t)$ the SFR. In the renewed thin-disc formalism, the key local quantity for a given age subcomponent $i$ and simulated volume element $dV$ is the mass reservoir,

$mass~reservoir = dV \times \rho(x,y,z,i),$

which is then partitioned across IMF intervals and converted into stars with assigned mass, age, metallicity, and evolutionary state (Mason et al., 25 Nov 2025).

For the thin disc, the IMF is implemented as a piecewise power law,

$\phi(m) = m^{-(1+x)},$

with continuity coefficients across mass breaks. Stars are generated by drawing mass from the IMF, age uniformly within the chosen age-bin interval, and metallicity from the star’s own age through an adopted age–metallicity relation rather than from the mean metallicity of the age bin. The code then interpolates in stellar evolutionary tracks to place the star in the HR diagram; if no living solution exists, the object is sent to the remnant pool. A notable numerical correction addresses very small local mass reservoirs: nearby cells are temporarily enlarged by a factor of 50 to avoid a biased realized IMF at high mass, and only the appropriate fraction of generated stars is then retained (Czekaj et al., 2014).

Binary-star treatment became a major conceptual extension of the renewed BGM. Every newly created living star is probabilistically assigned as single or as the primary of a binary system, with empirically based probabilities and secondary-star distributions. Secondaries inherit the primary’s age and metallicity, are placed on the same tracks, and are subtracted from the same mass reservoir, so multiplicity is generated within a fixed mass budget. The model deliberately does not assume that secondaries follow the same IMF as single stars; rather, it uses observed binary statistics, accepting a small distortion of the global IMF of all stars. In the same renewal, white dwarfs remained an explicit limitation, since they were still produced following the old model instead of the new self-consistent scheme (Czekaj et al., 2014).

The atmosphere and track libraries are integral to the observable-space forward model. The preferred 2014 track package combined Bertelli et al. for intermediate and high masses with Chabrier et al. below $0.7\,M_\odot$, while BaSeL 3.1 was preferred over BaSeL 2.2 because it moved the giant/red peak in Tycho color distributions by more than 0.1 mag toward bluer colors, substantially improving agreement. Later BGM descriptions also emphasize BaSeL and NextGen as standard atmosphere or spectral libraries within the broader framework (Czekaj et al., 2014).

3. Dynamical self-consistency and distribution functions

A longstanding distinctiveness of BGM has been dynamical self-consistency. In the renewed thin-disc scheme, the total local stellar mass density is partitioned among age bins according to the SFR, living stars and remnants are separated by local sphere simulations, and the vertical structure is updated through the coupled Poisson and vertical-equilibrium equations,

$\nabla^{2} \Phi (r,z) = 4 \pi \rho(r,z) G$

and

$$\sigma_W^2 \ln \left(\frac{\rho(r,z)}{\rho(r,0)}\right) = -\Phi(r,z) + \Phi(r,0).$$

Any change in the IMF, SFR, local density, age–metallicity relation, or age–velocity relation therefore forces a recomputation of the mass model (Czekaj et al., 2014).

The 2018 axisymmetric dynamical version generalized this idea from local vertical consistency to a broadly self-consistent Galactic disc model. It assumed stationarity and axisymmetry, solved Poisson’s equation numerically on a cylindrical grid, and represented the stellar discs with generalized Shu distribution functions depending on the exact integrals $\Phi(m)$0 and $\Phi(m)$1 plus an approximate third integral derived from a Stäckel-potential approximation. The disc DF takes the form

$\Phi(m)$2

with exponential radial scalings for the DF density and dispersions. In that formulation, scale heights are no longer independently prescribed structural parameters; they emerge from the DF and the potential. A crucial validation result is that radial and vertical forces recovered through the Jeans equations agree with the model potential to better than one per cent over most of the Galactic volume (Bienaymé et al., 2018).

This DF-based reformulation also clarifies a known limitation of older Besançon-style Gaussian kinematics. When local Gaussian velocity ellipsoids were tested against RAVE and GCS data, they were found to describe the local GCS sample acceptably but to fail for the more spatially extended RAVE survey, especially in the $\Phi(m)$3 distribution and in the inferred radial decline of dispersions. A Shu DF handled non-circular orbits more accurately and provided a better fit, highlighting that purely Gaussian Besançon-style kinematic prescriptions are serviceable locally but inadequate for extended-disc dynamical inference (Sharma et al., 2014).

The Gaia-adjusted self-consistent dynamical BGM carried this program further by fitting disc DF parameters and the gravitational potential directly to Gaia eDR3 parallaxes and proper motions. In that model the disc DF is written schematically as

$\Phi(m)$4

with $\Phi(m)$5 for the thin-disc AVR. The resulting model yields a naturally flaring thin disc, a distinct thick disc with shorter radial scale length and different kinematics, and a nearly spherical dark halo, while remaining a smooth axisymmetric approximation rather than a model of the barred or spiral Galaxy (Robin et al., 2022).

4. Thin-disc renewal, IMF and SFR constraints, and fast inference

The 2014 renewal of the thin-disc BGM was centered on an all-sky comparison with Tycho-2 and the local luminosity function. Tycho-2 supplied a nearly all-sky bright-star catalogue complete to approximately $\Phi(m)$6, with 864,816 stars used after the completeness cut. Because its color distribution is bimodal—blue young BAF main-sequence stars and a red giant-dominated peak—it provided strong leverage on recent thin-disc star formation and on the intermediate- to high-mass IMF. The renewed BGM systematically explored thirteen ingredients or ingredient-groups, including IMF, SFR, thin-disc age, age–metallicity relation, age–velocity relation, local stellar and ISM densities, atmosphere library, tracks, binarity, thick-disc parameters, extinction, and radial scale length (Czekaj et al., 2014).

The principal quantitative conclusions were that the high-mass IMF slope is close to 3.0 and excludes a shallower slope such as Salpeter’s, that the thin-disc SFR is decreasing rather than constant, and that the model is compatible with a local dark matter density of about $\Phi(m)$7. Relative to the old model, the whole-sky star-count excess of roughly 50% was reduced dramatically; model A predicted about 8% fewer stars over the whole sky and model B about 4% more. Residual discrepancies remained largest in the plane, where extinction modeling and unmodeled structures such as nearby dust complexes, spiral structure, star-forming regions, and the Gould Belt were explicitly identified as limitations (Czekaj et al., 2014).

BGM FASt recast this renewed thin-disc framework into a big-data inference engine by reweighting a pre-sampled “Mother Simulation” rather than regenerating stars from scratch. Its core assumption is that the generated-star distribution can be approximated analytically so that star counts in a small region of parameter space are proportional to the stellar mass assigned there. The resulting fast approximation was reported to be about $\Phi(m)$8 times faster than standard BGM, enabling sequential Monte Carlo ABC inference on Tycho-2 color–magnitude diagrams. In its first application, BGM FASt again inferred a decreasing thin-disc SFH, with a present rate of $\Phi(m)$9, a local stellar mass density of $\Psi(t)$0, a local dark matter density of $\Psi(t)$1, and an intermediate-mass IMF slope $\Psi(t)$2 between $\Psi(t)$3 and $\Psi(t)$4; the high-mass slope remained strongly dependent on the extinction map (Mor et al., 2018).

5. Structural inference in the bulge, halo, and discs

In bulge studies, BGM has been used as a full forward model rather than as a background subtraction device. A 2MASS-based analysis of the central Galaxy fit triaxial Ferrers-like ellipsoids together with the inner disc and concluded that a single bulge ellipsoid is inadequate over the region $\Psi(t)$5, $\Psi(t)$6. The preferred model is a two-component structure: a main boxy bar with orientation about $\Psi(t)$7 relative to the Sun–Galactic-centre line and a second, thicker and longer ellipsoid. In that interpretation, the first component represents the main boxy bar, while the second may be either a classical bulge flattened by the bar potential or an inner thick-disc counterpart. The same framework also showed that a small flare of the bar can qualitatively reproduce the double red clump, and that a vertical metallicity gradient can arise naturally from changing proportions of a near-solar boxy bar and a more metal-poor thicker component (Robin et al., 2011).

For the stellar halo, BGM has been used as a photometric–spectroscopic forward model in the inner Galaxy. In an analysis combining 2MASS infrared color distributions with SDSS-III/APOGEE metallicities over 40 inner-Galaxy fields, the standard Besançon halo prescription was tested in the central few kiloparsecs using an axisymmetric oblate double power-law density law. The baseline model underpredicted stars with $\Psi(t)$8 dex, indicating that the outer-Galaxy-calibrated halo parameterization is inadequate when extrapolated inward and suggesting an extended inner metal-poor halo component or a needed revision of the halo’s chemical–spatial calibration (Fernández-Trincado et al., 2015).

BGM has also been used as a null model for chemo-kinematic interpretation of the Galactic discs. In a Gaia–APOGEE–Kepler analysis, mock BGM catalogues for the Kepler field were subjected to the same selection functions as the data. The model reproduced broad sample-selection effects and some global trends, but it failed to account for several salient features: the high-$\Psi(t)$9, metal-rich thick-disc sequence, the strong dependence of $i$0 on chemistry, and a maximum in $i$1 around 8 Gyr for the high-$i$2, metal-poor thick disc. Those discrepancies were interpreted as evidence that standard BGM assumptions omit important chemo-dynamical processes, especially mergers and radial migration (Lagarde et al., 2021).

6. Applications, performance domain, and limitations

BGM has repeatedly been used as an end-to-end survey forward model. In spiral-arm microlensing toward the EROS2 fields, it generated extinction- and efficiency-corrected synthetic source catalogues that reproduced observed color–magnitude diagrams and yielded optical depths and mean durations in reasonable agreement with the data, while also illustrating that catalog optical depth depends critically on the modeled source-distance distribution (Moniez et al., 2017). In MOA-II bulge microlensing, an updated Besançon model was extended to include M dwarfs and brown dwarfs down to $i$3; fitting the event timescale distribution required a brown-dwarf mass-function slope of $i$4, improved the average duration and the timescale-distribution shape, yet still left a serious low-latitude discrepancy, with the model providing only about 50% of the observed optical depth and event rate per star around the inner bulge at $i$5, which was interpreted as evidence for a missing inner stellar population or underestimated extinction and star counts (Awiphan et al., 2015).

The model’s scope has extended well beyond classical Galactic structure. A 2025 SETI application used BGM as a statistical completion engine for stellar “bycatch” inside radio-telescope beams, simulating 6,182,364 stellar objects across 1229 pointings out to 25 kpc and thereby replacing Gaia-only counts as a lower bound on the number of stellar systems probed. That work relied on the BGM web service and the model’s ability to supply stellar positions, distances, magnitudes, ages, masses, effective temperatures, and classifications, but it also highlighted concrete limits: BGM appeared to underpredict local counts within roughly 100 pc by about 30%, brown dwarfs were not modeled in that configuration, and extinction remained problematic in the Galactic Centre and spiral-arm tangents (Mason et al., 25 Nov 2025).

Several limitations recur across the literature. In the 2014 thin-disc renewal, white dwarfs were not yet fully integrated into the new generation scheme; extinction remained a dominant uncertainty in the Galactic plane; the thick disc was not re-optimized; and the parameter exploration, although systematic, was not a full global statistical inference. In the DF-based dynamical versions, the model is stationary and axisymmetric, so it does not represent the bar, spiral asymmetries, bending or breathing modes, or merger debris as dynamical structures. In bulge and halo applications, field-dependent extinction and the survey selection function remain major degeneracy sources. These limitations do not negate the model’s role; rather, they define the domain in which BGM is most effective: as a physically structured forward model of the smooth Milky Way, with explicit bridges between stellar populations, observables, and, increasingly, Galactic dynamics (Czekaj et al., 2014).