Towards a fully consistent Milky Way disk model V. The disk model for 4 – 14 kpc

,


Introduction
During the last few decades, extensive astrometric, photometric, and spectroscopic data on the Milky Way (MW) stellar populations have been collected.As a result, it became possible to build complex and realistic MW models and constrain their parameters by comparing their predictions to the observed stellar distributions, motions, chemical abundances, and ages.And though no exhaustive MW model exists at the moment, significant progress has been achieved towards understanding the MW evolutionary processes.
Relatively simple semi-analytic MW models, such as the TRI-dimensional modeL of thE GALaxy (TRILEGAL 1 ), a priori prescribe some observationally motivated spatial density distributions for the Galactic components.With the help of a stellar evolution library and additional assumptions about the MW 2 star-formation rate (SFR) history, initial mass function (IMF), and age-metallicity relation (AMR), these densities are release (DR2) of the Gaia mission (Gaia Collaboration 2016, 2018) and find evidence of an increase in star-formation (SF) activity in the last ∼2-3 Gyr of the MW disk history.However, the results of that work do not allow the authors to draw a firm conclusion about the SFR shape because of the limited age resolution of the model: there are only seven sub-populations representing the thin disk in the BGM, and this may be not enough to decipher with certainty the evolutionary complexity of the MW disk.
During the last decade, we have been developing a semianalytic model of the MW disk with a time resolution as fine as 25 Myr (Just & Jahreiß 2010, hereafter Paper I).So rather than having seven sub-populations as in the BGM, the Just-Jahreiß (JJ) model builds the thin disk from 520 mono-age subpopulations.This allows the disk structure and evolution to be studied in great detail, which also allows theoretical modelling to keep up with the ever-growing amount and quality of data accumulated for the MW stellar populations.
The JJ model is based on an iterative solving of the Poisson equation, and thus reconstructs a self-consistent pair of the overall vertical gravitational potential and density profile.The thin disk is modelled in the presence of gas, thick disk, stellar halo, and DM components.The disk evolution is governed by four input functions: the SFR, IMF, age-velocity dispersion relation (AVR), and AMR.So far, the JJ model has been tested many times and calibrated against different data in the solar neighbourhood.In Just et al. (2011, hereafter Paper II) we used star counts towards the north Galactic pole from the seventh data release (DR7) of the Sloan Digital Sky Survey (SDSS; Abazajian et al. 2009) to optimise the SFR and thick-disk parameters.Rybizki & Just (2015, hereafter Paper III) improved the IMF with the help of data provided by the HIgh Precision PARallax COllecting Satellite (Hipparcos; van Leeuwen 2007) complemented by an updated version of the Catalogue of Nearby Stars (CNS 4; Jahreiß et al. 1997).In Sysoliatina et al. (2018) we demonstrated a good model-to-data consistency (with a ∼7.3-% discrepancy in terms of star counts) in the local 1-kpc height cylinder using stars from the Radial Velocity Experiment survey (RAVE; Steinmetz et al. 2006;Kunder et al. 2017) and the Tycho-Gaia Astrometric Solution catalogue (TGAS; Michalik et al. 2015;Lindegren et al. 2016).Most recently, we used Gaia DR2 stars within 600 pc of the Sun and simultaneously optimised 22 model parameters within the Bayesian framework (Sysoliatina & Just 2021, hereafter Paper IV).We also calibrated the AMR against the local red clump (RC) sample from the Apache Point Observatory Galactic Evolution Experiment survey (APOGEE; Eisenstein et al. 2011;Bovy et al. 2014;Majewski et al. 2017).We have identified two recent SF burst episodes in the thin-disk evolutionary history, which happened ∼0.5 Gyr and ∼3 Gyr ago and left behind dynamically heated fossil populations.
In this paper we generalise our local JJ model for the range of galactocentric distances 4 kpc R 14 kpc by building the disk out of a set of independent radial annuli.We assume radially exponential gaseous, thin and thick disks, a spherical stellar halo, and a cored isothermal DM halo.We fix the disk thickness by assuming a constant-thickness thick disk and a constant-thickness or flaring thin disk and allow radial variations in the input functions SFR, AVR, and AMR, which correspond to the inside-out disk growth scenario.We present a publicly available tool, python package jjmodel, which can be used to perform stellar population synthesis within the JJ-model framework.Following Paper IV, we constrain the thin-and thick-disk AMR using the generalised JJ model with the radially extended We then propose a new analytic form for the thin-disk AMR, which is able to comprise the variation in its shape across the studied distances, 4 kpc < R < 14 kpc.Thus, we add another dimension to the fine JJ-model resolution: the numerous MW mono-age sub-populations can be now synthesised over a fine R-z grid.This paper has the following structure.In Sect. 2 we review the JJ-model components, explain how different ingredients of the local JJ model are now extended to other galactocentric distances, and summarise the overall scheme of mock sample modelling.Section 3 contains a description of a test model and its predictions for the MW disk properties.In Sect. 4 we reconstruct the radially dependent AMR using APOGEE RC data.Then, we discuss the obtained results for AMR, as well as the generalised JJ-model predictions and their possible applications, in Sect. 5. Finally, we summarise the main aspects of this work in Sect.6.

Generalised JJ model
In this section we introduce our key assumptions about the structure and properties of the Galaxy (Sect.2.1), recall stellar dynamic equations used for calculation of the Galactic potential (Sect.2.2), and review our treatment of the MW components in the local JJ model to generalise it for the whole MW disk ).We also briefly remind our approach to converting the modelled densities into mock samples (Sect.2.9) and present an overall scheme of building the MW disk within the JJ-model framework (Sect.2.10).

Basic assumptions and limitations
Before going into details, we briefly summarise the main assumptions on which the JJ model is based.We also mention limitations imposed by these assumptions on the model applicability (more in Sect.2.2).
First, the MW disk rests in a steady state.This statement is in harmony with the observational data, which show no evidence of fast evolutionary processes happening in our Galaxy at the present day (no major merger, no activity in the nucleus).All transformations of the MW populations' distribution functions occur slowly, as a result of collective interactions of the MW components (the so-called secular evolution).Thus, for our purpose of describing the MW at the present-day time snapshot, it is safe to treat the system as a steady one.
Second, the MW disk is axisymmetric and plane-symmetric (no bar, spiral arms, or warp are included).This implies that the JJ model cannot be used at small radii, R 4 kpc, where this assumption breaks because of a presence of the bar.Also, for comparison to the data, unless the model is used to reveal the underlying deviations from the smooth density distribution, the JJ-model predictions need to be averaged over large enough volumes to smear out small-scale density perturbations, such as overdensities in the spiral arms.Due to the absence of the warp, the predicted vertical structure can be reliable up to moderate radii, R 14 kpc.Third, the system is flattened (radial and vertical dynamics can be decoupled).Within the range of considered galactocentric distances, 4 kpc < R < 14 kpc, this holds up to |z| ≈ 2 kpc (see Sect. 2.2 for details).
Fourth, the MW stellar, gaseous, and DM components are represented by an individual or multiple isothermal subpopulations.Some of these sub-populations are assumed to be globally isothermal (stellar and DM halo), while others are isothermal only locally at fixed R (thin disk, thick disk, molecular and atomic gas).
Fifth, all sub-populations are in equilibrium with the total gravitational potential (Poisson and Jeans equations can be applied).This may not hold for very young populations, which did not have enough time to reach the equilibrium.Thus, the model cannot be applied to the youngest populations (O-and B-type stars).

Vertical profiles
At each galactocentric distance, the overall density in the JJ model is prescribed as a sum of i different MW components: stellar (thin disk, thick disk, stellar halo), gaseous (molecular and atomic gas), and DM4 : Here Φ is the total gravitational potential produced by this matter.For a single isothermal sub-population with the vertical velocity dispersion σ z,j , its density is linked to the total gravitational potential according to the Jeans equation: Here we dropped a time-derivative term assuming a steady state.
We also ignored two cross-terms containing σ Rz , essentially implying that the radial and vertical motions in the Galactic disk are decoupled.It can be shown that the ratio of these cross-terms to the rest of the terms in the Jeans equation is of the order of z/R • h j /R j , where h j and R j are the sub-population's scale height and scale length, respectively (Bahcall 1984;Bienayme et al. 1987).In a flattened MW-like system, h j R j , and therefore, within the range of radii 4 kpc < R < 14 kpc this ratio is small up to |z| ≈ 2 kpc.
The solution of Eq. ( 2) implies that the vertical density falloff of the isothermal sub-population is an exponential function of the total gravitational potential: Here potential is zero in the Galactic plane, and ρ j0 is the subpopulation's midplane density defined via its surface density Σ j and scale height h j : As a next step, we assume that all MW components in Eq. ( 1) can be represented by a single or multiple mono-age isothermal sub-populations, whose vertical density profiles follow Eq. ( 3).As a result, Eq. ( 1) can be turned into a specific function of the total potential with the vertical velocity dispersions and midplane densities of sub-populations as model parameters.Then, we can substitute this function into the Poisson equation and find a self-consistent potential-density pair.
For an axisymmetric system, the Poisson equation in cylindrical coordinates can be written as where G is the gravitational constant and F R = −∂Φ(R, z)/∂R is the radial force.It is common to assume that in flattened systems the second term on the right-hand side (the radial term) is not important in comparison to the first one (the density term), and thus, can be neglected (Binney & Tremaine 2008).Indeed, as a zero-order approximation, RF R can be replaced by its value in the plane, namely, by the squared circular velocity υ 2 c (R).As the MW rotation curve is nearly flat in a wide range of galactocentric distances, the radial term is vanishingly small and can be dropped.For the rotation curve reconstruction in Sysoliatina et al. (2018), we developed a more accurate treatment of the radial force that also includes dependence on z (see Eq. (4.19) in Sysoliatina 2018).Thus, we can estimate a ratio of the radial to density terms, in the framework of our MW mass model used in Sysoliatina et al. (2018).We find that in the range of galactocentric distances 4 kpc < R < 14 kpc, this ratio is ∼0.01 in the Galactic plane and may become important only at |z| ≈ 2 kpc in the inner disk, R 5 kpc, where it reaches ∼0.2-0.5.But even in this case, the radial term can still be neglected as it remains smaller than the density term uncertainty.The latter is dominated by the uncertainty of a poorly known DM density, which typically reaches tens of percent (Read 2014).
With these considerations taken into account, we neglect the radial term in Eq. ( 5) that significantly simplifies our problem.The Poisson equation without the radial term states that the vertical gravitational potential at some radius is produced only by matter at this radius, that is to say, the radial and vertical structure of the disk can be treated separately.This means that our approach of building the MW disk out of a set of independent radial annuli is justified.At each radius, the solution of the simplified Poisson equation reads as Here ρ(Φ) is the total density at this radius given by Eq. ( 1) with substitution from Eqs. (3) and (4).
In the following subsections, we describe in detail how each MW component is modelled over the considered range of galactocentric radii.

Thin disk
There are four analytic functions in the JJ model, which reflect main aspects of the MW thin-disk evolution and describe its A130, page 3 of 20 present-day state: IMF, SFR, AVR, and AMR.The thin disk has a present-day age of t p = 13 Gyr and consists of a set of locally isothermal sub-populations covering the whole age range with a fine resolution of t r = 25 Myr.This results into 520 thin-disk sub-populations, each characterised by its presentday mean age, metallicity, vertical velocity dispersion, and surface density.Below we explain how each of these quantities are prescribed , and how, together with the thickdisk and halo sub-populations, they are eventually converted into star counts (Sect.2.9).

IMF
For the IMF, we adopt a four-slope broken power law.Its up-todate parameters are listed in Paper IV (Table 4 and Sect.4.6.2) where we optimised them in consistency with the kinematic and SFR parameters using the local Gaia DR2 sample.This default JJ-model IMF is also applied during stellar population synthesis of the thick disk and halo (Sect.2.9).In the jjmodel package, the user can both change the IMF parameters and customise the model by defining another IMF law.
It is worth mentioning here that there is different observational evidence indicating that IMF is more top-heavy in the environments with low metallicity and high density (Bartko et al. 2010;Banerjee & Kroupa 2012;Kalari et al. 2018), and thus, also in the early Universe.On the other side, massive early-type galaxies seem to have a bottom-heavy IMF, and the same is reported for the MW low-metallicity stars in Hallakoun & Maoz (2021).A possibility for IMF varying in space and time is discussed in the literature (Kroupa et al. 2020) and tested with various methods -for example, in the cosmological simulations of the MW-like galaxies (Guszejnov et al. 2017) or with semi-analytic models of SF in pre-stellar cores of different metallicities (Tanaka et al. 2017(Tanaka et al. , 2018(Tanaka et al. , 2019)).However, most of the observed IMF variations are moderate and can be fitted into the universal IMF paradigm (Bastian et al. 2010;Offner et al. 2014).In our model we chose to stick to the conventional universal IMF for all disk sub-populations at all R and z.We do this for two reasons: (1) spatially and time-dependent IMF remains debatable, and (2) the regions we model have never belonged to the environments with extreme physical conditions, such as the Galactic centre, where the most significant deviations of the canonical IMF are observed.

AMR
The AMR function describes in a simple way the chemical evolution of the disk by linking the metallicity and age of its stellar sub-populations.In the previous JJ-model versions, we used for the thin-disk AMR a four-parameter logarithmic law (Eqs.( 21) and ( 22) in Paper IV).It describes a relatively quick chemical enrichment during the first ∼0.5-1Gyr of the disk evolution with a subsequent quasi-linear increase of metallicity with age in the solar neighbourhood.
There are two options for introducing the thin-disk AMR in the jjmodel code.More sophisticated and default AMR is described by a seven-parameter tanh-law, whose parameters are constants, linear or broken linear laws of R (Sect.4.2.4).By construction, this AMR is consistent with the observed metallicity distribution (MD) of the APOGEE RC stars across the disk, the reconstruction process of this AMR is described in Sect. 4.
The second and simpler way to prescribe the thin-disk AMR is analogous to how we extend the local thin-disk SFR to other radii (Sect.2.3.3):parameters of the four-parameter logarithmic AMR function with known values at R (Table 3 in Paper IV) are assumed to vary with radius according to power laws with predefined slopes.

SFR
Starting from Paper I, the default thin-disk SFR was a power-law function with a peak at ∼10 Gyr ago and a monotonous decline to the present day.In Paper IV, we found that the thin-disk SFR prescribed in such a way is inconsistent with the observed Gaia DR2 star counts of young populations in the local volume.The data-to-model consistency can be significantly improved by adding to the declining SFR continuum two Gaussian peaks centred at ages of ∼0.5 Gyr and ∼3 Gyr.This result is consistent with Mor et al. (2019) who find a single SF burst centred at the age of 2-3 Gyr.Also, Ruiz-Lara et al. ( 2020) report three SF enhancements 5.7 Gyr, 1.9 Gyr, and 1 Gyr ago, which were presumably triggered by recurrent passages of the Sagittarius dwarf galaxy.Sahlholdt et al. (2022) studied disk subgiants and dwarfs from the Third Data Release of the Galactic Archaeology with HER-MES (GALAH DR3; Buder et al. 2021) and Gaia Early Data Release 3 (EDR3; Gaia Collaboration 2021).The authors have also identified several enhancements in the disk SFR with the most recent SF peak at ages 2-4 Gyr, in agreement with the above mentioned studies.Motivated by these findings, we now define the thin-disk SFR as a standard smooth continuum that can be modified by an arbitrary number of overlaying Gaussian peaks, which allows the construction and testing of SFR functions of complicated shapes.
At a fixed radius, we define the thin-disk SFR as follows: It is a product of a normalisation constant and a time-dependent function; the latter consists of two terms that are related to the continuum and additional peaks (denoted by 'dc' and 'dp', respectively).The function is normalised to the overall presentday thin-disk surface density Σ d : A mass loss function g d (t) in Eq. ( 8) accounts for the cumulative effect of stellar evolution and is needed to convert the presentday stellar mass into the mass at the time of formation.In particular, g d (t) gives a fraction of a sub-population's mass survived to the present day in the form of stars (including remnants).This fraction depends on the IMF and the sub-population's metallicity (i.e., AMR; for g d (t) expression see Eq. ( 9) in Paper IV).We calculate these fractions with the default JJ-model IMF over the grid of metallicities and ages using the chemical evolution code Chempy5 (Rybizki et al. 2017).For each realisation of the JJ model, the mass loss function is then constructed in consistency with the chosen AMR by interpolation over the tabulated values.
As the mass loss function depends on the IMF very weakly, and an impact of the mass loss function itself on the modelled quantities is of a higher order than of SFR and AVR, we do not provide a tabulated g d (t) for IMF laws other than our standard four-slope broken power law.The continuum part of the thin-disk SFR is given by Eq. (10) in Paper IV, and for the sake of completeness we also reproduce it here: In this equation, parameters t d1 , t d2 , ζ, and η define the function's peak position and amplitude, as well as the steepness of the declining part.Following Paper IV, we define each extra SFR peak as a Gaussian weighted by the value of the SFR continuum at the peak's mean age and normalised to the continuum maximum value: where Here Σ p is related to a peak's amplitude at a given radius, and parameters τ p and dτ p give the mean age and the age dispersion of the peak, respectively.Galactic time t max corresponds to the peak of the SFR continuum Σ max (Eq.( 11) in Paper IV6 ).The summation is performed over all N p peaks.We note that this process of adding extra peaks to the SFR does not result in an increase of the overall present-day thin-disk surface density in the model.Instead, the SFR is re-normalised, such that the total surface density still equals to Σ d .We can also introduce quantities f k and f 0 defined as contributions of each extra SFR peak and the continuum part to the total SFR of a sub-population: Here s f r dp,k (t) corresponds to the kth term of Eq. ( 10).These fractions are later linked to the disk kinematics (Sect.2.3.4).
To introduce the radial variation of SFR, we allow some of its parameters to be functions of R. For the radial surface density profile, we assume an exponential law with a scale length R d laying in the range of ∼2.5-3 kpc (Fig. 1): where is the local surface density of the thin disk7 .Parameter t d1 is a time delay of SF and can be used to construct the model with non-overlapping in time thin-and thickdisk phases.Currently, we set this parameter to zero at all radii, this corresponds to the in-parallel formation of the thin and thick disk.Parameters t d2 , ζ, and η are now calibrated only at the solar radius R .As a first guess, we assumed for them a power-law radial dependence -a simple assumption that at the same time allows the inside-out disk growth to be mimicked.
We also assume that the thin-disk sub-populations corresponding to each of N p peaks form an overdensity normally distributed around some mean radii R p with dispersions dR p .Then, at an arbitrary radius, the amplitude parameter of the kth SFR peak can be calculated as where Σ p is the local value of the parameter.An example of the radially dependent SFR of the thin disk with two additional peaks is shown at the top panel of Fig. 2 (see also Sect.3.1).

AVR
The vertical kinematics of the thin-disk sub-populations is represented by a power-law AVR (Eq.(30) in Paper I).At a fixed radius, we prescribe it as a product of a normalisation constant, scaling parameter σ e , and an age-dependent function, such that σ W (τ) = σ e •avr(τ).We also link the present-day thin-disk velocity dispersion to kinematics of the molecular gas (Eq.( 14) in Paper IV).
Following Paper IV, we allow some of the thin-disk subpopulations to have kinematics different from the one prescribed by the AVR.More precisely, f k fractions of sub-populations associated with the kth SFR peak (Eq.( 11)) can have the vertical velocity dispersion σ p,k higher than the one prescribed by the AVR function for their age.The remaining f 0 fractions of the thin-disk sub-populations follow the AVR law.
To add the radial dependence, we allow one AVR parameter, its scaling factor σ e , to vary with R. Thus, the AVR shape is constant at all radii and is the same as in the solar neighbourhood.The value of σ e (R) is iteratively adapted at each radius to obtain the thin disk of a prescribed thickness (Sect.2.10).By default, the thin-disk thickness is constant at all radii, which is consistent with observations of edge-on galaxies (van der Kruit & Searle 1981;Bizyaev & Mitronova 2002).However, a flaring thin disk can be also tested with the jjmodel code, for that a simple exponential flaring model from Kalberla et al. (2014) is included.In this case, the thin-disk half-thickness is described by two parameters: radius at which the flaring starts, R f , and the flaring scale length R df .
If the thin-disk SFR is modelled with additional peaks and the peaks-related sub-populations are allowed to have special vertical kinematics, then at an arbitrary radius, velocity dispersion corresponding to the kth SFR peak is calculated as follows: where Here σ p ,k is the local value of the parameter, and σ ex,k represents the local velocity dispersion excess with respect to the AVR.The value of this excess is then quadratically added to the AVR velocity dispersion at each radius.
The bottom panel of Fig. 2 shows AVR corresponding to the JJ model with the thin disk of a constant thickness and SFR from the upper panel (Sect.3.1).Crosses mark velocity dispersion of the SFR peaks' populations.

Thick disk
The thick disk is modelled in the full analogy to the thin disk, but is described by a smaller number of parameters.
In the JJ model, the thick disk is constructed out of a set of mono-age locally isothermal sub-populations covering the first several Gyr of the MW evolution.The thick-disk SFR is a power-exponential function that peaks at old ages, ∼12.5 Gyr ago (Eqs.( 7) and ( 8) and Table 2 in Paper IV).Following Paper IV, we adopt here t t2 = 4 Gyr for the duration of the thick-disk formation, which corresponds to N t = 160 mono-age sub-populations.At a fixed radius, all thick-disk sub-populations have the same vertical velocity dispersion σ t , such that kinematically the thick disk behaves as a single isothermal subpopulation with the total surface density Σ t .The thick disk is also characterised by its own AMR function in a form of tanhlaw (Eqs.( 22) and ( 24) in Sect.4.2.4).
When modelling the thick disk at the different radii, we do not vary the shape of its SFR, but only normalise it to the corresponding surface density Σ t (R).For the radial profile of the thick-disk surface density we assume, by analogy to Eq. ( 12), an exponential decline with a scale length shorter than that of the thin disk, R t ≈ 2 kpc.The thick disk is prescribed to have a constant thickness, which we use as a constraint to calculate its velocity dispersion σ t (R) at the radii other than solar.The thickdisk AMR is also not varied with galactocentric distance in our model, which is consistent with the well-mixed thick-disk population observed in our Galaxy (e.g., Hayden et al. 2015Hayden et al. , 2020)).
When the JJ model is used to model volumes extending to heights larger than z max = 2 kpc (e.g., a cone perpendicular to the Galactic plane in Paper II and Paper IV), we extrapolate the derived thick-disk vertical profile with a sech α -law of a single isothermal sub-population.As we showed in Paper II, a simple sech α -law is still well-suited for the fitting of the thick-disk vertical profile calculated in the presence of DM and other stellar components.

Gas
Our current gas model for the solar neighbourhood is described in Paper IV, which we now extend to other galactocentric distances.
We consider two gas populations: a cold molecular (H 2 ) and a warm atomic (HI) component.The adopted local surface densities of the molecular and atomic gas are taken from McKee et al. (2015) and are Σ H 2 = 1.7 M pc −2 and Σ HI = 10.9 M pc −2 , respectively (Table 2 and Table 2 in Paper IV).For the surface density radial profiles of both gas components, we assume exponential laws.Kramer & Randall (2016) discuss the radial column and surface density profiles of the atomic and molecular gas as derived by different authors before 2006-2008 yr.An averaged surface density profile of the molecular gas peaks at ∼4 kpc and then declines with a short scale length of ∼2-3 kpc.We now adopt R H 2 = 3 kpc as a default value of the molecular gas scale length in the generalised JJ model.In the case of atomic gas, its surface density remains roughly constant in the wide distance range of 4-14 kpc and declines with a scale length of ∼4.5 kpc in the outer disk.Similarly, Kalberla & Dedes (2008) describe HI radial distribution: it is approximately constant in the inner disk and exponentially declining in the outer disk with a scale length of ∼3.75 kpc at 12.5 kpc < R < 30 kpc.We input R HI = 10 kpc to describe this slow decline of the atomic gas surface density with radius, which is broadly consistent with the mentioned observational data.
By analogy to the thin-disk mono-age populations, both gas components are treated as isothermal at each radius.Their vertical velocity dispersions are adapted to reproduce the prescribed scale heights of the cold and warm gas.We adopt the exponentially growing radial profiles of the H 2 and HI half-thickness from Nakanishi & Sofue (2016), which the authors measured using their three-dimensional gas maps constructed on the basis of HI and CO surveys data and kinematic distances to gas clouds.
In order to get a realistic enclosed mass of gas for the MW rotation curve, we add inner holes for both gas layers at R 0,H 2 and R 0,HI , which roughly reflects a sharp decrease of gas density in the inner disk.

Stellar halo
Stellar halo is modelled as a globally isothermal population of a single age t p .At each radius, its vertical density profile depends on two parameters, a universal vertical velocity dispersion σ sh and its surface density at this radius Σ sh (R).To get the radial surface density profile of the halo, we numerically integrate its power-law matter density profile up to z max .The default power-law slope is α in = −2.5 corresponding to the inner halo (Bland-Hawthorn & Gerhard 2016): Here Σ sh is the local surface density of the halo, and C is a normalisation constant, and Here we do not include halo flattening, as it will introduce only a small correction to Eq. ( 15).As the halo itself makes the smallest contribution to the total density among all other components, adding its flattening can be viewed as a secondorder correction to the model predictions, so we neglect it for now.Though, in some cases more comprehensive modelling of the halo is justified.In Paper IV, we simulated a vertically extended sample consisting of Gaia DR2 stars selected in a cone perpendicular to the Galactic plane.To get the right halo star counts, we extended the vertical halo profile that was derived self-consistently with potential only up to z max .For extending it to larger z, we used a broken power-law profile representing the inner and outer halo (Eq.( 15) in Paper IV).If necessary, flattening can be also added at this stage.
When creating mock samples, we assume for our single-age halo a Gaussian MD with a mean [Fe/H] sh and a dispersion d[Fe/H] sh .That is to say, the halo population is further subdivided into n sh = 5-9 sub-populations with different metallicities (see also Sect.2.9).

DM halo
The DM halo is treated in the JJ model similarly to the stellar halo.At each radius, DM component is added in the form a single isothermal population characterised by the vertical velocity dispersion σ dh , which is constant at all radii, and a radially dependent surface density Σ dh (R).
Assuming for the spatial DM density a cored isothermal sphere, we can express the radial DM surface density as follows: where Σ dh is the local DM surface density, a h is the DM scaling parameter, and r a h (R) = a 2 h + R 2 .The scaling parameter is adapted to have the local total circular velocity υ c, consistent with the observed proper motion of Sgr A * (Reid et al. 2005) at the assumed solar radius R , with the peculiar motion of the Sun V taken into account.

Bulge
Though we do not model the inner 4 kpc of the MW disk, the bulge presence needs to be taken into account for the calculation of the rotation curve.As we mentioned above in Sect.2.7, we use the predicted circular velocity at the solar radius to constrain the DM scaling parameter a h , which then enters Eq. ( 16).For this purpose, we include the bulge in a simple form of a point mass characterised by a single parameter, mass M b .As we do not add an inner hole in the thin-and thick-disk radial density profiles, M b gives a difference between the bulge mass and the enclosed mass of our exponential disks extrapolated to the inner MW region.

Creating mock samples
In order to convert modelled density profiles into the mock samples of some observable populations, we need to use a stellar evolution library.As we demonstrated in Sysoliatina et al. (2018) and Paper IV, the choice of the stellar library may have a significant impact on the location of different features on the colour-magnitude diagram (CMD) and on the overall predicted star counts introducing up to ∼10-% uncertainty.To estimate the impact of the stellar library choice on our modelling, we used in our previous works different sets of isochrones that we now include in the jjmodel package.
First, our default set of isochrones is generated with the PAdova and TRieste Stellar Evolution Code (PARSEC v.1.2S;Bressan et al. 2012) and the COLIBRI code (v.S_35; Marigo et al. 2017) for the thermally pulsing asymptotic giant branch phase 8 .Our isochrone grid includes photometric columns in the U BVRI JHK system and G BP , G RP , and G bands for the Gaia DR2 (Maíz Apellániz & Weiler 2018) and EDR3 (Riello et al. 2021).
Second, a complementary isochrone grid is based on the MESA (Modules and Experiments in Stellar Astrophysics, Paxton et al. 2011Paxton et al. , 2013Paxton et al. , 2015) ) Isochrones and Stellar Tracks (MIST 9 , v.1.2;Dotter 2016; Choi et al. 2016).For this stellar library, the following photometric bands are provided in the jjmodel code: standard U BVRI, 2 Micron All Sky Survey (2MASS; Skrutskie et al. 2006) JHK s , and three Gaia G-bands for DR2 and EDR3.
Third, another complementary grid is based on the Bag of Stellar Tracks and Isochrones (BaSTI 10 ; Hidalgo et al. 2018).The available photometric system is Gaia EDR3.
Following Paper IV, we use a wide-range metallicity grid −2.6 . . .+ 0.47 of 62 isochrone tables, each covering the whole modelled age interval 0 . . .t p with a linear step of 50 Myr.Each of the JJ-model mono-age isothermal sub-populations of the thin disk, thick disk, and stellar halo is characterised by an agemetallicity pair -according to the disks' AMR (Sects.2.3.2 and 2.4) and the halo MD (Sect.2.6).In the case of the thin and thick disk sub-populations, a spread in metallicity can be additionally allowed: for each age, we assumed a Gaussian MD centred at the metallicity prescribed for this age by the AMR function, [Fe/H](τ j ), and characterised by dispersion d[Fe/H] dt , which is the same for all ages.In practice, this means that we increase the number of modelled age-metallicity pairs by the number of sub-populations used to sample the assumed MD at each age, n dt .Though, this added complexity slows down the calculation process, it can help to mimic the effect of radial migration and produce more realistic Hess diagrams (Sysoliatina et al. 2018).In the end, the total number of the stellar monoage mono-metallicity sub-populations can be calculated as 8 http://stev.oapd.inaf.it/cgi-bin/cmd 9http://waps.cfa.harvard.edu/MIST/ 10http://basti-iac.oa-abruzzo.inaf.it/isocs.htmlA130, page 7 of 20 (N d + N t ) • n dt + n sh , where n dt = 1 if spread in metallicity for the disks is not allowed.
For each mono-age mono-metallicity sub-population, we select the best isochrone (the one closest in metallicity and age).This isochrone then provide us with a set of 'stellar assemblies' (SA) -populations of the same age, metallicity, and mass.The exact SA number, N SA , corresponding to a single isochrone depends on metallicity, age, and most of all, on the mass resolution of the stellar library.On average, N SA for a fixed radius is of the order of 10 5 -10 6 in the JJ model.
At the next step, we calculate for each SA its present-day surface number density according to the adopted IMF: Here m low l and m up l are the boundaries of the SA mass range, and M l is the overall mass converted into stars at the time interval corresponding to the SA age τ l .For the thin and thick disk, this mass is given by their SFR multiplied by the time resolution, M l = S FR(τ l ) t r ω l .Here ω l = 1 if there is no metallicity spread allowed for the mono-age sub-populations or equals to a corresponding Gaussian weight otherwise.For a single-age halo for which we do not have SFR, the initial mass converted into stars is calculated by dividing its present-day surface density by its mass loss factor, M l = Σ sh /g sh • ω l .The adopted fraction of the halo mass lost due to stellar evolution, g sh , corresponds to the adopted mean halo metallicity [Fe/H] sh and is driven from our pre-calculated mass-loss grid based on the Chempy code (Sect.2.3).By analogy to the thin and thick disk, weights ω l are Gaussian weights that originate from the halo MD sampling.
Following this scheme at each radius, we create SA tables for the thin disk, thick disk, and halo and store them for later use.As isochrones provide us with a set of useful stellar parameters (such as present-day stellar mass M f , luminosity log L, effective temperature log T eff , and surface gravity log g) and absolute magnitudes in different photometric bands, we can use this information to define observable samples (e.g., see the selection of the RC stars in Sect.4.2.2).For the SA subset corresponding to a mock population, we then calculate spatial number densities as a function of height at a fixed R: Here h(τ l ) and σ 2 W (τ l ) are the SA scale height determined during an iterative solving of the Poisson equation (Sect.2.10) and its vertical velocity dispersion.For the thin disk, σ W depends on age and is given by the AVR function.We note that if extra N p peaks are added to the thin-disk SFR, and the corresponding sub-populations are allowed to have special vertical kinematics (Sects.2.3.3 and 2.3.4), the right-hand side of Eq. ( 18) splits into N p + 1 terms.In the first term, which is weighted by the factor f 0 (Eq.( 11)), σ W is prescribed by the AVR.In the rest N p terms, weighted by factors f k , σ W is the velocity dispersion of the kth peak, σ p,k .The application of Eq. ( 18) to the thick-disk and halo SA is much simpler.In this case, the role of σ W play quantities σ t and σ sh .Also, velocity dispersions and scale heights of all thick-disk and halo SAs are independent of age by construction of our model.
With Eq. ( 18), prediction of star counts in some volume is straightforward.If the volume radial extension exceeds the adopted radial resolution of the JJ model, it can be split into several radial slices, in each of which Eq. ( 18) has to be applied separately.

Modelling scheme
At this point, we introduced all elements of the JJ model (Sects.2.1-2.9),so we can finally summarise the key steps of modelling process.To do so, we firstly rehearse physical parameters of the JJ model (Sect.2.10.1), and then list the sequence of steps required to create mock samples (Sect.2.10.2).

JJ-model parameters
The parameter file of the jjmodel package includes not only physical parameters of the model, but also some technical ones.The latter are boolean quantities that code special instructions (e.g., whether to add extra peaks to the thin-disk SFR or not, use a constant-thickness or flaring thin disk, build the local or radially extended model, etc.).The complete list of the JJ-model parameters with comments on their meaning, units, and usage can be found in the jjmodel documentation11 .Here we discuss only physical parameters, which we divide into several groups for convenience (see Table 2).
Solar parameters.We specify solar peculiar velocity V , its galactocentric distance R , and distance from the Galactic plane z .
R-z grid.The spatial grid of the JJ model is given by the centres of the inner and outer radial bins, R min and R max , bin width dR, maximum absolute height z max , and the vertical step dz.
Radial profiles.Another group of parameters prescribe the radial structure of the MW components.Scale lengths R d , R t , R H 2 , and R HI describe an exponential radial density fall-off of the thin disk, thick disk, molecular and atomic gas, respectively.For the gas components, we reserve two additional parameters R 0,H 2 and R 0,HI corresponding to the radii of the inner holes in their density profiles.Stellar halo radial surface density depends on the inner-halo power-law slope α in , and the point-mass bulge is described by its total mass M b .
Local normalisations.The overall radial density profiles of the modelled MW components are normalised to the local surface densities: Σ d , Σ t , Σ H 2 , Σ HI , Σ dh , and Σ sh .
SFR.The universal shape of thick-disk SFR is described by four parameterst t1 , t t2 , γ, and β.The continuum part of the local thin-disk SFR is given by t d1 , t d2 , ζ, and η.The corresponding power-law indices k t d2 , k ζ , and k η describe an extension of the local thin-disk SFR to other galactocentric distances.If there are N p additional Gaussian peaks added to the standard SFR of the thin disk, then each peak is characterised by the following five parameters: amplitude-related parameter Σ p , mean age τ p , age dispersion dτ p , mean radius R p , and radial dispersion dR p .
IMF.Our four-slope broken power-law IMF is defined by seven parameters -four slopes, α 0 , α 1 , α 2 , and α 3 , and three break-point masses, m 0 , m 1 , and m 2 .In order to work with a different IMF, one can introduce another set of IMF parameters in the jjmodel code.
and ( 23), an alternative set of parameters must be given (first row in Table 3 W-velocities.The thin-disk local kinematics is given by the AVR power index α and scaling factor σ e .Special kinematics of the sub-populations associated with each of the additional SFR peaks is described by W-velocity dispersion σ p .The local dispersions of the thick disk, DM, and stellar halo are given by parameters σ t , σ dh , and σ sh , respectively.

Building the disk
To construct the MW disk in the framework of the JJ model, we do the following sequence of steps.
Then, we calculate the input quantities, such as the radial surface density profiles of the MW components, IMF, local AVR, thin-and thick-disk AMR and SFR functions.At this step, we also optimise the DM scaling parameter a h using the predicted circular speed at the solar radius.
Finally, we proceed with the iterative solving of the Poisson equation at the solar radius.The solving procedure can be described in terms of three sub-steps.
First, the scale heights of isothermal sub-populations h j are initialised with realistic values (100 pc for the gas components, 400 pc and 1 kpc for the thin and thick disk, respectively, and 2 kpc for the DM and stellar halo).Additionally, we assume σ H 2 = 3 km s −1 and σ HI = 10 km s −1 for the vertical velocity dispersion of the molecular and atomic gas.
Second, we solve Eq. ( 6) with substitution of Eqs. ( 3) and (4) up to the maximum value of potential, Φ max , which is predetermined empirically at each R. The obtained dependence z(Φ) is then numerically inverted, and the result is fitted with a polynomial to get a smooth function Φ(z).
Third, using the derived vertical potential, we update all scale heights by integrating the normalised density profiles up to infinite z (Eq.( 8) in Paper I).In the case of the gas components, we also compare their derived scale heights to the shale heights from Nakanishi & Sofue (2016), which we aim to reproduce; the calculated difference is used as a weight for modifying σ H 2 and σ HI .Then the solving of the Poisson equation at the previous step is repeated with the updated scale heights and gas velocity dispersions as an input.The iterations continue until all differences between the old and new scale heights, as well as between the adopted and calculated scale heights of the molecular and atomic gas, become less than the vertical resolution of the model.With dz = 2 pc, the procedure converges after ∼5 iterations.
Using the obtained local density profiles, we calculate the overall thin-disk half-thickness at R .Then we prescribe its value also at the non-solar galactocentric distances (depending A130, page 9 of 20 on the model parameters, it can be a constant or flaring).Similarly, we fix the half-thickness of the thick disk at other radii by setting it everywhere to its local value.
At the next step, we solve the vertical Poisson equation at other radii.At R R , the solving procedure is applied with some modifications.Together with the gas velocity dispersions, we also assume the AVR scaling parameter σ e and the thickdisk velocity dispersion σ t .We also check the consistency of the derived thin-and thick-disk half-thickness with their prescribed values for this radius, and use the discrepancy values to modify σ e and σ t .For the iterations to stop, also the values of these discrepancies must become less than the vertical resolution dz.
All subsequent calculations can be classified as postprocessing of the obtained self-consistent vertical potential and density profiles.For example, we can calculate the vertical and radial density profiles of the mono-age disk sub-populations, which with the help of the AMR can be transformed into the profiles of the mono-abundance populations.We can study age, metallicity, or velocity distributions at the different R and z.Finally, IMF and isochrones can be used to create mock samples comparable to the real observable MW populations (Sect.2.9).

Results
To illustrate the predictive capability of the generalised JJ model, we present its test realisation, which is consistent with the local calibration from Paper IV.Here we must note that not all of the parameters defining the MW radial structure are calibrated against the data or are adopted from the literature; for several parameters we only assume plausible radial trends, which need to be tested with some observational sample in the future.

Test model
For the test model (TM1 hereafter), we adopt the following parameters (Tables 2 and 3).
The solar galactocentric distance is set to R = 8.17 kpc (GRAVITY Collaboration 2019), and its height above the Galactic plane is z = 20 pc (consistent with Bennett & Bovy 2019).For the solar peculiar velocity we take the classical value from Schönrich et al. (2010), V = 12.25 km s −1 .
Vertically we go up to 2 kpc with a fine 2-pc resolution.Radially, we consider galactocentric distances 4-14 kpc with 0.5-kpc bin width.This large radial step is sufficient for the illustration purposes of this section, but there are no physical or technical limitations for dR: it can be chosen to be as fine as required by the sample's modelling.
We model the thin disk of a constant thickness.The thinand thick-disk scale lengths are set to 2.5 kpc and 2.0 kpc, respectively, which is consistent with the different observational results summarised in Bland-Hawthorn & Gerhard (2016).Our choice of scale lengths of two gas components is described in Sect.2.5.The radial surface density profiles of both molecular and atomic gas have a cut-off at R 0,H 2 = R 0,HI = 4 kpc.For the slope of the stellar halo density profile we adopt α in = −2.5 (Bland-Hawthorn & Gerhard 2016), and the bulge mass is set to M b = 0.8 • 10 10 M .Our bulge component has a relatively small mass, because part of the b/p bulge connected to the bar is already included in the disk components.
The local surface densities, IMF, kinematic, thick-disk and local thin-disk SFR parameters are adopted from Paper IV. Parameters k td 2 , k ζ , and k η , which prescribe the radial variation of the thin-disk SFR continuum, are chosen such that the insideout disk growth is mimicked: the SF peak shifts from older to younger ages as we go from the inner to outer disk (upper panel of Fig. 2).Two extra Gaussian peaks at ages 0.5 Gyr and 3 Gyr are assumed to be centred at the radii of 7 kpc and 9 kpc with the radial dispersions of 0.5 kpc and 1 kpc, respectively.
The halo MD is characterised by two parameters, ([Fe/H] sh , d[Fe/H] sh ) = (−1.5,0.4) (consistent with Zuo et al. 2017).Finally, we adopt a new treatment of the thin-disk AMR based on the model calibration against the APOGEE RC sample (Sect.4).The thin-and thick-disk AMR parameters adopted for TM1 are listed in Table 3.We assume no scatter in metallicity for the disk sub-populations.

Predictions
We concentrate on the following model predictions: vertical and radial density profiles, scale heights of the mono-age and monoabundance thin-disk populations, and W-velocity distribution functions of several well-defined mock populations, such as G dwarfs and RC stars.

Radial and vertical density profiles
The top panel of Fig. 3 shows the vertical and radial thin-disk density profiles predicted for 0.25-Gyr age bins.The vertical profiles are calculated for the solar radius.The radial profiles correspond to the midplane densities.
The vertical density profiles for ages 4 Gyr demonstrate a monotonous decline with z.Though the youngest age bins are characterised by a more complex vertical density fall-off, we can identify several regimes with the different slopes.This behaviour is a direct consequence of the model input.In TM1, the monoage sub-population fractions associated with the extra SFR peaks are more dynamically heated than their remaining parts following AVR (see the bottom panel of Fig. 2).Thus, the youngest age bins consist of several kinematically different sub-components, and the final shape of the vertical profile depends on these subcomponents' fractions and W-velocity dispersions.We observed this multi-slope shape of the vertical density profile for the local A-type star sample selected from the Gaia DR2 (Paper IV).
As to the radial profiles, we see that though the overall thindisk surface density law is assumed to be exponential (Eq.( 12) and Fig. 1), the profiles of individual age bins can significantly deviate from this shape.The oldest thin-disk populations have quasi-exponential radial profiles, though two-slope exponential law can provide a more accurate fit in this case.Starting from ∼8-9 Gyr, profiles have a peak at the radius, which shifts to the outer disk as the age decreases.This reflects the lack of young populations in the inner disk due to the chosen SFR: young thindisk stars formed mostly at large radii.The density profiles of the youngest populations have additional overdensities located around the extra SFR peaks' mean radii R p .
This result is in a good agreement with the observed radial surface density trends in the MW disk.Bovy et al. (2016) and Mackereth et al. (2017) used APOGEE data to study the radial disk structure in terms of mono-age and mono-abundance populations.They find that the radial surface density profiles of both mono-age and mono-abundance low-α populations (which correspond to the chemically defined thin disk) are at best described by a broken exponential law.High-α populations representing the thick disk show no significant change of the slope over the studied distance range of 3-15 kpc.Similar radial structure follows from the JJ-model predictions (see the upper-right panel of Fig. 1, which shows radial surface density profiles of the mono-age thindisk sub-populations).As to the thick disk in the JJ model, its mono-age and mono-abundance sub-populations are all single exponential profiles with a slope of the overall thick-disk radial  profile shown in Fig. 1.Amôres et al. (2017) used BGM in combination with the 2MASS data to study the disk properties such as warp, flare, and populations' scale lengths.They find that the thindisk scale length is age dependent, with the youngest populations characterised by the largest scale lengths.The JJ model also recreates this trend which naturally follows from the inside-out disk growth scenario implied by the shape of our adopted SFR (Fig. 2).
The same negative correlation of scale length with age and doubleslope radial profiles of the thin-disk mono-age sub-populations are predicted by hydrodynamic simulations of the MW-like galaxies in Minchev et al. (2015).

Scale heights
The lower panel of Fig. 3 displays the radial change of the thindisk scale heights calculated for different age and metallicity bins.One of the interesting results of our modelling is that the mono-age populations flare even though the overall thickness of the thin disk is the same at all radii (dashed blue line in the lowerleft panel of Fig. 3).The same effect of mono-age sub-populations flaring when the overall disk thickness is approximately constant is reproduced in Minchev et al. (2015).The effect is generated by an interplay of two factors: (1) the inner disk is on average more dynamically heated than the outer disk, and (2) the relative fraction of old-to-young stars decreases with increasing radius.Similarly, Mackereth et al. (2017) find flaring for essentially all mono-age sub-populations of the low-α APOGEE stars.As follows from Fig. 3, for the youngest and oldest age bins in the range of 4-14 kpc scale heights change from approximately 80 pc to 180 pc and from 440 pc to 700 pc, respectively.Thus, flaring scale lengths predicted by this model lay in the range of ∼12.5-17 kpc.However, Mackereth et al. (2017) report a trend for the flaring with age that is the oppositive of that given by the JJ model: according to the APOGEE data, the youngest populations flare the most.Mackereth et al. (2017) suggest that at the early stage of the disk formation flaring could be suppressed by the active accretion.This topic may be deeper investigated in our future works.Another feature is the presence of a peak in the scaleheight profiles of the young populations due to the contribution of kinematically hot SFR peaks' sub-populations.Also, we see an opposite 'response' in the scale heights of older ages around 9 kpc: here scale heights grow slower than we would expect from the trends up to this radius.We can attribute this behaviour to the SFR re-normalisation.As explained in Sect.2.3.3,we renormalise the total thin-disk SFR to the same Σ d after adding extra peaks.In the case of TM1 we have two extra peaks at young ages, that is, we increase the contribution of the young sub-populations by decreasing contributions of all other ages.And though these young sub-populations are partly more heated than prescribed by AVR, the modified age balance may in some age bins result in a more plane-concentrated matter distribution than we could expect without extra SFR peaks.
As mono-age populations are not easily observable, it is useful to convert them to mono-abundance populations using the AMR function (Sect.4.2.4).At the lower right panel of Fig. 3 we show the radial variation of the thin-disk scale heights A130, page 11 of 20 calculated in 10 metallicity bins on the interval from −0.7 to 0.3 dex with a 0.1-dex bin width.We see that the scale-height radial profiles of the metal-poor populations closely follow the corresponding profiles of the oldest age bins.In other words, mono-abundance thin-disk populations can be also considered as mono-age when the metallicity is low.To the same conclusion arrived Minchev et al. (2017) who investigated the relationship between mono-age and mono-abundance sub-populations using the chemodynamical model of the MW-like galaxy with inside-out disk growth.As the metallicity increases, monoabundance populations begin to include more ages, and the corresponding age distributions depend on radius.For metallicities >−0.2 dex, profiles do not span over the whole range of galactocentric distances, but end at some radius, and this stop point moves to the inner disk as the metallicity increases.This happens because a given metallicity ceases to be present starting from some radius, according to our simple chemical evolution model.Due to the multi-age composition of most of the mono-abundance populations, the flaring trend is observed only for most metal-poor bins, while for other mono-abundance subpopulations it is suppressed.Also, in Minchev et al. (2017) the authors find that flaring is not prominent in most of the metallicity bins because of their negative radial age gradient.On the other side, radial surface density profiles of APOGEE monoabundance sub-populations show flaring in almost all metallicity bins (Bovy et al. 2016;Mackereth et al. 2017).This can be viewed as a lack of flaring in the overall thin-disk profile, which can be easily included into the model.Also, this can point to the lack of some fundamental physical process, such as radial migration, that would relocate some fraction of stars from the inner to the outer disk and thus lead to an extra flaring.
As to the thick disk, our non-flaring disk is consistent with the observed nearly constant thickness of the high-α APOGEE mono-abundance sub-populations (Bovy et al. 2016;Mackereth et al. 2017), properties of the thick disk constrained with the SDSS data in the BGM framework (Robin et al. 2014), and the thick disk simulated in Minchev et al. (2015).

Vertical kinematics
Following a routine described in Sect.2.9, we converted the modelled densities into SA using PARSEC isochrones.Then we created several mock samples of the well-defined observable populations, such as A, F, and RC stars and G and K dwarfs.All of them, except RC stars, were selected by simple cuts applied to the effective temperature and surface gravity (in the case of dwarfs).Our RC modelling routine is described in Sect.4.2.2.
Figure 4 shows W-velocity distribution functions of these five mock thin-disk samples calculated for the solar radius for |z| < 2 kpc.As expected, A-and F-star samples consisting exclusively or predominantly from the young populations, have f (|W|) with a prominent core.The main sequence (MS) stars, G and K dwarfs, have velocity distributions with a much more significant contribution of the dynamically heated stars.Interestingly, the RC sample consisting of a mixture of different ages with a significant fraction of young populations, has f (|W|) essentially indistinguishable from the one of the MS samples (also see Figs. 8 and 10 in Paper IV).For the current model realisation, this can be understood as follows: though the RC sample contains more young stars than G and K dwarfs, young populations of the SFR peaks also partly contribute to the tail of the velocity distribution, such that the resulting RC and G-and K-dwarf f (|W|) are very similar.Figure 4 is also representative for other modelled radii.

AMR from APOGEE RC stars
In this section, we take the first step towards calibration of the generalised JJ model across the disk.We use the developed machinery to constrain the AMR function with the APOGEE data covering distances 4 kpc < R < 14 kpc.Motivated by the obtained results, we then develop a new analytic representation of the thin-disk AMR.

Sample selection
To constrain the AMR across a wide range of galactocentric distances, we use the latest data release, DR16, of the APOGEE RC catalogue (Bovy et al. 2014;Ahumada et al. 2020).The APOGEE RC catalogue comprises a kinematically unbiased sample of 39 675 bright stars spanning over a large fraction of the Galactic disk.As the luminosity of RC stars depends very little on their chemical composition and age, they can play the role of standard candles.Thus, for the RC stars, it is possible to obtain robust photometric distances -distance uncertainties of the APOGEE RC stars are reported to be less than ∼5% (Bovy et al. 2014).Using heliocentric distances provided in the RC catalogue, we calculated Galactic cylindrical coordinates R and z (the adopted position of the Sun (R , z ) is the same as in Table 2).
It is well known that when projected onto [α/Fe]-[Fe/H] chemical abundance plane, the MW disk stars form two distinct sequences, which are associated with the thin-and thick-disk populations.We used a simple criterion to separate low-and high-α sequences, similar to the one adopted in Paper IV12 (Fig. 5, top panel): To eliminate halo stars, we applied a cut [Fe/H] > −1.From the remaining subsample, we took stars with |z| < 2 kpc, both for the high-and low-α population.Finally, we defined the radial range  19) is shown with dashed black line.Bottom: selected low-α and high-α RC populations in the R-z plane.Dashed lines mark the adopted limits for the galactocentric distance and height.The radial binning of the low-α population corresponds to run0 (see text).The adopted position of the Sun is marked with a yellow dot.
of the samples.Here we keep in mind that the high-α population is known to be chemically homogeneous, while the low-α MD systematically change with galactocentric distance.Therefore, in case of high-α sample, we selected stars with galactocentric distances 4 kpc < R < 14 kpc and did not apply radial binning (red dashed box at the lower panel of Fig. 5).For the low-α stars, we used a grid of 1-kpc bins centred at 4.5, 5.5, . . ., 13.5 kpc (blue dashed lines).We later refer to this grid as run0.In order to eliminate an impact of binning on the reconstructed AMR, we used four additional grids for the low-α sample (Sect.4.2): three grids shifted with respect to the initial one by +0.25, +0.5, and +0.75 kpc (run1, run2, and run3, respectively) and a finer grid with 0.5-kpc bins centred at 4.25, 4.75, . . ., 13.75 kpc (run4).Also, we found that bins of these grids that partially lay at R > 14 kpc and contain fewer than ∼300 stars cannot provide enough data to set a robust constraint on the thin-disk AMR parameters.For this reason, we excluded such bins from our analysis.
In the selected samples, 95% of stars have a signal-to-noise ratio (S /N) 50, so we did not apply any quality cuts.Our final RC samples contain 4444 (high-α) and 32 917 stars (low-α, run0).
It is also worthwhile to address the problem of extinction, which is known to be very significant in the Galactic plane at the large distances reached by our sample.As explained in Majewski et al. (2011), the APOGEE photometry is corrected for extinction with the help of the Rayleigh-Jeans colour excess (RJCE) method.The RJCE technique is based on sampling the spectral energy distribution (SED) in the Rayleigh-Jeans regime where the shape of stellar SED is essentially independent of its spectral type.Thus, extinction for each star can be robustly estimated by measuring the colour excess from its photometry in near-and mid-infrared bands.As the RC distances are based on the APOGEE photometry corrected for extinction, they are not biased.As reported in Bovy et al. (2014), uncertainties in extinction determination make only a minor contribution to the total error budget of the RC distances.The stellar parameters that were used to select the RC stars from the full APOGEE catalogue, such as metallicity, log g, T eff , were derived with the APOGEE Stellar Parameter and Chemical Abundances Pipeline (ASPCAP).The effect of reddening is suppressed in ASPCAP by working with the observed and template spectra that are continuum-normalised (García Pérez et al. 2016).Taking all this into account, we can conclude that no special measures need to be taken in our analysis to account for the extinction, and we can compare the APOGEE RC sample to our mock RC population simulated with the stellar parameters and distances not affected by extinction.

AMR reconstruction
Our approach for the AMR reconstruction is based on a direct comparison of the observed MDs to the modelled age distributions of some population and was explained in detail in Paper IV.
Here we summarise its main steps as applied to this study.
First, we calculated the observed cumulative metallicity distribution functions (CMDFs) of the low-and high-α stars using the selected RC data samples.In the case of the low-α population, CMDF was obtained separately in the different radial bins.
Then assuming AMR for the thick disk (same at all radii) and for the thin disk (radially dependent) and working in the framework of the generalised JJ model, we derived the total vertical gravitational potential and the vertical density profile of the RC population at each radius.The density profile was then used to calculate the cumulative age distribution functions (CADFs) of the thin-and thick-disk RC stars within the selected R-z limits.
At the next step we compared the observed CMDFs to the predicted CADFs: we put in correspondence stellar metallicities and ages, and thus, arrived at the updated AMR of the thin-and thick-disk (radially dependent and the same for all R, respectively).To achieve self-consistency, we iterated over the last two steps until the obtained thin-and thick-disk AMR converged to a constant shape (∼3-5 iterations).
After that we fit the obtained thin-and thick-disk AMR with analytic functions.In the case of the thin disk, we investigated the radial trends of the AMR best-fit parameters and also fit them with simple laws.As a result, we obtained a 'generalised' thindisk AMR fit, which is well-consistent with the originally reconstructed 'raw' AMR.
Finally, we performed a consistency test: calculated the thinand thick-disk MDs of the RC stars using the updated AMR and compared these predictions to the data.
Each of these steps is described in more detail in the successive subsections.

Cumulative metallicity distributions
The CMDFs of the selected low-and high-α APOGEE RC stars are shown in Fig. 6.There are three distinct regimes that can be identified for all CMDFs: a rapid non-linear increase at low metallicities, linear part, and a non-linear deceleration in growth at high metallicities.Following Paper IV, we make a simplifying assumption: at each radius, chemical enrichment process was monotonous in time (i.e. the AMR of interest is a monotonously growing function of Galactic time).This implies that low-and high-metallicity tails of CMDFs correspond to old and young populations, respectively.Then the non-linear growth of a CMDF at low metallicities reflects a quick chemical enrichment at the beginning of the MW disk evolution when the SF processes, and therefore the enrichment of ISM with metals, were most intense.The non-linear shape of CMDF at the highmetallicity end is, however, puzzling in this framework, as we do not expect a significant increase in the ISM chemical enrichment rate at the present day and do not have observational data that imply that.So, in order to obtain the AMR that does not show any special behaviour at the present day, we represented the upper part of each CMDF as a convolution of the linear trend (monotonous chemical enrichment) with a Gaussian kernel (error model).The analytic form of such a convolution was given in Paper IV (Eq.( 17)), and for the sake of completeness of narration, we show it here: where Parameters k and b are slope and intercept of the linear part of CMDF, and σ [Fe/H] is a dispersion of the Gaussian kernel.
For each CMDF, we selected the linear part and obtained parameters k and b (0.4 < N cum /N tot < 0.8 for the low-α stars with bin centres at R < 7 kpc and 0.3 < N cum /N tot < 0.7 at all other radii and for the high-α population).Then, with k and b known, we fit the upper parts of CMDFs with Eq. ( 20) and thus determined σ [Fe/H] -dispersion required to reproduce the high-metallicity CMDF tail with the Gaussian error model.The corresponding values of σ [Fe/H] for the different R bins of the low-α population are displayed at the upper panel of Fig. 6.Later we used these values to model the high-metallicity MD parts.The 'deconvolved' upper parts of CMDFs are shown with dotted lines in Fig. 6.
The error model can represent different aspects (or be a result of their interplay): (1) systematic errors at the metal-rich end or (2) the presence of the real stars that have migrated to a given volume from the inner radii and therefore have an imprint of the different environment and probably do not belong to the youngest population.Which of these interpretations is correct, and to which extent, we do not investigate within our approach.

RC star modelling
TM0 framework.For the AMR reconstruction, we used a simpler version of the test realisation of the generalised JJ model presented in Sect.3.1.The most uncertain ingredient of TM1 is the extra SFR peaks at R R .In order not to propagate this uncertainty to our results and also following Paper IV, we omitted the additional SFR peaks for the RC sample modelling.To initialise AMR at the second step of the AMR reconstruction scheme, we used the thin-and thick-disk AMR from Paper IV calibrated against the local APOGEE RC sample (Eqs.( 21) and ( 22)).The thin-disk AMR was extended to other radii by assuming power-law radial variations of its four parameters.The power indices were chosen such that the radially dependent AMR predicted quicker and more intense chemical enrichment in the inner disk: RC mock sample.To select the RC stars in TM0, we used the following criteria for the effective temperature, surface gravity, and luminosity: 4250 K < T eff < 5250 K, log g < 2.75, 1.65 < log (L/L ) < 1.85.These cuts determine a sample of RC contaminated by red giant branch (RGB) stars.To remove RGB stars, Bovy et al. (2014) use an additional sloping cut in log g, a linear function of the effective temperature and metallicity: log g < f (T eff , [Fe/H]).As we have already found in Paper IV (Eqs.( 18)-( 20)), the exact values of this cut's slope and intercept depend on the stellar evolution library, and also on the population's MD (i.e., on the assumed AMR).When visualised in (log g, log T eff , log L)-space, the pre-selected sample forms a fork-like structure, with the two hands being RC and RGB sequences.The f (T eff , [Fe/H]) criterion defines a plane that separates one sequence from another.In order to avoid the need to manually adapt the f (T eff , [Fe/H]) criterion each time we change the model parameters or isochrones, we introduced a numerical routine to separate the RC and RGB stars (included into the jjmodel code).After applying our routine to the pre-selected sample, we got a clean RC population.At the subsequent steps, we modified the vertical spatial densities of A130, page 14 of 20 this complete sample to account for the geometry of the observational data (Sect.4.2.3) and then used these vertical density profiles for the AMR reconstruction (Sect.4.2.4).

Cumulative age distributions
To perform realistic modelling of some stellar population, it is important to properly take into account possible selection effects, so we address this question before modelling RC CADFs.
The R-z plot of our APOGEE samples in Fig. 5 (bottom) shows that the selected RC stars are distributed in space very unevenly.Indeed, the APOGEE survey does not provide a full sky coverage but targets a set of 'grid' and 'non-grid' pointings, whose selection is motivated by location of objects of special interest and goals of the different scientific projects, as well as some technical limitations (Zasowski et al. 2013(Zasowski et al. , 2017)).To observe spectra, fibre-plugged plates are used, such that each fibre corresponds to a certain area in the sky.As a result, space coverage of the APOGEE consists of so-called pencil beams probing the MW stellar populations within the selected directions (Wilson et al. 2019).
We can formulate this spatial selection in the following way: at each radius, the observed vertical number density law deviates from the real one, and the residual reflects the sample geometry (assuming that there are no other losses).In the JJ model, on the contrary, we by default work with complete populations whose vertical number density declines smoothly with z.So we need to understand whether the spatial irregularity of the APOGEE RC sample is relevant for the calculation of CADFs, and in which form it can be taken into account.
The thick disk is known to be a well-mixed populationquite homogeneous in age, chemical composition, and kinematics (Lee et al. 2011;Haywood et al. 2013;Hayden et al. 2015Hayden et al. , 2020)).Within the JJ-model framework, it is represented with a single isothermal population at each R.This implies that the predicted thick-disk age and MDs are the same at all heights at a given radius.Thus, even if the predicted thick-disk vertical number density law (complete sample) does not exactly follow the observed vertical distribution ('pencil beams'), this cannot introduce any bias to the modelled age distribution, and therefore, no correction for the sample geometry is needed.
In the case of the thin disk, the situation is different.As we represent the thin disk with multiple isothermal populations, its kinematic, chemical, and age properties significantly change with z at a fixed R. Therefore, to model the realistic age distributions of the thin-disk APOGEE RC sample, it is important to properly account for the sample's spatial geometry.
To do so, we use the following approach.At each radius, we calculated N data (|z|), a distribution of the observed thindisk RC stars over |z| up to 2 kpc with a step of 50 pc.In the model, we do not model the real star counts, but predict the vertical number density of RC stars according to Eq. ( 18), N V model (|z|).To be able to compare modelled and observed quantities, we normalised both vertical laws to their midplane values.As N V model (|z|)/N V model (0) = N model (|z|)/N model (0), we arrived at the comparable quantities.The observed and predicted normalised vertical profiles are shown at the top panel of Fig. 7.We calculated the ratio: At each |z|, quantity ω z shows by how much the observed vertical distribution of stars deviates from the vertical density law of a complete population.Assuming that TM0 predictions are entirely correct, this deviation reflects the observed sample's incompleteness and geometry.Thus, we can use ω z to modify the modelled RC age distributions to mimic the geometry and incompleteness of the APOGEE RC sample.
For the reasons explained below in Sect.4.2.4,we find it useful to test an impact of the warp on our modelling.We adopted a simple linear warp model from Chen et al. (2019).The normalised vertical number density laws and the relative weights ω z corresponding to TM0 with and without the warp are shown in Fig. 7.The mock RC CADFs calculated with TM0 without the warp and weights ω z taken into account are shown in Fig. 8.The displayed CADFs correspond to the 5th iteration (consistent with the reconstructed AMR in Sect.4.2.4).

AMR fitting process
In Paper I-Paper IV, we used a logarithmic function to describe a monotonous increase of the thin-disk sub-populations' metallicity with time, and in Paper IV we proposed a tanh-law for the AMR of the thick disk.In this study, we find that the shape of the thin-disk AMR is very different at the different R, ranging from a function with a steep increase at small t in the inner disk to a quasi-linear law in the outer disk (Fig. 9).Though our old loglaw is well-suited for reproducing the local thin-disk AMR, it is unable to describe such a family of curves.Therefore, for the radially dependent AMR of the thin disk, we introduce a new seven-parameter tanh-function, which can be viewed as a more complex form of the thick-disk tanh-law.As before, the metallicity as a function of time for both thin and thick disk is expressed in the following form: where functions f (t) for the thin and thick disk are given by Eqs. ( 23) and ( 24) 13 , respectively: Figure 9 shows the thin-and thick-disk AMR (solid curves), which we reconstructed following the procedure described in the previous sections.We note that AMR not only of the thick disk, whose SF ceased in the model ∼9 Gyr ago, but also of the thin disk does not extend to the present day.In the case of the thin disk, this reflects the age span of the modelled RC sample: depending on the radius, the last ∼0.5-5Gyr of the thin-disk evolution are not represented by our sample (see age distributions in Fig. 8).We recovered the missing ages by fitting the reconstructed raw AMR and extrapolating these fits to the present day.
Fitting the thick-disk AMR with Eqs. ( 22) and ( 24) is straightforward.The obtained parameters are listed in Table 3, and the fitted AMR is shown in Fig. 9.
The analytic thin-disk AMR has to be a function of time t and galactocentric distance R. The presented thin-disk AMR equation, Eqs. ( 22) and ( 23), is a function of time only, as it describes the AMR shape at some fixed R. To add dependence   22)-( 24), parameters from Table 3).
Next, we fit parameter α ω with a linear profile and updated the remaining (t 01 , t 02 ).
Finally, we adopted fits for parameters t 01 and t 02 .Two bottom panels at Fig. 10 show their best-fit values at different R. It is clear that both of these parameters change with radius nonmonotonously.At the same time, the value of metallicity predicted according to Eqs. ( 22) and ( 23) is very sensitive to t 01 and t 02 .As a result, it is not possible to fit the obtained data points with a constant or linear trend as we did with [Fe/H] 0 and α ω and obtain a good fit of the raw AMR.For this reason, for both t 01 and t 02 we used three-slope broken linear laws.
Figure 10 shows the best values of five parameters of the thin-disk AMR, which we obtained using such a step-by-step fitting routine.We repeated this fitting procedure for our five R grids.Using different binning results in a scatter in the parameters' radial trends, but overall the data points of all runs are well consistent with one another; therefore the generalised AMR that we defined based on these profiles is not biased by the choice of bin edges.
There is a distinct change of slope in both t 01 and t 02 trends around ∼10 kpc.In order to check whether this feature partially or fully originates from the absence of the warp in our model, we tested TM0 with the warp.It is already clear from Fig. 7 that the impact of warp starts to be noticeable only at R ≈ 13 kpc and is very small.The best-fit values corresponding to the AMR reconstruction with TM0 with the warp are essentially indistinguishable from the values obtained for the runs without the warp at Fig. 10.This essentially does not change if we use the model with the more pronounced warp from Poggio et al. (2020).From this we can at least conclude that the non-monotonous behaviour of t 01 and t 02 radial profiles in the outer disk is not related to the lack of warp in our model.
As a wrap-over, we applied the step-by-step fitting to the parameter values of run0-run4 collectively and obtained the final radial profiles of the AMR parameters (green lines on Fig. 10).At this stage, we also optimised the break-point positions for t 01 and t 02 by minimising a mean relative deviation between the real and fitted best parameter values: Here f t 0i are the analytic profiles with some assumed break-point positions, σ t 0i is the dispersion.The calculation was performed in the R grid with bin centres at 4.5, 5.5, . . ., 14.5 kpc.We find that the optimal position of the second break is 9.75 kpc for both parameters, while the first break lies at 7.5 kpc and 6.25 kpc for t 01 and t 02 , respectively.The break point at 6.25 kpc is very close to R br1 = 6 kpc, which we adopted for [Fe/H] p , as this parameter is strongly correlated with t 02 (both affect the AMR shape at young ages).In order not to further complexify our analytic representation of the thin-disk AMR, we set the first break of t 02 profile also to 6.0 kpc.The obtained fits for the AMR parameters are listed in Table 3.The resulting generalised thin-disk AMR is plotted in Fig. 9 with dashed coloured curves.We see that the final AMR fit has a very good consistency with the raw reconstructed AMR.

Consistency test
To demonstrate the consistency between the APOGEE data and the JJ model after the AMR calibration, we predicted the RC MDs using the derived AMR and compared them to the observed metallicities.The MD modelling was performed in the full consistency with the approach used to reconstruct the AMR function.For the thick disk, MD was predicted for the mock RC population within the range of 4-14 kpc and at |z| < 2 kpc in the TM0 framework.For the thin disk, MDs were calculated separately in each radial bin, and the vertical densities of the modelled RC populations were weighted with ω z .We also added an error model as a post-processing step: each MD was converted into cumulative form, convolved with a Gaussian kernel of σ [Fe/H] as derived from the data (Fig. 6), and then transformed back to non-cumulative distribution.The resulting thinand thick-disk MDs are shown in Figs.11 and 12. Two types of step histograms correspond to TM0 predictions with the raw and generalised analytic AMR.
We see that the model and data are well consistent in terms of the overall shapes and peak positions of the distributions.The maximum shift that we report between the modelled and observed MD is only ∼0.05 dex.

Discussion
The goal of the JJ model is to provide a flexible modelling tool for the main body of the MW disk in terms of stellar populations and their properties.We model the axisymmetric disk in dynamical equilibrium.In order to allow the investigation of the A130, page 17 of 20 Galaxy evolution over the full age range, we represent the stellar content of the disk by sequences of stellar sub-populations with a high time resolution of 25 Myr.This ansatz is complementary to the BGM model, which tries to incorporate all MW features including the bar, spiral structure, and the warp at the cost of a restricted time resolution.The input JJ-model functions SFR, AVR, and AMR are relatively simple functions of galactocentric distance and time (i.e.age).The IMF is assumed to be universal, but this could easily be made a function of metallicity in the framework of the JJ model (see Sect. 2.3.1).

Model properties and limitations
The JJ model provides predictions for the density distribution, vertical velocity distribution function, age distributions, and MD function for stellar subsamples (along the MS, RC stars, etc.).The parameters of the input functions are constrained by comparison of these predictions with corresponding observed samples.We did not use observed samples with individual age determinations, because they are still very restricted and may have significant biases.Therefore, direct age determinations can be used as independent tests of the model (or our model can be used as a tool for deriving age distributions of stellar populations).We only remark here that to perform a robust model-todata comparison, the impact of dust must be taken into account: a three-dimensional extinction map can be used to redden our synthetic photometry to compare to the data (forward modelling, Sysoliatina et al. 2018), or the data can be preliminary de-reddened and then compared to the model (as in Paper IV). Also, the selection function of the sample needs to be thoroughly scrutinised to include the corresponding selection effects to the model.Since the JJ model contains only axisymmetric populations in dynamical equilibrium, it can be used to characterise the properties of model-to-data deviations, such as the spiral arms and the warp.Also, due to the modular design of the code, simple spiral arm or warp models can be implemented.As an example, we tested the impact of the warp on the AMR derivation and found no significant differences (see Sect. 4.2.4).
The SFR, AVR, and AMR describe the present-day properties of stellar sub-populations at the current galactocentric distance.All three functions depend on each other for the derivation of the present-day distribution functions to be consistent with observations.Also the stellar birth places do not appear explicitly in our model.Only in the limiting case of negligible radial migration (blurring and churning), the time-and radially dependent SFR is the star formation history at radius R, the AVR represents the heating function, and the AMR corresponds to the metal enrichment law.In the case of a significant radial migration impact, these input functions can be viewed as the results of a convolution of the true SFR, AVR, and AMR with the radial migration model.In Paper IV, we addressed this topic when discussing physical meaning of the reconstructed local SFR of the thin disk.Assuming that the impact of migration is a function of time, we can formulate this the other way: specify a range of radii that are represented by a particular age interval of the SFR, AVR, or AMR.
Additionally, the JJ model can be combined with the Chempy code to derive the chemical enrichment history at each radius including the gas infall and abundances of other elements.

AMR and radial migration
Abundance distributions of different elements play a fundamental role in determining the evolution of galaxies.Here we used only the iron abundance and the α enhancement in a simple way to separate the thin and thick disk and derive the AMRs.Many other elements are available to test the local enrichment model without radial migration.If radial migration is relevant, then the interpretation of the AMR and the other input functions is more complex as we explained above.The unique AMR, which we derived, corresponds to the mean metallicity as a function of age, and the underlying real AMR is a two-dimensional distribution as derived in many publications based on individual stellar ages (e.g., Bergemann et al. 2014).Our findings for the AMR are in general comparable to the picture presented in Feuillet et al. (2019) where the disk AMR was reconstructed at different R and z based on the Gaia DR2 and APOGEE data.They also see a significant variation of the AMR across the disk with the outer disk being more metal-poor.Their spatial age distribution gives evidence for disk flaring, which we also reproduce in our model.
Local mixing (blurring) produced by the non-circularity of the orbits can be taken into account by replacing the current radius by the guiding radius of the orbit.For the population properties at a given current radius a distribution of guiding radii has to be derived.In a similar way, radial migration working over larger distances via resonances (churning) can be modelled by replacing the guiding radii by distributions of birth radii, which define the efficiency of the migration process.In both cases, the main challenge is to take into account the correlations between different quantities due to the radial gradients in ages, densities, and metallicities.
In order to explain the large spread of abundances in the [α/Fe]-[Fe/H] plane, local enrichment models can be combined with radial migration (Schönrich & Binney 2009;Kubryk et al. 2015;Hayden et al. 2015;Johnson et al. 2021).Minchev et al. (2018) used N-body simulation model for a MW-like galaxy extended by a local chemical enrichment model to derive the evolution and enrichment history of the stellar disk based on a simple radial migration model.Recently, Xiang & Rix (2022) investigated the early MW evolution based on the ages and abundances of about 250 000 subgiant stars.They find a well-mixed α-enhanced thick disk with a similar early formation and enrichment as used in our model, but with an even more extended formation history of 4 Gyr and a higher maximum metallicity.Our thin-disk model corresponds to their dominating low-α evolution sequence.For a more detailed comparison the spatial distribution as a function of galactocentric radii needs to be analysed.

Conclusions
We have presented a publicly available update of our semianalytic MW disk model, which is based on an iterative solving of the Poisson equation and recovers a self-consistent potentialdensity pair independently at each radius14 .We generalised the local JJ model previously calibrated against the Gaia DR2 sample (Paper IV) by introducing radial variation to such input functions as the thin-disk SFR, AVR, and AMR.Our test model, which implies inside-out disk growth, is able to mimic known trends in the MW disk structure, such as broken exponential radial profiles and flaring of the mono-age sub-populations.
As a first step towards the global calibration of the JJ model, we reconstructed the AMR using the extended APOGEE RC sample.We find that the shape of the thin-disk AMR changes a lot within the considered range of galactocentric distances, 4-14 kpc, so we introduced a new analytic function that can reproduce this variation, a seven-parameter double-tanh law.
Applicable to a large volume, 4 kpc R 14 kpc and |z| 2 kpc, the model can test different SFR shapes, disk flaring, IMFs, and AMRs.The JJ model is characterised by a fine time resolution of 25 Myr which can also be combined with a fine spatial resolution within the allowed range of R and z.Thus, the JJ model complemented by the PARSEC, MIST, and BaSTI isochrones can be used for recovering disk evolution, creating MW mock catalogues, or revealing the spiral structure and any other density deviations from the smooth axisymmetric disk.

Fig. 1 .
Fig. 1.Radial density profiles of the JJ-model components.The dashed grey line marks the solar position.

Fig. 2 .
Fig. 2. Thin-disk SFR (top) and AVR (bottom) as functions of radius.Dashed black curves correspond to the solar neighbourhood.Crosses give the vertical velocity dispersion corresponding to the fractions of sub-populations associated with the extra SFR peaks.For the details of this model realisation, see Sect.3.1.
) in Paper IV, there are four parameters of the local AMR: [Fe/H] d0 , [Fe/H] dp , r d , and q.Power-law slopes k [Fe/H] d0 , k [Fe/H] dp , k r d , and k q then prescribe the radial change of the thindisk AMR.When the thin-disk AMR is given by Eqs.(22)

Fig. 3 .
Fig. 3. Radial and vertical structure of the thin disk according to the JJ model realisation TM1 (Sect.3.1).Top: vertical and radial density profiles of the thin-disk mono-age sub-populations.The vertical density profiles correspond to the solar radius.Bottom: radial variation in scale heights of the thin-disk mono-age and mono-abundance sub-populations.

Fig. 4 .
Fig. 4. Normalised W-velocity distribution functions of the different mock populations at the solar radius.

Fig. 5 .
Fig. 5. Selection of the APOGEE RC samples.Top: low-α (green) and high-α (orange) RC sequences.The separation border corresponding to Eq. (19) is shown with dashed black line.Bottom: selected low-α and high-α RC populations in the R-z plane.Dashed lines mark the adopted limits for the galactocentric distance and height.The radial binning of the low-α population corresponds to run0 (see text).The adopted position of the Sun is marked with a yellow dot.

Fig. 6 .
Fig. 6.CMDFs of the APOGEE RC samples.Top: low-α CMDFs calculated in the different radial bins (run0 grid).Dotted coloured curves are CMDFs with extrapolated linear parts.Dashed black curves are C([Fe/H]), convolutions of the 'deconvolved' CMDFs with the Gaussian kernel (Eq.(20)).The kernel dispersions σ [Fe/H] are shown on the plot with the same colour coding.Bottom: same as the top panel, but for the high-α stars and no radial binning.

Fig. 12 .
Fig.12.Same as in Fig.11, but for the thick disk.Both observed and predicted MDs correspond to the range of galactocentric distances 4-14 kpc.

Table 1 .
List of abbreviations used in this paper.

Table 2 .
Parameters of the test generalised JJ model TM1.The thin-and thick-disk AMR parameters are listed separately in Table3.