Nowcasting aerosol particle size distributions from PM2.5 with physics-regularized flow matching

# Nowcasting aerosol particle size distributions from PM2.5 with physics-regularized flow matching
Source: [https://arxiv.org/html/2610.33557](https://arxiv.org/html/2610.33557)
\\Author

\[1,2,3\]\[zhisong\.liu@lut\.fi\]Zhi\-SongLiu\\Author\[1,2,3\]JenniKöykkä\\Author\[2,3\]WenqingPeng\\Author\[1,2,3\]AmmarKheder\\Author\[2,3\]PetriClusius\\Author\[1,2,3\]MichaelBoy

1\]Department of Computational Engineering, LUT University, Finland 2\]Atmospheric Modelling Center, Lahti \(AMC\-Lahti\), Finland 3\]Institute for Atmospheric and Earth System Research \(INAR\), University of Helsinki, Finland

\\pubdiscuss\\published

###### Abstract

The mass concentration of particles smaller than 2\.5μ\\mum \(PM2\.5\) is cheap to measure and widely monitored, whereas the particle number size distribution \(PSD\) requires costly instruments and is measured at only a few research stations\. Yet many aerosol effects, such as vapor uptake, cloud droplet formation and deposition in the human respiratory tract, depend on the PSD rather than on mass\. We explore whether recent generative deep learning models can nowcast the PSD: from the PM2\.5record and eight meteorological parameters over a 24 h window, a model generates an ensemble of plausible hourly size distributions for the same window, reflecting that one mass can arise from many size distributions\. Rather than proposing a single best model, we use two contrasting Finnish stations, the boreal forest site Hyytiälä and the urban site Kumpula, as a controlled testbed\. On it, we benchmark six generative models against deterministic and climatological baselines, locally and after transfer between sites, and test physical constraints for neural network training\. Flow matching performs best overall; its guided variant has the lowest continuous ranked probability score \(CRPS\) at both sites\. A constraint based on a reduced aerosol general dynamic equation \(GDE\) lowers the CRPS by 5\.5–11\.5 % at the training site and reduces errors after transfer between sites\. Transfer remains harder than local estimation, identifying training on large multi\-station data sets as the next step towards PSD estimates at unmonitored locations\. This benchmark establishes physics\-informed generative models as a promising route from routine air\-quality monitoring to size\-resolved aerosol information\.

††firstpage:1\\introduction

Particle number and particle mass describe different aspects of an aerosol population\. Newly formed or freshly emitted ultrafine particles can dominate number concentration while contributing negligibly to mass, whereas far fewer accumulation\-mode particles can dominate PM2\.5, the mass concentration of particles smaller than 2\.5μ\\mum in aerodynamic diameter\. The particle number size distribution \(PSD\) is the state on which aerosol processes act: condensational growth, coagulation and deposition change it, and its shape in turn determines how fast vapors condense onto the particles, how many particles can activate as cloud droplets and where they deposit in the respiratory tract\([Seinfeld and Pandis, 2016](https://arxiv.org/html/2610.33557#bib.bib23);[Kulmala et al\., 2001](https://arxiv.org/html/2610.33557#bib.bib15)\)\. Air\-quality networks routinely report PM2\.5, whereas size\-resolved measurements require differential mobility particle sizers \(DMPS\) with specialist operation and quality control\([Kulmala et al\., 2012](https://arxiv.org/html/2610.33557#bib.bib16);[Wiedensohler et al\., 2012](https://arxiv.org/html/2610.33557#bib.bib27)\)\. The practical question is therefore how much of the PSD can be recovered from routine PM2\.5and meteorological observations\.

Two families of methods approach this question from opposite directions\. Process\-based models such as SOSAA\([Boy et al\., 2011](https://arxiv.org/html/2610.33557#bib.bib3);[Clusius et al\., 2026](https://arxiv.org/html/2610.33557#bib.bib5)\)solve the*forward*problem: given emissions, meteorology and chemistry, they simulate how the PSD evolves by resolving gas\-phase chemistry, boundary\-layer dynamics and aerosol microphysics\. Their inputs, however, include emissions, vertical meteorology and gas\-phase concentrations that are not normally available at a monitoring station\. Recovering the PSD from PM2\.5is instead an*inverse*problem: the size distribution must be inferred from a single mass integral of it, together with the weather\. Machine\-learning methods suit this inverse problem since they learn the mapping from bulk quantities to size information directly from data\.[Vartiainen et al\. \(2025\)](https://arxiv.org/html/2610.33557#bib.bib25)downscaled climate\-model aerosol output to daily number concentrations in three size ranges, and[Barahona et al\. \(2026\)](https://arxiv.org/html/2610.33557#bib.bib2)inferred modal size distributions and mixing state from bulk aerosol mass and meteorology\. Together, these works establish that bulk properties contain recoverable size information\.

However, both studies return a single deterministic estimate, which cannot represent the central difficulty of this inverse problem: the same PM2\.5can arise from a few large particles or many small ones\. We refer to this as the*mass–number ambiguity*\. A regression trained on squared error returns the conditional mean of the plausible spectra, which smooths the hour\-to\-hour and size\-resolved structure that matters for number\-sensitive applications\. In our data, a gradient\-boosting regression reproduces only 58 % of the observed between\-window variability of the PSD, against 83 % for the generative model introduced below \(Appendix[D\.2](https://arxiv.org/html/2610.33557#A4.SS2)\)\. Two further questions are rarely tested\. First, a good full\-spectrum score need not imply accurate derived quantities, because different quantities weight particle diameter differently\. The number of ultrafine particles \(smaller than 100 nm\) is dominated by the smallest particles and is relevant for health and new particle formation\([World Health Organization, 2021](https://arxiv.org/html/2610.33557#bib.bib28)\), whereas the condensation sink, the rate at which condensable vapors are lost to existing particles, is weighted towards larger particles\([Kulmala et al\., 2001](https://arxiv.org/html/2610.33557#bib.bib15)\)\. A model can reproduce one of these quantities and miss the other\. Second, a model fitted at one station may not transfer to a station with a different aerosol regime\.

Conditional generative models\([Lipman et al\., 2023](https://arxiv.org/html/2610.33557#bib.bib17);[Liu et al\., 2023](https://arxiv.org/html/2610.33557#bib.bib18);[Rombach et al\., 2021](https://arxiv.org/html/2610.33557#bib.bib22);[Lovelace et al\., 2023](https://arxiv.org/html/2610.33557#bib.bib20);[Liu et al\., 2025](https://arxiv.org/html/2610.33557#bib.bib19)\)learn a distribution of outputs for a fixed input, and therefore offer a direct way to represent this ambiguity\. Among them, flow matching \(FM\) learns a smooth transport from Gaussian noise to the data by regressing a velocity field, which makes training a simple regression problem and sampling a fixed number of integration steps\([Lipman et al\., 2023](https://arxiv.org/html/2610.33557#bib.bib17);[Liu et al\., 2023](https://arxiv.org/html/2610.33557#bib.bib18)\)\. In this work, we use FM to generate an ensemble of hourly PSDs for a 24 h window from the concurrent PM2\.5and meteorological record \(Fig\.[1](https://arxiv.org/html/2610.33557#S0.F1)\)\. We call this task*PSD nowcasting*: the model reconstructs the PSD of the observed window based on the observed PM2\.5, and receives no PSD measurement at any time\. Such a nowcast could fill gaps in existing size\-distribution records or provide an estimated PSD at stations that measure only PM2\.5\. This also changes the spatial scale of the available information: a size\-resolved measurement describes the aerosol at a single research station, whereas a computational estimate can be produced wherever PM2\.5is monitored, for example at every station of a municipal air\-quality network, and thus at a much higher spatial density\. Since a generative objective alone does not enforce aerosol physics, we add soft penalties of two kinds:*integral constraints*on resolved mass and number, and a*dynamical constraint*based on a reduced general dynamic equation \(GDE\)\.

Since generative models have not yet been applied to this problem, our aim is exploratory\. Rather than proposing a single best model, we benchmark the main families of generative models for PSD nowcasting and test whether physical constraints can guide neural network training\. A benchmark of this kind needs controlled conditions more than scale: every model must see the same inputs, training period and evaluation windows, so that differences can be attributed to the modelling choices\. We therefore use two well\-characterized stations with long, continuous PSD records as a testbed: the Hyytiälä boreal forest station and the Kumpula urban station in Finland\. Their contrasting aerosol regimes also make the pair a deliberately demanding test of transfer between environments\. Within this testbed we address four research questions:

- •RQ1\.Does conditional generation reconstruct the full spectrum better than deterministic, climatological, and alternative generative methods?
- •RQ2\.Do physical constraints improve the generated spectra, and does the form of the constraint matter?
- •RQ3\.For two derived quantities with contrasting size weighting, ultrafine particle number and the condensation sink, are the generated spectra more accurate than past PSDs observed in the same season and at the same hour of day?
- •RQ4\.How much of this information survives transfer between the boreal and the urban site?

Our contributions are a benchmark and a set of diagnostics rather than a single optimized model, and they map onto the three gaps above\. First, we formulate PSD nowcasting as conditional generation, so that the mass–number ambiguity appears as ensemble spread rather than being averaged away, and compare six generative models on equal terms\. Second, we compare integral and dynamical physical constraints under one training objective\. Third, we evaluate the spectra on derived quantities and on transfer between sites\. As the reference, we sample complete PSD sequences observed in the training period in the same season and at the same hour of day; this historical ensemble shows what could be estimated without any knowledge of the current PM2\.5and weather, so improving on it shows what those observations add\. Our main findings are: \(1\) all six generative models have a lower continuous ranked probability score \(CRPS\) than gradient boosting at both sites; \(2\) the dynamical constraint lowers full\-spectrum CRPS by 11\.5 % and 5\.5 % at the two sites, whereas either integral constraint alone raises it; \(3\) the GDE\-constrained model has the lowest CRPS of the four learned models for both derived quantities at both sites; and \(4\) transferability depends on both the quantity and the direction of transfer, and boreal\-to\-urban transfer fails for ultrafine particle number\. These findings are a proof of concept on a controlled two\-station testbed\. They show which modelling choices matter and where models trained on a single station fail, which is the information needed before scaling the approach to large, multi\-station training data\. Section[1](https://arxiv.org/html/2610.33557#S1)describes the observations and the task, Sect\.[2](https://arxiv.org/html/2610.33557#S2)the model and evaluation, and Sect\.[3](https://arxiv.org/html/2610.33557#S3)the results, organized by research question\. Section[4](https://arxiv.org/html/2610.33557#S4)interprets the results, lists limitations and outlines future work\.

![Refer to caption](https://arxiv.org/html/2610.33557v1/figures/teaser.jpg)Figure 1:PSD nowcasting as conditional generation\.The same 24 h of PM2\.5and meteorology that condition the transformation of Gaussian noise also cover the generated hourly PSD: the target window is concurrent with the input record\. Repeated noise draws give an ensemble for a fixed input\. Generative integration time is distinct from the physical \(clock\) time of the measurements\.## 1Observations and nowcasting task

### 1\.1A boreal–urban station pair

We use SMEAR II at Hyytiälä \(61\.85∘N, 24\.29∘E\) and SMEAR III at Kumpula, Helsinki \(60\.20∘N, 24\.96∘E\)\. Hyytiälä is a boreal forest station where secondary particle formation and growth are prominent features of the particle population\([Hari and Kulmala, 2005](https://arxiv.org/html/2610.33557#bib.bib9);[Dal Maso et al\., 2005](https://arxiv.org/html/2610.33557#bib.bib6)\)\. Kumpula is an urban background station near a major road, with local traffic contributions alongside regional aerosol\([Järvi et al\., 2009](https://arxiv.org/html/2610.33557#bib.bib14)\)\. The two stations therefore represent contrasting situations in which emissions and meteorology shape the observed particle mass and number concentrations\. We chose them as a testbed, not as a representative sample of aerosol environments: their seven\-year common record allows identical training and evaluation periods for all models, and their contrast provides an out\-of\-distribution test in each direction\.

### 1\.2Data

The data set used in this study covers 1 January 2018–31 December 2024 at both stations, at hourly resolution\. DMPS observations are averaged from native 10 min measurements\. We retain the 50 archived diameter bins between 2\.82 and 794 nm, without diameter interpolation\. This grid resolves part of the nucleation range and the Aitken and accumulation ranges, and its limits define the support of every predicted PSD and derived quantity\.

The upper limit of 794 nm has a direct consequence for mass\. PM2\.5includes particles up to 2\.5μ\\mum in aerodynamic diameter, which corresponds to a mobility diameter of about 2μ\\mum for a density of 1\.5 g cm\-3\(Appendix[C](https://arxiv.org/html/2610.33557#A3)\)\. Even the full DMPS range, which ends at 1000 nm, does not reach this size\. The mass integrated from a predicted spectrum is therefore a resolved, partial mass that cannot be required to equal PM2\.5; Sect\.[2\.2](https://arxiv.org/html/2610.33557#S2.SS2)explains how we use it only as a soft reference\.

The model inputs are in\-situ PM2\.5and eight ERA5 variables from the nearest grid cell: 2 m temperature and dewpoint, 10 m wind components, surface pressure, boundary\-layer height, total precipitation and downward solar radiation\([Hersbach et al\., 2020](https://arxiv.org/html/2610.33557#bib.bib10)\)\. Additionally, a condensation particle counter \(CPC\) is used for one training constraint\. PSDs are log\-transformed and normalized per bin\. The transformations, data checks and the unverified fitting period of the normalization are in Appendix[A](https://arxiv.org/html/2610.33557#A1)\.

### 1\.3Nowcasting task

The task is a nowcast, not a forecast: the model estimates the PSD for the same hours as its inputs\. For a window ending at hourτ\\tau, the model receives the PM2\.5and ERA5 record of those 24 hours,𝐜τ=\(PM2\.5,ERA5\)τ−23:τ\\mathbf\{c\}\_\{\\tau\}=\(\\mathrm\{PM\}\_\{2\.5\},\\mathrm\{ERA5\}\)\_\{\\tau\-23:\\tau\}, and generates the PSD for the same 24 hours,𝐱τ=PSDτ−23:τ∈ℝ50×24\\mathbf\{x\}\_\{\\tau\}=\\mathrm\{PSD\}\_\{\\tau\-23:\\tau\}\\in\\mathbb\{R\}^\{50\\times 24\}\(50 diameter bins by 24 hours\)\. No PSD measurement is given to the model at any time\. A 24 h window covers one full diurnal cycle of boundary\-layer development, radiation and emissions, which strongly shape the PSD\. Each example therefore requires 24 consecutive valid hours\. Models are trained separately at each site and evaluated on a later period at that site and, without retraining, at the other site\.

## 2Methods

### 2\.1Conditional generation with flow matching

Flow Matching \(FM\)\([Lipman et al\., 2023](https://arxiv.org/html/2610.33557#bib.bib17)\)learns a velocity field𝐯θ\\mathbf\{v\}\_\{\\theta\}that transports Gaussian noise𝐱0\\mathbf\{x\}\_\{0\}to a PSD sequence𝐱1\\mathbf\{x\}\_\{1\}, conditional on𝐜\\mathbf\{c\}\. During training,ttis sampled uniformly on\[0,1\]\[0,1\], the intermediate state is𝐱t=\(1−t\)​𝐱0\+t​𝐱1\\mathbf\{x\}\_\{t\}=\(1\-t\)\\mathbf\{x\}\_\{0\}\+t\\mathbf\{x\}\_\{1\}, and the target velocity is𝐯∗=𝐱1−𝐱0\\mathbf\{v\}^\{\*\}=\\mathbf\{x\}\_\{1\}\-\\mathbf\{x\}\_\{0\}:

ℒFM=𝔼t,𝐱0,\(𝐱1,𝐜\)​‖𝐯θ​\(𝐱t,t,𝐜\)−𝐯∗‖2\.\\mathcal\{L\}\_\{\\mathrm\{FM\}\}=\\mathbb\{E\}\_\{t,\\mathbf\{x\}\_\{0\},\(\\mathbf\{x\}\_\{1\},\\mathbf\{c\}\)\}\\left\\\|\\mathbf\{v\}\_\{\\theta\}\(\\mathbf\{x\}\_\{t\},t,\\mathbf\{c\}\)\-\\mathbf\{v\}^\{\\ast\}\\right\\\|^\{2\}\.\(1\)Herettis a generative time that indexes the integration from noise to data; it is distinct from the physical \(clock\) timeτ\\tauof the measurements\. At inference, 50 Euler steps integrate the learned field from a fresh Gaussian draw\. Different draws give different PSD sequences for the same input, and eight draws form the ensemble\. The network𝐯θ\\mathbf\{v\}\_\{\\theta\}is a Transformer\([Vaswani et al\., 2017](https://arxiv.org/html/2610.33557#bib.bib26)\)that treats the 24 hourly spectra of a window as a sequence: each hour is embedded as one element, and self\-attention lets the estimate for every hour draw on the inputs and spectra of all other hours\. We chose this architecture because the PSD evolves through the day, so the spectrum at one hour is informative about its neighbours\. The network has approximately1\.24×1061\.24\\times 10^\{6\}parameters \(Appendix[B](https://arxiv.org/html/2610.33557#A2)\), and Fig\.[8](https://arxiv.org/html/2610.33557#A2.F8)illustrates the generation path\.

### 2\.2Physical constraints

A generative objective does not by itself make the generated spectra obey aerosol physics\. We therefore add soft physical penalties to the training objective\. Physical relationships are defined for clean spectra, not for the partly noised intermediate states𝐱t\\mathbf\{x\}\_\{t\}, so we evaluate them on the one\-step clean\-sample estimate𝐱^1=𝐱t\+\(1−t\)​𝐯θ\\hat\{\\mathbf\{x\}\}\_\{1\}=\\mathbf\{x\}\_\{t\}\+\(1\-t\)\\mathbf\{v\}\_\{\\theta\}, obtained by inverting the interpolation, after conversion to physical units\. This estimate reuses the network output, so the penalties require no additional network evaluation\. The total objective is

ℒ=ℒFM\+λm​t2​ℒmass\+λn​t2​ℒCN\+λg​t3​ℒGDE\.\\mathcal\{L\}=\\mathcal\{L\}\_\{\\mathrm\{FM\}\}\+\\lambda\_\{\\mathrm\{m\}\}t^\{2\}\\mathcal\{L\}\_\{\\mathrm\{mass\}\}\+\\lambda\_\{\\mathrm\{n\}\}t^\{2\}\\mathcal\{L\}\_\{\\mathrm\{CN\}\}\+\\lambda\_\{\\mathrm\{g\}\}t^\{3\}\\mathcal\{L\}\_\{\\mathrm\{GDE\}\}\.\(2\)Near the noise end of the generative path \(t→0t\\to 0\),𝐱^1\\hat\{\\mathbf\{x\}\}\_\{1\}is little more than a guess, and penalizing it would mainly add noise to training\. The factorst2t^\{2\}andt3t^\{3\}therefore down\-weight the penalties there\. The GDE term receives the steeper weight because it involves derivatives in time and size, which amplify any remaining noise\. The two integral constraints \(ℒmass\\mathcal\{L\}\_\{\\mathrm\{mass\}\}andℒCN\\mathcal\{L\}\_\{\\mathrm\{CN\}\}\) restrict spectral integrals; the dynamical constraint \(ℒGDE\\mathcal\{L\}\_\{\\mathrm\{GDE\}\}\) links neighbouring sizes and hours\. Coefficients and numerical details are given in Appendix[C](https://arxiv.org/html/2610.33557#A3)\.

#### Integral constraint 1: pseudo\-mass consistency

Letn^h,i=d​N/d​log10​Dp,i\\hat\{n\}\_\{h,i\}=\\mathrm\{d\}N/\\mathrm\{d\}\\log\_\{10\}D\_\{p,i\}be the estimated number concentration in diameter biniiat hourhh, normalized by the logarithmic bin widthΔi\\Delta\_\{i\}, so thatn^h,i​Δi\\hat\{n\}\_\{h,i\}\\Delta\_\{i\}is the number concentration of particles in that bin\. Assuming spherical particles of densityρ=1\.5\\rho=1\.5g cm\-3, the mass resolved by the spectrum is

M^res,h=ρ​π6​∑i=1BDp,i3​n^h,i​Δi,B=50,\\hat\{M\}\_\{\\mathrm\{res\},h\}=\\frac\{\\rho\\pi\}\{6\}\\sum\_\{i=1\}^\{B\}D\_\{p,i\}^\{3\}\\hat\{n\}\_\{h,i\}\\Delta\_\{i\},\\qquad B=50,\(3\)which is compared with the observed PM2\.5,mhobsm\_\{h\}^\{\\rm obs\}, through a robust log\-space discrepancy:

ℒmass=1T​∑h=1THuber0\.5⁡\[ln⁡\(1\+M^res,h/M0\)−ln⁡\(1\+max⁡\(mhobs,0\)/M0\)\],M0=1​μ​g​m−3,\\mathcal\{L\}\_\{\\rm mass\}=\\frac\{1\}\{T\}\\sum\_\{h=1\}^\{T\}\\operatorname\{Huber\}\_\{0\.5\}\\\!\\left\[\\ln\(1\+\\hat\{M\}\_\{\\mathrm\{res\},h\}/M\_\{0\}\)\-\\ln\(1\+\\max\(m\_\{h\}^\{\\rm obs\},0\)/M\_\{0\}\)\\right\],\\quad M\_\{0\}=1\\,\\mu\\mathrm\{g\}\\,\\mathrm\{m\}^\{\-3\},\(4\)withT=24T=24hours\. The predicted spectrum covers only part of the size range sampled by the PM2\.5instrument, so exact agreement could only be reached by adding spurious mass within the resolved range\. We therefore compare the two masses in logarithmic space, where a constant support mismatch becomes a constant offset, and use the Huber function, which is quadratic for small discrepancies and linear for large ones, so that this systematic mismatch cannot dominate the loss\. Together with a finite weight, this imposes a soft reference rather than an equality\. We call this term*pseudo\-mass consistency*and label it “mass” in the tables\. In the nowcast, the target window coincides with the input window, somhobsm\_\{h\}^\{\\rm obs\}is the same PM2\.5record that the model receives as input\. The term is thus a self\-consistency penalty between input and output\.

#### Integral constraint 2: aerosol\-number consistency

Mass constrains small particles only weakly, whereas total number is dominated by them\. An independent measurement of total number can therefore constrain the part of the spectrum that mass resolves poorly\. Letchobsc\_\{h\}^\{\\rm obs\}be the total number concentration measured by the CPC at hourhh, and letℋCN⊆\{1,…,T\}\\mathcal\{H\}\_\{\\rm CN\}\\subseteq\\\{1,\\ldots,T\\\}be the set of hours in the window with a valid CPC measurement\. We compare the resolved numberCN^res,h\\hat\{\\mathrm\{CN\}\}\_\{\\mathrm\{res\},h\}with the CPC record:

CN^res,h=∑i=1Bn^h,i​Δi,ℒCN=1\|ℋCN\|​∑h∈ℋCNHuber0\.5⁡\[ln⁡\(1\+CN^res,h/N0\)−ln⁡\(1\+max⁡\(chobs,0\)/N0\)\],N0=1​cm−3\.\\hat\{\\mathrm\{CN\}\}\_\{\\mathrm\{res\},h\}=\\sum\_\{i=1\}^\{B\}\\hat\{n\}\_\{h,i\}\\Delta\_\{i\},\\qquad\\mathcal\{L\}\_\{\\rm CN\}=\\frac\{1\}\{\|\\mathcal\{H\}\_\{\\rm CN\}\|\}\\sum\_\{h\\in\\mathcal\{H\}\_\{\\rm CN\}\}\\operatorname\{Huber\}\_\{0\.5\}\\\!\\left\[\\ln\(1\+\\hat\{\\mathrm\{CN\}\}\_\{\\mathrm\{res\},h\}/N\_\{0\}\)\-\\ln\(1\+\\max\(c\_\{h\}^\{\\rm obs\},0\)/N\_\{0\}\)\\right\],\\quad N\_\{0\}=1\\,\\mathrm\{cm\}^\{\-3\}\.\(5\)The CPC and DMPS have different detection ranges and responses, so their totals are closely related but not identical \(log\-space correlations of 0\.988 and 0\.963, median DMPS/CPC ratios 0\.87 and 1\.03 at the two sites; Appendix[C](https://arxiv.org/html/2610.33557#A3)\)\.

#### Dynamical constraint: reduced general dynamic equation \(GDE\)

Aerosol dynamical processes, namely new particle formation, condensational growth, coagulation and deposition, change the particle population over time\. Growth and coagulation in particular link neighbouring sizes and successive hours, a structure that the integral constraints cannot express\. Atmospheric mixing is not an aerosol\-dynamical process, but it also changes the population observed at a fixed station, so we include boundary\-layer dilution as an additional term\. We penalize the relative residual of a reduced GDE,

rh,i=1n^h,i\+ϵ​\[∂n^∂τ\+∂∂ln⁡Dp​\(GR^Dp​n^\)−S^\+CoagS^​n^\+λdil​n^\]h,i,r\_\{h,i\}=\\frac\{1\}\{\\hat\{n\}\_\{h,i\}\+\\epsilon\}\\left\[\\frac\{\\partial\\hat\{n\}\}\{\\partial\\tau\}\+\\frac\{\\partial\}\{\\partial\\ln D\_\{p\}\}\\left\(\\frac\{\\widehat\{\\mathrm\{GR\}\}\}\{D\_\{p\}\}\\hat\{n\}\\right\)\-\\hat\{S\}\+\\widehat\{\\mathrm\{CoagS\}\}\\,\\hat\{n\}\+\\lambda\_\{\\mathrm\{dil\}\}\\hat\{n\}\\right\]\_\{h,i\},\(6\)whereτ\\tauis physical \(clock\) time in hours andϵ\\epsilonstabilizes the division\. The growth rateGR^\\widehat\{\\mathrm\{GR\}\}moves particles through size space, the coagulation sinkCoagS^\\widehat\{\\mathrm\{CoagS\}\}is evaluated on the predicted PSD,S^\\hat\{S\}absorbs unresolved sources and transport, andλdil=\(d​H/d​τ\)\+/H\\lambda\_\{\\mathrm\{dil\}\}=\(\\mathrm\{d\}H/\\mathrm\{d\}\\tau\)\_\{\+\}/His the dilution rate from boundary\-layer heightHH\. An auxiliary head reads only the conditioning and predicts growth and source terms together with a bounded log residual variancesh,is\_\{h,i\}\. The variance is needed because a single station cannot close the aerosol budget: on some days, advection of a different air mass dominates every local term, and the model should be able to down\-weight such hours rather than be forced to balance them\. The penalty is therefore the heteroscedastic Gaussian negative log\-likelihood

ℒGDE=1T​B​∑h=1T∑i=1B\[exp⁡\(−sh,i\)​rh,i2\+sh,i\],\\mathcal\{L\}\_\{\\mathrm\{GDE\}\}=\\frac\{1\}\{TB\}\\sum\_\{h=1\}^\{T\}\\sum\_\{i=1\}^\{B\}\\left\[\\exp\(\-s\_\{h,i\}\)r\_\{h,i\}^\{2\}\+s\_\{h,i\}\\right\],\(7\)in which the second term prevents the model from lowering the loss merely by inflating the variance\. The equation is reduced: new particle formation enters only as a growth flux into the smallest resolved bin and throughS^\\hat\{S\}, deposition and the coagulation gain term are not represented explicitly, andS^\\hat\{S\}can absorb these and other omitted processes\. It therefore regularizes the generated evolution, but it does not close the aerosol budget at an open measurement site\.

### 2\.3Comparison models and references

Generative models\.The five comparison models span the main ways of learning a conditional distribution: iterative denoising, alternative transport couplings, guidance, one\-step generation and energy\-based modelling\. A denoising diffusion model\([Ho et al\., 2020](https://arxiv.org/html/2610.33557#bib.bib12)\)reuses the FM backbone as a noise predictor\. FM\+OT replaces independent noise–data pairing with minibatch optimal\-transport coupling\([Tong et al\., 2024](https://arxiv.org/html/2610.33557#bib.bib24)\)\. FM\+CFG uses classifier\-free guidance\([Ho and Salimans, 2021](https://arxiv.org/html/2610.33557#bib.bib11)\), with the tabulated operating point at guidance scalew=1w=1\. The Drifting Model\([Deng et al\., 2026](https://arxiv.org/html/2610.33557#bib.bib7)\)generates with one network evaluation, and Energy Matching\([Balcerak et al\., 2025](https://arxiv.org/html/2610.33557#bib.bib1)\)learns a scalar potential\.

Deterministic models\.Gradient boosting fits one histogram regressor per diameter bin, and a multilayer perceptron \(MLP\) predicts all24×5024\\times 50values jointly\. Both return one deterministic PSD per input, and so show what is gained or lost when the ambiguity is not represented\.

Non\-learning references\.A month/hour*climatology*is the reference for the full spectrum\. For the derived quantities, we use a deterministic*seasonal template*\(mean PSD by season and UTC hour\) and a probabilistic historical ensemble\. Each member of the historical ensemble is a complete 24 h PSD sequence from the training period, matched to the target season and starting hour, so that it retains realistic dependence across bins and hours\. It uses no PM2\.5or meteorology and therefore measures directly what the conditioning adds\. All references are built from source\-site training data only\.

### 2\.4Evaluation design and fairness

Training ends at 2 November 2024, 00:00 UTC\. The temporal holdout covers 2 November–31 December 2024 and contains 1417 hourly\-stride windows per site\. Table[1](https://arxiv.org/html/2610.33557#S2.T1)lists the three evaluation populations\. In the paired transfer test, a model trained at the other station, applied without retraining, and a locally trained model are scored on identical target\-site windows, so their ratio isolates the cost of transfer\. For a fair comparison, all learned models share the same inputs, target grid, training period and evaluation windows\. Since the purpose is a controlled comparison of model families and constraints, we keep these settings fixed rather than tuning each model to its best achievable score on the test data; the reported scores therefore indicate relative performance under identical conditions, not the upper limit of any method\.

Table 1:Evaluation populations\. H and K denote Hyytiälä and Kumpula\. All predictions cover the same 24 hours as their input record\.\\tophlineTestScoring populationQuestion\\middlehlineSame\-site temporal holdout1417 windows at each site; 2 November–31 December 2024Local skill on a later period \(RQ1–RQ3\)Paired transferH\-trained and K\-trained models on the identical target\-site holdoutCost of transfer for each quantity \(RQ4\)Extended full\-spectrum transferAll valid windows at the other siteRobustness over a longer period \(Appendix[E](https://arxiv.org/html/2610.33557#A5)\)\\bottomhline
### 2\.5Metrics

Table[2](https://arxiv.org/html/2610.33557#S2.T2)summarizes the metrics\. Each ensemble member is converted to physical units before scoring\. For evaluated bin–hour pairsj=\(w,h,i\)j=\(w,h,i\), indexed by windowww, the ensemble mean overM=8M=8members isn¯j\\bar\{n\}\_\{j\}, and the root mean square error overQQpairs is

n¯j=1M​∑m=1Mn^j\(m\),RMSE=\[1Q​∑j\(n¯j−njobs\)2\]1/2\.\\bar\{n\}\_\{j\}=\\frac\{1\}\{M\}\\sum\_\{m=1\}^\{M\}\\hat\{n\}\_\{j\}^\{\(m\)\},\\qquad\\mathrm\{RMSE\}=\\left\[\\frac\{1\}\{Q\}\\sum\_\{j\}\(\\bar\{n\}\_\{j\}\-n\_\{j\}^\{\\rm obs\}\)^\{2\}\\right\]^\{1/2\}\.\(8\)The CRPS\([Gneiting and Raftery, 2007](https://arxiv.org/html/2610.33557#bib.bib8)\)scores the whole ensemble; we use its finite\-ensemble \(fair\) form

CRPSj=1M​∑m=1M\|zm−y\|−12​M​\(M−1\)​∑m=1M∑k≠m\|zm−zk\|,\\mathrm\{CRPS\}\_\{j\}=\\frac\{1\}\{M\}\\sum\_\{m=1\}^\{M\}\|z\_\{m\}\-y\|\-\\frac\{1\}\{2M\(M\-1\)\}\\sum\_\{m=1\}^\{M\}\\sum\_\{k\\neq m\}\|z\_\{m\}\-z\_\{k\}\|,\(9\)withy=njobsy=n\_\{j\}^\{\\rm obs\}andzm=n^j\(m\)z\_\{m\}=\\hat\{n\}\_\{j\}^\{\(m\)\}, averaged over all pairs\. The first term penalizes distance from the observation; the second credits spread, so that spread alone cannot improve the score\. For a deterministic prediction, CRPS equals the absolute error\. The skill score relative to the historical ensemble is

CRPSS=1−CRPS¯model/CRPS¯historical\.\\mathrm\{CRPSS\}=1\-\\overline\{\\mathrm\{CRPS\}\}\_\{\\rm model\}/\\overline\{\\mathrm\{CRPS\}\}\_\{\\rm historical\}\.\(10\)Interval reliability uses the empirical 10th and 90th percentilesq0\.1,jq\_\{0\.1,j\}andq0\.9,jq\_\{0\.9,j\}:

C80=1Q∑j𝟏\{q0\.1,j≤yj≤q0\.9,j\},W80=1Q∑j\(q0\.9,j−q0\.1,j\)\.C\_\{80\}=\\frac\{1\}\{Q\}\\sum\_\{j\}\\mathbf\{1\}\\\{q\_\{0\.1,j\}\\leq y\_\{j\}\\leq q\_\{0\.9,j\}\\\},\\qquad W\_\{80\}=\\frac\{1\}\{Q\}\\sum\_\{j\}\(q\_\{0\.9,j\}\-q\_\{0\.1,j\}\)\.\(11\)The indicator function𝟏​\{⋅\}\\mathbf\{1\}\\\{\\cdot\\\}counts an observation when it lies inside the interval\. A valueC80=0\.60C\_\{80\}=0\.60, for example, means that only 60 % of observations fall inside a nominal 80 % interval\. Block\-bootstrap intervals resample 24 contiguous windows \(Appendix[D](https://arxiv.org/html/2610.33557#A4)\)\.

Table 2:Metrics, their interpretation and optimal values\. RMSE and CRPS carry the units of the scored quantity\.\\tophlineMetricWhat it measuresOptimum\\middlehlineRMSE, Eq\. \([8](https://arxiv.org/html/2610.33557#S2.E8)\)Point error of the ensemble mean, weighting large concentration errors0 \(lower is better\)corr\.Pearson correlation of log\-spectra across bins and hours within a window, averaged over windows\.1CRPS, Eq\. \([9](https://arxiv.org/html/2610.33557#S2.E9)\)Accuracy of the whole predictive distribution0 \(lower is better\)CRPSS, Eq\. \([10](https://arxiv.org/html/2610.33557#S2.E10)\)Improvement over the historical ensemble1; 0 equals the referenceC80C\_\{80\}, Eq\. \([11](https://arxiv.org/html/2610.33557#S2.E11)\)Fraction of observations inside the central 80 % interval0\.80W80W\_\{80\}, Eq\. \([11](https://arxiv.org/html/2610.33557#S2.E11)\)Mean width of that interval; read withC80C\_\{80\}Narrowest width atC80=0\.80C\_\{80\}=0\.80Brier scoreMean squared error of an event probability \(Appendix[F\.1](https://arxiv.org/html/2610.33557#A6.SS1)\)0\\bottomhline#### Derived atmospheric quantities

As motivated in the Introduction, we evaluate two quantities with contrasting diameter weighting\. The first is the number concentration of ultrafine particles \(UFP number\) over the observed range 2\.82–100 nm\. Its upper bound follows the conventional definition of ultrafine particles as those smaller than 100 nm\([World Health Organization, 2021](https://arxiv.org/html/2610.33557#bib.bib28)\); its lower bound is set by the DMPS grid\. The second is the resolved condensation sink\([Kulmala et al\., 2001](https://arxiv.org/html/2610.33557#bib.bib15)\),

CS=2​π​D​∑iβm​\(Dp,i\)​Dp,i​Ni,\\mathrm\{CS\}=2\\pi D\\sum\_\{i\}\\beta\_\{m\}\(D\_\{p,i\}\)\\,D\_\{p,i\}\\,N\_\{i\},\(12\)whereNi=ni​Δ​log10​Dp,iN\_\{i\}=n\_\{i\}\\Delta\\log\_\{10\}D\_\{p,i\},DDis the vapor diffusion coefficient andβm\\beta\_\{m\}the Fuchs–Sutugin correction\. We use sulfuric\-acid vapor properties, unit accommodation, 293\.15 K and 101325 Pa for both predicted and observed PSDs, so that errors reflect the size distribution only\. Number counts every particle equally, whereas the sink weights each particle by its diameter and a transition\-regime correction, so that larger particles contribute more\. Each ensemble member is integrated separately before scoring, which preserves the dependence across sizes\.

![Refer to caption](https://arxiv.org/html/2610.33557v1/fig_grid_extrap.png)Figure 2:Generative models reproduce hour\-to\-hour spectral structure that deterministic and climatological estimates smooth out\.Nowcast 24 h PSDs for one temporal\-holdout day at Hyytiälä, compared with the observation \(top left\)\. Each generative panel shows a single ensemble member\. Color islog10⁡d​N/d​log⁡Dp\\log\_\{10\}\\,\\mathrm\{d\}N/\\mathrm\{d\}\\log D\_\{p\}; the vertical axis is mobility diameter and the horizontal axis is UTC hour\.

## 3Results

We first inspect generated spectra and full\-spectrum scores \(RQ1\), then the effect of physical constraints \(RQ2\), the derived quantities \(RQ3\) and transfer \(RQ4\)\. Finally, we place the spectra in the context of a process model\. Throughout, the aim is to compare model families and constraints under identical conditions and to identify where they fail\.

### 3\.1RQ1: Generative ensembles give the best distributional skill, while gradient boosting remains competitive on point error

Before scoring, we check whether the generated spectra look like observed spectra\. On a representative holdout day \(Fig\.[2](https://arxiv.org/html/2610.33557#S2.F2)\), the FM variants reproduce the hour\-to\-hour movement of the dominant band and the abrupt change in the late afternoon\. Gradient boosting and climatology return smooth bands that barely change through the day\. The Drifting Model shows hour\-to\-hour striping, and the MLP reproduces the structure but exaggerates the small\-particle band\. Aggregated over the holdout, FM retains 83 % of the observed between\-window standard deviation, whereas gradient boosting retains 58 % \(Appendix[D\.2](https://arxiv.org/html/2610.33557#A4.SS2)\)\.

Figure 3:Generative ensembles have lower CRPS than gradient boosting at both sites, while point errors are similar\.Top: RMSE of the ensemble mean \(deterministic models supply a single prediction\)\. Bottom: CRPS with 95 % block\-bootstrap intervals\. Lower is better\.Table[3](https://arxiv.org/html/2610.33557#S3.T3)and Fig\.[3](https://arxiv.org/html/2610.33557#S3.F3)give the full scores\. All six generative models have lower CRPS than gradient boosting and climatology at both sites\. FM\+CFG has the lowest CRPS among the unconstrained models, 126 cm\-3at Hyytiälä and 421 cm\-3at Kumpula, compared with 142 cm\-3and 453 cm\-3for plain FM and 172 cm\-3and 663 cm\-3for gradient boosting, respectively\.

Point scores tell a more mixed story\. At Hyytiälä, gradient boosting has the lowest RMSE \(354 cm\-3, against 389 cm\-3for FM\) and the highest correlation \(0\.862, against 0\.820\)\. This advantage is expected: gradient boosting minimizes the squared error in each bin and therefore returns the conditional mean, which is exactly what RMSE rewards\. At Kumpula, FM\+OT has the lowest RMSE \(1427 cm\-3\), and five of the six generative models are below gradient boosting \(1571 cm\-3\)\.

For the constraint and transfer experiments that follow, we retain plain FM as the base model\. Guidance and optimal\-transport coupling change the sampling and the training pairs, respectively, independently of the physical penalties, so they could be combined with the constraints later; using the plain model keeps the effect of each constraint isolated\. The guidance trade\-off is analysed in Appendix[D\.4](https://arxiv.org/html/2610.33557#A4.SS4)\.

Table 3:Full\-spectrum nowcast scores on the same\-site temporal holdouts\.corr\. uses log\-transformed spectra; CRPS and RMSE are in cm\-3, with 95 % block\-bootstrap intervals for CRPS\. Best values are in bold\.\(a\) Hyytiälä

\\tophlineModelcorr\.↑\\uparrowCRPS↓\\downarrowRMSE↓\\downarrow\\middlehlineClimatology0\.776234 \(195–274\)442Gradient Boosting0\.862172 \(144–204\)354MLP0\.783246 \(213–281\)513Diffusion Model0\.771132 \(110–158\)382Drifting Model0\.779152 \(127–180\)425Energy Matching0\.821129 \(110–152\)376Flow Matching0\.820142 \(118–171\)389FM \+ OT0\.822146 \(122–175\)383FM \+ CFG0\.826126 \(104–153\)357\\bottomhline
\(b\) Kumpula

\\tophlineModelcorr\.↑\\uparrowCRPS↓\\downarrowRMSE↓\\downarrow\\middlehlineClimatology0\.889782 \(667–906\)1710Gradient Boosting0\.915663 \(575–762\)1571MLP0\.887856 \(774–946\)1993Diffusion Model0\.926511 \(412–628\)1568Drifting Model0\.888535 \(475–606\)2059Energy Matching0\.925437 \(375–508\)1496Flow Matching0\.929453 \(391–529\)1509FM \+ OT0\.932442 \(389–505\)1427FM \+ CFG0\.932421 \(366–488\)1452\\bottomhline

### 3\.2RQ2: A dynamical constraint improves the spectra; integral constraints alone do not

Table[4](https://arxiv.org/html/2610.33557#S3.T4)shows that the form of the constraint matters more than its presence\. The GDE constraint lowers CRPS from 141\.9 cm\-3to 125\.7 cm\-3\(−11\.5\-11\.5%\) at Hyytiälä and from 455\.3 cm\-3to 430\.5 cm\-3\(−5\.5\-5\.5%\) at Kumpula, and it also lowers the PSD RMSE and raises the correlation at both sites\. The integral constraints have the opposite effect when used alone: pseudo\-mass consistency raises CRPS to 157\.8 cm\-3\(\+11\+11%\) at Hyytiälä and 528\.3 cm\-3\(\+16\+16%\) at Kumpula, and number consistency raises it by 20 % at both sites\. Combining them with GDE recovers part of this loss\. Mass\+GDE gives the lowest CRPS and RMSE at Kumpula \(428\.4 cm\-3and 1430\.0 cm\-3, respectively\), whereas GDE alone remains best at Hyytiälä, and number\+GDE does not reach the unconstrained CRPS at either site\.

Within these formulations, therefore, constraining an integral is not enough to recover how the corresponding quantity is distributed across size and time\. A mismatch in measurement support also explains why the integral constraints need care\. The predicted spectrum ends at 794 nm, whereas the PM2\.5instrument also samples larger particles, so forcing agreement between the two can push the spectrum in the wrong direction\. At Hyytiälä, pseudo\-mass consistency alone even raises the mass RMSE from 0\.19μ\\mug m\-3to 0\.22μ\\mug m\-3, while GDE alone lowers it to 0\.18μ\\mug m\-3\. Since the mass reference is the input PM2\.5itself \(Sect\.[2\.2](https://arxiv.org/html/2610.33557#S2.SS2)\), these mass\-RMSE values measure input–output consistency rather than independent accuracy\. Number is better matched: CPC number correlates closely with DMPS\-integrated number, although the median DMPS/CPC ratios of 0\.87 and 1\.03 show that the two instruments still differ \(Sect\.[2\.2](https://arxiv.org/html/2610.33557#S2.SS2)\)\.

Table 4:Physical\-constraint ablation on Flow Matching for the same\-site temporal holdouts\.CRPS and PSD RMSE are in cm\-3; mass RMSE is inμ\\mug m\-3and measures consistency with the input PM2\.5\. Best values are in bold\.\(a\) Hyytiälä

\(b\) Kumpula

### 3\.3RQ3: Routine observations add information beyond historical spectra

A model can improve a bin\-wise score without improving a derived quantity\. We therefore compare ultrafine particle number and the condensation sink with the historical ensemble, which already contains realistic spectra and their seasonal variation \(Fig\.[4](https://arxiv.org/html/2610.33557#S3.F4), Table[3\.3](https://arxiv.org/html/2610.33557#S3.SS3)\)\. Every learned model matches or improves on this reference at both sites \(CRPSS≥0\\geq 0; Table[F](https://arxiv.org/html/2610.33557#A6)\); the weakest case is gradient boosting, which only matches it for ultrafine particle number at Hyytiälä \(CRPSS 0\.00\)\. FM\+GDE has the lowest CRPS of the four learned models in all four site–quantity comparisons\. Its CRPSS is 0\.27 and 0\.42 for ultrafine particle number and 0\.83 and 0\.82 for the condensation sink at Hyytiälä and Kumpula, respectively\.

Most of the improvement comes from conditioning itself\. For ultrafine particle number, the CRPS falls from 250 cm\-3for the historical ensemble to 211 cm\-3for FM at Hyytiälä, and from 1028 cm\-3to 739 cm\-3at Kumpula\. GDE then lowers it further, to 181 cm\-3\(−14\-14%\) and 599 cm\-3\(−19\-19%\)\. The PM2\.5and meteorological record therefore carries size information beyond the season and hour of day\.

The intervals are informative but too narrow\. For ultrafine particle number, FM\+GDE’s mean 80 % interval is 391 cm\-3wide at Hyytiälä and 1583 cm\-3at Kumpula, against 1158 cm\-3and 4070 cm\-3for the historical ensemble, but it contains only 49 % and 52 % of the observations; for the condensation sink, coverage is 50 % and 49 %\. Notably, the GDE constraint widens plain FM’s intervals and raises its coverage for ultrafine particle number from 28 % to 49 % at Hyytiälä and from 29 % to 52 % at Kumpula\. The constraint thus improves calibration as well as accuracy, although no learned model reaches the nominal 80 %\.

Figure 4:Locally trained models add skill beyond the historical ensemble for both quantities, but their intervals remain too narrow\.Top: CRPSS for ultrafine particle number and the resolved condensation sink; zero marks the historical ensemble and positive values indicate improvement\. Bottom: central\-80 % coverage, with the nominal 0\.80 marked\. Circles denote Hyytiälä and triangles Kumpula\.Table 5:Derived\-quantity scores for local and transferred models on identical target\-site windows\.CRPS and central\-80 % coverage \(C80C\_\{80\}, nominal 0\.80\) are given for ultrafine particle \(UFP\) number and the condensation sink\. Interval widths and CRPSS are in Table[F](https://arxiv.org/html/2610.33557#A6)\.\(a\) Hyytiälä target

\(b\) Kumpula target

### 3\.4RQ4: Transferability depends on both the quantity and the direction of transfer

A station without its own PSD archive would have to rely on a model trained elsewhere\. Transfer between our two stations is the most demanding test in this benchmark, because each model has seen only one aerosol regime during training\. We use it as a stress test that reveals where single\-station training breaks down, not as an estimate of what a model trained on many stations could achieve\. To measure the cost of transfer, we apply each model, without retraining, to the other site and divide its CRPS by that of the corresponding locally trained model on the same windows \(Fig\.[5](https://arxiv.org/html/2610.33557#S3.F5)\)\. A ratio of one means that transfer costs nothing\.

Urban to boreal\.From Kumpula to Hyytiälä, the condensation sink degrades less than ultrafine particle number for all four learned models, with CRPS ratios of 1\.2–2\.5 against 2\.7–6\.1\. Diffusion transfers best: its CRPS for the sink rises 1\.48\-fold, from0\.90×10−40\.90\\times 10^\{\-4\}s\-1to1\.33×10−41\.33\\times 10^\{\-4\}s\-1, while its CRPS for ultrafine particle number rises 2\.74\-fold, from 194 cm\-3to 532 cm\-3\. Plain FM transfers worse on both counts, with a 6\.07\-fold degradation for ultrafine particle number and a 2\.54\-fold degradation for the sink\. Good point accuracy after transfer does not imply good calibration: diffusion’s transferred coverage for the sink falls to 0\.28, from 0\.53 locally\.

Boreal to urban\.From Hyytiälä to Kumpula, the two quantities degrade by comparable amounts, with no consistent ordering between them, and this direction shows the more severe failure\. For ultrafine particle number, every Hyytiälä\-trained model is less skilful than the Hyytiälä\-derived historical ensemble itself \(CRPSS−0\.07\-0\.07to−0\.36\-0\.36; Table[F](https://arxiv.org/html/2610.33557#A6)c\), and coverage collapses to 0\.07–0\.10\. Plain FM’s CRPS for ultrafine particle number rises from 739 cm\-3with local training to 1896 cm\-3after transfer\. In the example of Fig\.[15](https://arxiv.org/html/2610.33557#A5.F15), no transferred model reproduces the observed daytime enhancement of small particles at Kumpula\. Condensation sink fares better: transferred FM\+GDE has a CRPSS of\+0\.39\+0\.39, although its CRPS \(3\.26×10−43\.26\\times 10^\{\-4\}s\-1\) remains more than twice that of local FM\+GDE \(1\.49×10−41\.49\\times 10^\{\-4\}s\-1\)\.

Effect of the constraint\.GDE lowers FM’s transferred CRPS in all four direction–quantity cases, by 4–30 %, but it does not close the gap to local training in any of them\. Skill scores against a transferred reference must also be read with care\. For ultrafine particle number, the Kumpula\-derived historical ensemble scores 2010 cm\-3at Hyytiälä, against 250 cm\-3for the local ensemble, so the large positive CRPSS in the urban\-to\-boreal direction partly reflects a weak reference rather than successful transfer\. The paired ratio against local training avoids this ambiguity\. Full\-period spectra, fitted\-component ratios and complete quantity diagnostics are given in Appendices[E](https://arxiv.org/html/2610.33557#A5)–[G](https://arxiv.org/html/2610.33557#A7)\.

Figure 5:Transfer costs more skill for ultrafine particle number than for the condensation sink from Kumpula to Hyytiälä, but not in the reverse direction\.Points show transferred CRPS divided by the same model’s locally trained CRPS on identical target\-site windows\. The dashed line at one denotes no loss; larger ratios indicate degradation\. The horizontal axis is logarithmic\.
### 3\.5Comparison with the SOSAA\-FP process model

To place the data\-driven spectra in context, we compare them with SOSAA\-FP, a one\-dimensional model that explicitly represents chemistry and aerosol microphysics\([Boy et al\., 2011](https://arxiv.org/html/2610.33557#bib.bib3);[Clusius et al\., 2026](https://arxiv.org/html/2610.33557#bib.bib5)\)\. We use the SOSAA\-FP runs of[Clusius et al\. \(2026\)](https://arxiv.org/html/2610.33557#bib.bib5)for 1 March 2018, 02:00 UTC to 1 November 2018, 01:00 UTC \(5880 hourly rows\)\. Their 60\-bin, 1\.07–2000 nm grid contains our 50 diameters, so the spectra can be matched without interpolation\. The matched windows are drawn from the stratified interpolation holdout, which shares observations with neighbouring training windows\.

On 645 matched windows, FM has a correlation of 0\.938 and an RMSE of 1002 cm\-3, against a correlation of 0\.715 and an RMSE of 2105 cm\-3for SOSAA\-FP \(Table[6](https://arxiv.org/html/2610.33557#S3.T6)\)\. The MLP has the lowest RMSE \(867 cm\-3\)\. In two matched cases, SOSAA\-FP underestimates concentrations across most sizes by about an order of magnitude, whereas FM follows the observed amplitude above about 10 nm \(Figs\.[6](https://arxiv.org/html/2610.33557#S3.F6)and[7](https://arxiv.org/html/2610.33557#S3.F7)\)\. This comparison comes with several caveats\. The data\-driven models were fitted to local DMPS records, whereas SOSAA\-FP is a simulation driven by emissions, chemistry and meteorology\. It can therefore diagnose processes and sources that the learned models cannot\. The result shows only that routine observations reproduce the local spectrum more closely in this matched setting\. The running\-mean and seasonal diagnostic of[Clusius et al\. \(2026\)](https://arxiv.org/html/2610.33557#bib.bib5)is given in Appendix[H](https://arxiv.org/html/2610.33557#A8)\.

Table 6:On matched 2018 windows at Hyytiälä, the learned models agree more closely with DMPS observations than SOSAA\-FP does\.Metrics use 645 windows and all 50 common diameter bins \(2\.82–794 nm\); RMSE is in cm\-3\. The “vs SOSAA\-FP” rows measure disagreement between models and are not skill scores\.\\belowtable

SOSAA\-FP is a process\-model simulation; the data\-driven models were fitted to the local measurement record\. Agreement with SOSAA\-FP is therefore not a measure of accuracy\.

Figure 6:In two matched cases, Flow Matching follows the observed spectral amplitude, whereas SOSAA\-FP is about an order of magnitude lower\.Hyytiälä, \(a\) 10 June 2018, 23:00 UTC, and \(b\) 12 August 2018, 10:00 UTC\. Black: observed DMPS spectrum; salmon dashed: Flow Matching; red dashed: SOSAA\-FP\. Both use the common 2\.82–794 nm interval\.![Refer to caption](https://arxiv.org/html/2610.33557v1/fig_sosaa_heatmaps.png)Figure 7:The same contrast holds through the day\.Time\-resolved PSDs for the cases of Fig\.[6](https://arxiv.org/html/2610.33557#S3.F6): \(a–c\) 10 June 2018 and \(d–f\) 12 August 2018\. Within each row, panels show the observed DMPS distribution, SOSAA\-FP and Flow Matching\. Color giveslog10⁡\(d​N/d​log⁡Dp\)\\log\_\{10\}\(\\mathrm\{d\}N/\\mathrm\{d\}\\log D\_\{p\}\)on a common scale within each row\.
## 4Discussion

### 4\.1Interpretation

In summary, conditioning a flow\-matching model on PM2\.5and meteorology yields spectra whose distribution, not only whose mean, tracks the observations, and a reduced GDE residual improves them further\. Since this study is a benchmark on a two\-station testbed, the value of these results lies less in the absolute scores than in what they reveal about model families, constraints and failure modes\. Three results call for interpretation: the different verdicts of RMSE and CRPS, the contrast between constraint types, and the dependence of transfer on quantity and direction\.

Why ensembles beat regression on CRPS but not always on RMSE\.Under the mass–number ambiguity, several distinct spectra are consistent with the same input\. A squared\-error regression returns their average, which RMSE rewards but which lacks the structure and spread of any single realization \(Fig\.[2](https://arxiv.org/html/2610.33557#S2.F2)\)\. A generative ensemble spreads probability across the plausible spectra, which CRPS rewards\. The two metrics therefore answer different questions\. Gradient boosting’s lower RMSE at Hyytiälä and the generative models’ lower CRPS at both sites are consistent: the former is the better point estimate at that site, the latter the better description of what remains uncertain\.

Why the dynamical constraint helps and the integral constraints do not\.An integral constraint fixes one summary number and leaves many spectra that satisfy it\. In the nowcast, the mass reference also duplicates information already given as input\. When the resolved spectrum and the instrument do not sample the same size range, forcing agreement can distort the spectrum\. The GDE residual instead constrains how neighbouring bins and hours relate to each other, and it also widens the too\-narrow FM ensemble\. These observations are consistent with a regularizing effect of the dynamics, but the evidence is indirect: the constrained runs also have more parameters and a longer training budget\.

Why transferability depends on quantity and direction\.Number and condensation sink weight the same spectrum differently\. The sink is weighted towards larger particles, which tend to be regionally mixed and therefore shared between the two sites, whereas ultrafine particle number is dominated by the smallest particles, which reflect more local sources such as traffic at Kumpula and particle formation at Hyytiälä\. We hypothesize that this explains why the sink transfers better from Kumpula to Hyytiälä\. The boreal\-to\-urban failure is consistent with a missing source: the Hyytiälä record contains no analogue of the urban daytime small\-particle enhancement, and a soft dynamical constraint cannot supply one\. The present two\-site comparison cannot separate the roles of emissions, formation, growth, transport and mixing, so these explanations remain hypotheses\. A practical implication follows regardless: a PSD model should be judged by the quantity and uncertainty that its application needs, because full\-spectrum correlation can remain high while the relevant integral degrades\.

### 4\.2Limitations

Several limitations follow from the deliberate choice of a small, controlled testbed, and should be read in that light\. The evidence rests on two stations, one late\-year holdout of about 60 days and one training run per configuration, so the bootstrap intervals do not include training\-run variability; repeated training and holdouts from other seasons would address this\. Two stations are enough to compare model families under identical conditions and to expose failure under transfer, but not to show how far neural networks can generalize when trained on many aerosol regimes\. That question requires large, multi\-station training data and is the natural next step\. Every learned ensemble is too narrow, with 80 % intervals containing about half of the observations locally and fewer after transfer, and eight members limit how finely quantiles can be resolved\. Post\-hoc calibration, larger ensembles and guidance selected on validation data rather than the exploratory test\-set sweep \(Appendix[D\.4](https://arxiv.org/html/2610.33557#A4.SS4)\) are natural remedies\. The fitting period of the shared PSD normalization has not been verified, so a scaling leak from test data is not yet ruled out, and quality flags were not applied \(Appendix[A](https://arxiv.org/html/2610.33557#A1)\)\. The high\-number episode analysis contains only 0, 3, 29 and 0 qualifying days in the four settings, too few for reliable probability verification \(Appendix[F\.1](https://arxiv.org/html/2610.33557#A6.SS1)\)\.

The physical interpretation is also incomplete\. The GDE gain cannot yet be attributed to the residual itself, because the constrained runs have an auxiliary head and a longer training budget, and the sulfuric\-acid proxy in the growth term depends on SO2, which is not among the stated inputs \(Appendix[C](https://arxiv.org/html/2610.33557#A3)\); capacity\- and budget\-matched controls are needed\. The integral constraints use references with different measurement support, and in the nowcast the pseudo\-mass reference is the input itself, so a size\-resolved conversion between DMPS, PM2\.5and CPC support would be required to make them physically meaningful\. Finally, the derived quantities are computed from the same DMPS spectra used for evaluation and so do not validate process physics independently, and the SOSAA\-FP comparison uses windows that overlap training data and a single model configuration, so it cannot rank data\-driven and process\-based modelling in general\.

\\conclusions

We asked how much of the aerosol particle size distribution can be recovered from the concurrent PM2\.5and meteorological record\. Treating this nowcast as conditional generation, we benchmarked six generative models and three physical constraints on a controlled two\-station testbed, with the aim of identifying which modelling choices matter rather than proposing a single best model\. Four findings stand out\. \(1\) Generative ensembles give better distributional skill than deterministic and climatological estimates: all six generative models have lower CRPS than gradient boosting at both sites, although gradient boosting remains competitive on point error\. \(2\) The form of a physical constraint matters: a reduced GDE lowers CRPS by 11\.5 % and 5\.5 % and widens the ensemble, whereas integral constraints on mass or number alone degrade the spectra\. \(3\) The routine record adds information beyond historical spectra: FM\+GDE has the lowest CRPS of the learned models for ultrafine particle number and the condensation sink at both sites, with skill scores of 0\.27–0\.83 against the historical ensemble, but its intervals contain only about half of the observations\. \(4\) Transferability depends on the quantity and the direction: the condensation sink transfers better than ultrafine particle number from the urban to the boreal site, whereas boreal\-to\-urban transfer fails for ultrafine particle number\. Local PSD nowcasting from routine observations is therefore feasible today, and the testbed has shown which ingredients are worth carrying forward\.

These results raise the question of what would be required to nowcast PSDs at a station that measures only PM2\.5\. First, and most importantly, the method should be scaled up: trained on large data sets from many stations that cover contrasting aerosol regimes, neural networks can be tested where their main strength lies, in learning from diverse data and generalizing to new domains and out\-of\-distribution conditions\. The boreal\-to\-urban failure observed here, after training on a single station, is exactly the situation that multi\-station training is designed to overcome\. Second, calibrated uncertainty is a prerequisite for any use of the intervals\. Third, training\-derived climatology could be represented explicitly, for example by starting generation from a historical PSD template instead of pure noise\. Fourth, source\-sensitive predictors such as traffic or emission indicators could reduce the mass–number ambiguity directly\. Finally, capacity\-matched ablations would clarify how much of the constraint’s benefit is physical\. These improvements are complementary to process models: data\-driven nowcasts can extend the observational record, while process models remain necessary to explain it\. The evaluation framework presented here can be used to measure progress along each of these axes\.

## Appendix AData preparation and evaluation splits

### A\.1Size grid and data assembly

The DMPS archives at both stations report 61 diameter bins between 1\.00 and 1000 nm\. We use a fixed 50\-bin subset between 2\.82 and 794 nm, excluding the nine smallest archived bins \(1\.00–2\.51 nm\) and the two largest \(891 and 1000 nm\), and we keep the archived bin centres without interpolation\. Hourly values are means of the native 10 min measurements, each weighted equally\. ERA5 variables are taken from the grid cell nearest to each station\. After assembly, each station has 61 368 hourly records with complete PM2\.5, PSD and ERA5 values on a common time axis\.

### A\.2Transformations

Particle concentrations span several orders of magnitude across bins and over time, so the model operates on log\-transformed values\. For each bin, the concentrationnnis transformed toℓ=log⁡\(1\+n\)\\ell=\\log\(1\+n\)and then scaled linearly to the interval\[−1,1\]\[\-1,1\],

ℓ~=2​ℓ−ℓminℓmax−ℓmin−1,\\tilde\{\\ell\}=2\\,\\frac\{\\ell\-\\ell\_\{\\min\}\}\{\\ell\_\{\\max\}\-\\ell\_\{\\min\}\}\-1,\(13\)with fixed bin\-wise constantsℓmin\\ell\_\{\\min\}andℓmax\\ell\_\{\\max\}shared by both stations\. The same scaling therefore applies at the training site and at the transfer target\. Values outside the scaled range are kept rather than clipped\. To return to physical units, the scaling is inverted andn=exp⁡\(ℓ\)−1n=\\exp\(\\ell\)\-1is applied after negative values ofℓ\\ellare set to zero, which guarantees non\-negative concentrations\. PM2\.5is transformed withlog⁡\(1\+x\)\\log\(1\+x\), and positive, strongly skewed meteorological variables are also log\-transformed\.

### A\.3Examples and splits

An example is a 24 h window of inputs together with the PSD of the same 24 hours \(Sect\.[1](https://arxiv.org/html/2610.33557#S1)\)\. Windows start at every hour, so consecutive windows overlap by 23 hours, and a window is used only if all 24 hours are valid\. Missing hours are kept as gaps on the time axis rather than removed, so that every window spans exactly 24 consecutive clock hours\.

The data are split in time\. All training windows end before 2 November 2024, 00:00 UTC, and the temporal holdout consists of the 1417 windows that start on or after this cutoff \(60 days, 2 November–31 December 2024\)\. Before the cutoff, 10 % of the eligible windows \(6048 per station\), stratified by month and hour of day, form an*interpolation*split, and the remaining 53 857 windows per station are used for training\. Because the interpolation windows overlap with neighbouring training windows by up to 23 hours, the interpolation split measures how well a model reproduces the training period, not its skill on unseen periods\. We use it only for the diagnostics in Appendix[D\.2](https://arxiv.org/html/2610.33557#A4.SS2)and for the matched comparison with SOSAA\-FP \(Sect\.[3\.5](https://arxiv.org/html/2610.33557#S3.SS5)\); all skill statements in the main text refer to the temporal holdout\. No separate validation set or early stopping is used\.

## Appendix BGenerative models and baselines

### B\.1Flow matching

Flow matching connects a sample𝐱0\\mathbf\{x\}\_\{0\}of standard Gaussian noise and a data sample𝐱1\\mathbf\{x\}\_\{1\}by the straight path

𝐱t=\(1−t\)​𝐱0\+t​𝐱1,𝐯∗=d​𝐱td​t=𝐱1−𝐱0,\\mathbf\{x\}\_\{t\}=\(1\-t\)\\,\\mathbf\{x\}\_\{0\}\+t\\,\\mathbf\{x\}\_\{1\},\\hskip 18\.49988pt\\mathbf\{v\}^\{\\ast\}=\\frac\{\\mathrm\{d\}\\mathbf\{x\}\_\{t\}\}\{\\mathrm\{d\}t\}=\\mathbf\{x\}\_\{1\}\-\\mathbf\{x\}\_\{0\},\(14\)and trains the network𝐯θ​\(𝐱t,t,𝐜\)\\mathbf\{v\}\_\{\\theta\}\(\\mathbf\{x\}\_\{t\},t,\\mathbf\{c\}\)to predict the velocity𝐯∗\\mathbf\{v\}^\{\\ast\}from the intermediate state \(Eq\.[1](https://arxiv.org/html/2610.33557#S2.E1)\)\. Averaged over all pairs that pass through a given state, the regression target becomes the velocity of a transport that carries the whole noise distribution onto the conditional distribution of PSDs\([Lipman et al\., 2023](https://arxiv.org/html/2610.33557#bib.bib17);[Liu et al\., 2023](https://arxiv.org/html/2610.33557#bib.bib18)\)\. A new sample is generated by integrating this velocity field from a fresh noise draw\. We use the explicit Euler scheme withKKsteps of sizeΔ​t=1/K\\Delta t=1/K,

𝐱\(0\)∼𝒩\(0,I\),𝐱\(k\+1\)=𝐱\(k\)\+Δt𝐯θ\(𝐱\(k\),tk,𝐜\),tk=kΔt,k=0,…,K−1,\\mathbf\{x\}^\{\(0\)\}\\sim\\mathcal\{N\}\(0,I\),\\hskip 18\.49988pt\\mathbf\{x\}^\{\(k\+1\)\}=\\mathbf\{x\}^\{\(k\)\}\+\\Delta t\\,\\mathbf\{v\}\_\{\\theta\}\\\!\\left\(\\mathbf\{x\}^\{\(k\)\},t\_\{k\},\\mathbf\{c\}\\right\),\\hskip 9\.24994ptt\_\{k\}=k\\Delta t,\\hskip 9\.24994ptk=0,\\ldots,K\-1,\(15\)withK=50K=50\(Appendix[I](https://arxiv.org/html/2610.33557#A9)\)\. The integration is deterministic for a given starting point, so all variability within an ensemble comes from the different noise draws𝐱\(0\)\\mathbf\{x\}^\{\(0\)\}\. Figure[8](https://arxiv.org/html/2610.33557#A2.F8)shows how an initially random field is resolved into a coherent spectrum along this path\.

### B\.2Network architecture

The velocity network is a Transformer\([Vaswani et al\., 2017](https://arxiv.org/html/2610.33557#bib.bib26)\)that represents a 24 h window as a sequence of 24 hourly elements\. For each hour, the current 50\-bin state𝐱t\\mathbf\{x\}\_\{t\}and the inputs of that hour are projected into a common feature space and added, and a learned position embedding identifies the hour within the window\. Self\-attention layers then exchange information between all hours\. The generative timettenters through adaptive layer normalization, which rescales the features of every layer according tott; the network can therefore behave differently near the noise end and near the data end of the path\. A final projection maps each hour back to 50 bins\. The network has about1\.24×1061\.24\\times 10^\{6\}trainable parameters\. The diffusion model, the classifier\-free\-guidance variant and the optimal\-transport variant use the same network, so their differences reflect the generative objective rather than the architecture\. Training settings are provided with the code \(see Code availability\); for interpreting the comparisons, the relevant difference is that plain FM, diffusion, FM\+CFG and the MLP were trained for 200 epochs, whereas FM\+OT, the Drifting Model, Energy Matching and all constrained FM variants were trained for 400 epochs\.

### B\.3Comparison generative models

Diffusion\.A denoising diffusion model\([Ho et al\., 2020](https://arxiv.org/html/2610.33557#bib.bib12)\)corrupts a data sample𝐲\\mathbf\{y\}gradually with Gaussian noise,𝐱s=α¯s​𝐲\+1−α¯s​ϵ\\mathbf\{x\}\_\{s\}=\\sqrt\{\\bar\{\\alpha\}\_\{s\}\}\\,\\mathbf\{y\}\+\\sqrt\{1\-\\bar\{\\alpha\}\_\{s\}\}\\,\\boldsymbol\{\\epsilon\}, whereα¯s\\bar\{\\alpha\}\_\{s\}decreases from one to nearly zero over the noising stepsssfollowing a linear schedule\. The network is trained to predict the added noiseϵ\\boldsymbol\{\\epsilon\}, and samples are generated by reversing the corruption step by step, adding fresh noise at each step\. Unlike FM, the sampling path is stochastic\.

Optimal\-transport coupling \(FM\+OT\)\.Standard FM pairs noise and data samples at random, so the straight paths of different pairs can cross\. FM\+OT instead pairs the samples within each training batch by solving an assignment problem that minimizes the total squared distance between paired samples\([Tong et al\., 2024](https://arxiv.org/html/2610.33557#bib.bib24)\)\. The resulting paths are shorter and cross less often, which can make the learned velocity field smoother\.

Classifier\-free guidance \(FM\+CFG\)\.During training, the conditioning𝐜\\mathbf\{c\}is replaced by a learned null input∅\\emptysetfor 15 % of the examples, so that one network learns both the conditional and the unconditional velocity\([Ho and Salimans, 2021](https://arxiv.org/html/2610.33557#bib.bib11)\)\. At sampling time, the two are combined as

𝐯~=𝐯θ​\(𝐱t,t,∅\)\+w⁡\[𝐯θ​\(𝐱t,t,𝐜\)−𝐯θ​\(𝐱t,t,∅\)\]\.\\tilde\{\\mathbf\{v\}\}=\\mathbf\{v\}\_\{\\theta\}\(\\mathbf\{x\}\_\{t\},t,\\emptyset\)\+w\\left\[\\mathbf\{v\}\_\{\\theta\}\(\\mathbf\{x\}\_\{t\},t,\\mathbf\{c\}\)\-\\mathbf\{v\}\_\{\\theta\}\(\\mathbf\{x\}\_\{t\},t,\\emptyset\)\\right\]\.\(16\)A guidance scalew=1w=1recovers ordinary conditional sampling,w\>1w\>1strengthens the influence of the inputs, andw<1w<1moves the samples towards the unconditional distribution, which widens the ensemble\. The tables usew=1w=1; the effect ofwwis shown in Appendix[D\.4](https://arxiv.org/html/2610.33557#A4.SS4)\.

Drifting Model\.The Drifting Model\([Deng et al\., 2026](https://arxiv.org/html/2610.33557#bib.bib7)\)learns a generator that maps noise to a sample in a single network evaluation\. During training, the distribution of generated samples is moved step by step towards the data distribution, so that no integration is needed at inference\.

Energy Matching\.Energy Matching\([Balcerak et al\., 2025](https://arxiv.org/html/2610.33557#bib.bib1)\)learns a scalar potential whose negative gradient defines the transport from noise to data\. Because the potential does not depend on the generative time, the network does not receivettas an input\. We train it with its flow\-based objective only; the additional contrastive refinement stage of the method was not used, because it was unstable on these data\.

### B\.4Deterministic baselines and non\-learning references

Gradient boosting fits one regression model per diameter bin\. Each hour of a training window is treated as a separate sample, with the hour of day encoded by its sine and cosine, so the model can adapt its estimate to the time of day but does not represent dependence between bins or between hours\. The MLP is a fully connected network that maps the complete 24 h input window to all24×5024\\times 50outputs at once\. Both are trained on squared error and return a single PSD per input\.

The month/hour climatology is the mean training\-period PSD for each of the 288 combinations of calendar month and hour of day; each target hour is assigned the value of its own month and hour\. The seasonal template is defined in the same way using seasons instead of months\. The historical ensemble consists of eight complete 24 h PSD sequences drawn from the training period of the source station, each starting in the same season and at the same hour of day as the target window\. Because whole sequences are used, every member retains the observed dependence across bins and hours\. The historical ensemble ignores the current inputs and therefore shows what can be estimated from the training archive alone\.

![Refer to caption](https://arxiv.org/html/2610.33557v1/fig_fm_trajectory.png)Figure 8:Flow Matching resolves Gaussian noise into a coherent 24 h spectrum during integration\.Ensemble\-mean state along the Euler integration ofd​𝐱/d​t=𝐯θ\\mathrm\{d\}\\mathbf\{x\}/\\mathrm\{d\}t=\\mathbf\{v\}\_\{\\theta\}for one temporal\-holdout sample \(10 December 2024\), from an almost purely Gaussian field at step 10 to the generated distribution at step 50, alongside the observation\. Color islog10⁡d​N/d​log⁡Dp\\log\_\{10\}\\,\\mathrm\{d\}N/\\mathrm\{d\}\\log D\_\{p\}; the vertical axis is mobility diameter\.

## Appendix CPhysical constraints

### C\.1Resolved mass and number

The resolved mass and number \(Eqs\.[3](https://arxiv.org/html/2610.33557#S2.E3)and[5](https://arxiv.org/html/2610.33557#S2.E5)\) are sums over the 50 bins of the predicted PSD, using logarithmic bin widthsΔi=Δ​log10​Dp,i\\Delta\_\{i\}=\\Delta\\log\_\{10\}D\_\{p,i\}and a constant particle densityρ=1\.5\\rho=1\.5g cm\-3\. PM2\.5is defined by the aerodynamic diameter, whereas the DMPS sizes particles by their electrical mobility\. For spherical particles, the two diameters are approximately related byDm≈Da\(ρ/ρ0\)−1/2D\_\{\\rm m\}\\approx D\_\{\\rm a\}\(\\rho/\\rho\_\{0\}\)^\{\-1/2\}withρ0=1\\rho\_\{0\}=1g cm\-3, neglecting shape and slip corrections\. The PM2\.5cut therefore corresponds to a mobility diameter of about2\.5μm×1\.5−1/2≈2\.0μ2\.5\\,\\mu\\mathrm\{m\}\\times 1\.5^\{\-1/2\}\\approx 2\.0\\,\\mum, well above the 794 nm upper limit of the grid\. All predicted bins thus lie inside the PM2\.5size range, but the mass between 794 nm and about 2μ\\mum is missing from the resolved mass, which is why the mass term is used only as a soft reference\.

Both integral constraints compare logarithms of the resolved and the measured quantity through the Huber function,

Huberδ⁡\(u\)=\{u2/2,\|u\|≤δ,δ⁡\(\|u\|−δ/2\),\|u\|\>δ,δ=0\.5,\\operatorname\{Huber\}\_\{\\delta\}\(u\)=\\begin\{cases\}u^\{2\}/2,&\|u\|\\leq\\delta,\\\\ \\delta\(\|u\|\-\\delta/2\),&\|u\|\>\\delta,\\end\{cases\}\\hskip 18\.49988pt\\delta=0\.5,\(17\)whereuuis the difference of the two logarithms\. In logarithmic space, a constant multiplicative mismatch, such as the missing mass above 794 nm, becomes a constant offset, and the linear branch of the Huber function prevents such systematic offsets from dominating the gradient\. The reference scalesM0M\_\{0\}andN0N\_\{0\}in Eqs\. \([4](https://arxiv.org/html/2610.33557#S2.E4)\) and \([5](https://arxiv.org/html/2610.33557#S2.E5)\) make the logarithms well defined at zero concentration, and negative reference values, which can arise from instrument noise, are set to zero\. The number constraint is evaluated only in hours with a valid CPC measurement\.

Each physical term is weighted so that its initial magnitude is comparable to that of the FM loss; for the mass term, this givesλm=0\.05\\lambda\_\{\\rm m\}=0\.05, and the number and GDE weights were set by the same rule\. A mass weight of 1\.0 led to a degenerate solution that satisfied the mass constraint at the expense of a realistic spectrum, which illustrates why the mass term must remain soft\.

The CPC and the DMPS measure closely related but not identical quantities\. Over the full record, their log\-transformed total number concentrations correlate with Pearsonr=0\.988r=0\.988at Hyytiälä \(59 689 hours with both measurements\) andr=0\.963r=0\.963at Kumpula \(57 424 hours\), with median DMPS/CPC ratios of 0\.87 and 1\.03, respectively\. Possible contributors to the differences include particles below the DMPS lower limit, which the CPC counts, and differences in detection efficiency and sampling losses\([Wiedensohler et al\., 2012](https://arxiv.org/html/2610.33557#bib.bib27)\), which the number constraint does not model\.

### C\.2Condensation and coagulation sinks

The condensation sink \(Eq\.[12](https://arxiv.org/html/2610.33557#S2.E12)\) is the rate at which a condensable vapor is lost to the particle population\([Kulmala et al\., 2001](https://arxiv.org/html/2610.33557#bib.bib15)\)\. The Fuchs–Sutugin correction accounts for the transition between the continuum regime, where vapor diffuses to large particles, and the free\-molecular regime of small particles,

βm​\(Dp\)=1\+Kn1\+\(43​α\+0\.377\)​Kn\+43​α​Kn2,Kn=2​λvDp,\\beta\_\{m\}\(D\_\{p\}\)=\\frac\{1\+\\mathrm\{Kn\}\}\{1\+\\left\(\\frac\{4\}\{3\\alpha\}\+0\.377\\right\)\\mathrm\{Kn\}\+\\frac\{4\}\{3\\alpha\}\\mathrm\{Kn\}^\{2\}\},\\hskip 18\.49988pt\\mathrm\{Kn\}=\\frac\{2\\lambda\_\{v\}\}\{D\_\{p\}\},\(18\)whereλv\\lambda\_\{v\}is the mean free path of the vapor molecules andα\\alphathe accommodation coefficient, set to one\. For small particles,βm​Dp\\beta\_\{m\}D\_\{p\}grows approximately withDp2D\_\{p\}^\{2\}, and for large particles approximately withDpD\_\{p\}, so the sink is weighted towards larger particles\. The coagulation sink of a particle of diameterdd,

CoagS⁡\(d\)=∑iK⁡\(d,Dp,i\)​Ni,\\mathrm\{CoagS\}\(d\)=\\sum\_\{i\}K\(d,D\_\{p,i\}\)\\,N\_\{i\},\(19\)is its loss rate by coagulation with the existing particles, whereKKis the Brownian coagulation kernel\([Kulmala et al\., 2001](https://arxiv.org/html/2610.33557#bib.bib15);[Seinfeld and Pandis, 2016](https://arxiv.org/html/2610.33557#bib.bib23)\)andNi=ni​Δ​log10​Dp,iN\_\{i\}=n\_\{i\}\\Delta\\log\_\{10\}D\_\{p,i\}is the number concentration in binii\. Both sinks are evaluated at fixed reference conditions \(293\.15 K, 101325 Pa\) for predicted and observed spectra alike, so that differences between them reflect the size distribution only\.

### C\.3Reduced general dynamic equation

For a size distributionn=d​N/d​log10​Dpn=\\mathrm\{d\}N/\\mathrm\{d\}\\log\_\{10\}D\_\{p\}at a fixed station, the continuity equation underlying the GDE constraint is

∂n∂τ\+∂∂ln⁡Dp​\[GRDp​n\]=S−CoagS​n−λdil​n,\\frac\{\\partial n\}\{\\partial\\tau\}\+\\frac\{\\partial\}\{\\partial\\ln D\_\{p\}\}\\left\[\\frac\{\\mathrm\{GR\}\}\{D\_\{p\}\}n\\right\]=S\-\\mathrm\{CoagS\}\\,n\-\\lambda\_\{\\mathrm\{dil\}\}n,\(20\)whereτ\\tauis physical \(clock\) time\. The second term on the left describes condensational growth as a flux through size space with velocityGR/Dp\\mathrm\{GR\}/D\_\{p\}inln⁡Dp\\ln D\_\{p\}\. On the right,SScollects sources and transport, the coagulation term removes particles that collide with the existing population, andλdil=\(d​H/d​τ\)\+/H\\lambda\_\{\\mathrm\{dil\}\}=\(\\mathrm\{d\}H/\\mathrm\{d\}\\tau\)\_\{\+\}/Hdilutes the population when the boundary layer of heightHHgrows\. The equation is reduced in three respects: coagulation enters only as a loss, not as a gain at larger sizes; deposition is not represented explicitly; and new particle formation appears only as the growth flux into the smallest resolved bin\. The source termSSis unconstrained in sign and can absorb these omitted processes, as well as the advection of different air masses past the station\.

The residual of Eq\. \([20](https://arxiv.org/html/2610.33557#A3.E20)\) is evaluated on the predicted spectrum \(Eq\.[6](https://arxiv.org/html/2610.33557#S2.E6)\)\. Time derivatives are computed by finite differences over the 24 hourly values, central in the interior and one\-sided at the first and last hour, and size derivatives by upwind differences inln⁡Dp\\ln D\_\{p\}, which follow the direction of growth\. At the lower edge of the grid, the incoming flux is set by the growth velocity of the first bin; it represents growth into the observed range rather than a nucleation rate\. All rates are expressed per hour\. The growth rate and the source term are predicted by an auxiliary network head from the conditioning only, so that they do not inherit noise from the intermediate generative state\. The growth rate is constrained to be non\-negative and is informed by a sulfuric\-acid proxy proportional to SO2concentration times radiation divided by the condensation sink\([Petäjä et al\., 2009](https://arxiv.org/html/2610.33557#bib.bib21)\); the dependence of this proxy on SO2, which is not among the model inputs, is discussed in Sect\.[4](https://arxiv.org/html/2610.33557#S4)\. The head also predicts the log residual variancesh,is\_\{h,i\}of Eq\. \([7](https://arxiv.org/html/2610.33557#S2.E7)\), which is bounded to\[−5,5\]\[\-5,5\], i\.e\. to variances between6\.7×10−36\.7\\times 10^\{\-3\}and1\.5×1021\.5\\times 10^\{2\}, so that the model can neither ignore a residual entirely nor overweight it without limit\. BecauseSSandsh,is\_\{h,i\}are flexible, a low residual shows consistency with Eq\. \([20](https://arxiv.org/html/2610.33557#A3.E20)\) but does not uniquely identify the growth or source rates\.

### C\.4Coupling the constraints to flow matching

During training, the network sees partly noised states𝐱t\\mathbf\{x\}\_\{t\}, whereas the physical relationships apply to clean spectra\. Inverting the straight path of Eq\. \([14](https://arxiv.org/html/2610.33557#A2.E14)\) gives the clean\-sample estimate implied by the current velocity prediction,

𝐱^1=𝐱t\+\(1−t\)​𝐯θ,\\hat\{\\mathbf\{x\}\}\_\{1\}=\\mathbf\{x\}\_\{t\}\+\(1\-t\)\\,\\mathbf\{v\}\_\{\\theta\},\(21\)which equals𝐱1\\mathbf\{x\}\_\{1\}exactly when the predicted velocity is correct\. The constraints are evaluated on𝐱^1\\hat\{\\mathbf\{x\}\}\_\{1\}after conversion to physical units, and their gradients flow back into𝐯θ\\mathbf\{v\}\_\{\\theta\}\. Near the noise end \(t→0t\\to 0\),𝐱^1\\hat\{\\mathbf\{x\}\}\_\{1\}is an extrapolation from an almost random state, so the weightst2t^\{2\}andt3t^\{3\}of Eq\. \([2](https://arxiv.org/html/2610.33557#S2.E2)\) suppress the constraints there; the GDE term, which contains derivatives that amplify noise, receives the steeper weight\. The constraints act only on the training objective; generation at inference proceeds exactly as for unconstrained FM\. All physical ablations use FM as the base model, and every constrained variant uses this time weighting\.

## Appendix DEvaluation metrics and diagnostics

### D\.1Scores and their aggregation

For generative models, point scores \(RMSE and correlation\) use the ensemble mean\. The correlation is computed for each window on the log\-transformed spectrum, treating its50×2450\\times 24values as one vector, and then averaged over windows\. It measures agreement in spectral shape and temporal evolution, and it can be high simply because the typical shape of the PSD persists from window to window\.

The CRPS of Eq\. \([9](https://arxiv.org/html/2610.33557#S2.E9)\) uses the fair estimator for a finite ensemble ofM=8M=8members, which corrects the tendency of the ordinary estimator to penalize small ensembles; for deterministic predictions \(M=1M=1\), the pairwise term vanishes and the CRPS reduces to the absolute error\. The CRPS scores each bin and hour separately and therefore does not test whether the dependence between bins and hours is realistic\.

Coverage is computed from the empirical 10th and 90th percentiles of the eight members\. With eight members, these percentiles are coarse, so coverage values should be read as approximate\. Coverage can be aggregated over bins in two ways: Table[7](https://arxiv.org/html/2610.33557#A4.T7)reports the median of the per\-bin coverage \(0\.454 for FM\), whereas the guidance sweep in Appendix[D\.4](https://arxiv.org/html/2610.33557#A4.SS4)reports the mean \(0\.528 for FM\)\. The rank of an observation within its ensemble counts the members that lie strictly below it; for a calibrated ensemble, all ranks are equally likely\. The spread–skill ratio is the mean ensemble standard deviation divided by the RMSE of the ensemble mean; values near one indicate an ensemble whose spread matches its error, and values below one indicate an ensemble that is too narrow\.

Confidence intervals for the mean CRPS use a block bootstrap with 1000 replicates\. Because consecutive windows overlap, the data are resampled in blocks of 24 consecutive windows, and the 2\.5th and 97\.5th percentiles of the resampled means form the 95 % interval\. These intervals describe sampling uncertainty over windows, not variation between training runs\. Ensemble members are drawn anew in each evaluation, so two evaluations of the same trained model differ slightly, as for FM in Tables[3](https://arxiv.org/html/2610.33557#S3.T3)and[4](https://arxiv.org/html/2610.33557#S3.T4)\(CRPS of 142 and 141\.9 cm\-3\)\.

### D\.2Interpolation, variability and calibration

Three diagnostics describe how much of the observed variability a model reproduces\. The standard\-deviation ratio compares the standard deviation of the predicted and observed spectra across windows\. The texture ratio compares their bin\-to\-bin variability within a window, which measures how much fine spectral structure is retained\. The temporal correlation is the correlation between the predicted and observed time series of total number concentration\. For generative models, all three use the ensemble mean, which is expected to be smoother than any single observation, so ratios below one partly reflect averaging rather than model deficiency\.

On the overlapping interpolation split, the MLP reaches a correlation of 0\.95, a value that reflects the overlap with training windows rather than skill on new periods \(Figs\.[9](https://arxiv.org/html/2610.33557#A4.F9)and[10](https://arxiv.org/html/2610.33557#A4.F10)\)\. On the temporal holdout, FM retains standard\-deviation and texture ratios of 0\.83 and 0\.75, compared with 0\.58 and 0\.50 for gradient boosting \(Table[7](https://arxiv.org/html/2610.33557#A4.T7), Fig\.[11](https://arxiv.org/html/2610.33557#A4.F11)\)\. The rank histograms of all generative models depart from uniformity \(Fig\.[12](https://arxiv.org/html/2610.33557#A4.F12)\), indicating ensembles that are too narrow, conditionally biased, or both\.

![Refer to caption](https://arxiv.org/html/2610.33557v1/fig_grid_interp.png)Figure 9:On the overlapping interpolation split, all learned models reproduce the observed structure closely\.Illustrative Hyytiälä example\. The first panel is the observation; subsequent panels show one draw from each generative model or the deterministic prediction\. Color islog10⁡d​N/d​log⁡Dp\\log\_\{10\}\\,\\mathrm\{d\}N/\\mathrm\{d\}\\log D\_\{p\}; the vertical axis is mobility diameter and the horizontal axis is UTC hour\.Figure 10:Scores on the overlapping interpolation split are not estimates of temporal skill\.Hyytiälä\. \(a\) Mean correlation of log\-transformed spectra on the interpolation split and temporal holdout; lines connect the same model\. \(b\) Temporal\-holdout CRPS with 95 % block\-bootstrap intervals; the dashed line denotes climatology\.Figure 11:Flow Matching preserves more of the observed PSD variability than gradient boosting\.Temporal holdout, Hyytiälä\. \(a\) Ratio of predicted to observed standard deviation across windows and \(b\) ratio of predicted to observed within\-window spectral texture; unity is ideal\. \(c\) Taylor diagram combining the standard\-deviation ratio \(radius\) with temporal correlation \(angle\); the black star marks the observation\. Generative models use the ensemble mean\.Table 7:Distribution diagnostics for the Hyytiälä temporal holdout\. Variability ratios use the ensemble mean\. Coverage is the median across bins of central\-80 % interval coverage \(nominal 0\.80\); definitions are in Appendix[D](https://arxiv.org/html/2610.33557#A4)\.![Refer to caption](https://arxiv.org/html/2610.33557v1/figures/rank_histograms.png)Figure 12:No generative ensemble is calibrated\.Ensemble rank histograms on the temporal holdout\. Uniform ranks are expected for exchangeable observations and members \(dashed line\)\. A U\-shape can reflect under\-dispersion or a mixture of conditional biases; a dome indicates over\-dispersion\. Coverage is given in Table[7](https://arxiv.org/html/2610.33557#A4.T7)\.
### D\.3Generative objectives and small\-particle errors

The generative objectives trade different aspects of skill\. At Hyytiälä, optimal\-transport coupling slightly improves the correlation \(from 0\.820 to 0\.822\) and the RMSE \(from 389 cm\-3to 383 cm\-3\) of FM but worsens its CRPS \(from 142 cm\-3to 146 cm\-3\), and the one\-step Drifting Model reaches a CRPS of 152 cm\-3and an RMSE of 425 cm\-3at a fraction of the sampling cost\. No objective is best in accuracy, coverage and cost at once\.

Errors in the smallest bins are of particular interest because these bins carry the signature of new particle formation\. We therefore examine the ten smallest bins \(2\.82–7\.94 nm\) separately \(Fig\.[13](https://arxiv.org/html/2610.33557#A4.F13), Table[8](https://arxiv.org/html/2610.33557#A4.T8)\)\. Windows whose maximum observed concentration in these bins reaches at least the 75th percentile of the holdout form a peak group \(356 of 1417 windows, 25\.1 %\); the remaining windows form the background group\. Errors concentrate between 11 and 17 UTC, a period selected by inspecting the holdout errors, so this analysis is descriptive\. FM has the lowest RMSE in the peak group \(188\.1 cm\-3, against 189\.9 cm\-3for gradient boosting\), whereas gradient boosting has the lowest RMSE in the background group \(63\.6 cm\-3\)\. The peak group is defined by concentration alone and does not identify new particle formation events\.

Figure 13:Errors in the smallest bins concentrate around midday for every model\.Exploratory diagnostics for the ten smallest resolved bins \(2\.82–7\.94 nm\) at Hyytiälä\. \(a\) Signed bias by UTC hour; the 11–17 UTC window was selected by inspecting these errors\. \(b\) Background and peak\-window RMSE; peak windows reach at least the holdout 75th\-percentile amplitude \(n=356n=356of 1417\)\. Values are in Table[8](https://arxiv.org/html/2610.33557#A4.T8)\.Table 8:Diurnal and smallest\-bin peak\-window errors at Hyytiälä for the ten smallest resolved bins \(2\.82–7\.94 nm\), averaged over 11–17 UTC\. Peak windows are the top quartile by smallest\-bin amplitude \(n=356n=356of 1417\)\. Bias and RMSE are in cm\-3\.
### D\.4Guidance trade\-off

The guidance scalewwof FM\+CFG controls how strongly the samples follow the inputs, and therefore how wide the ensemble is \(Fig\.[14](https://arxiv.org/html/2610.33557#A4.F14)\)\. Atw=1w=1, FM\+CFG improves on plain FM \(CRPS of 126 cm\-3against 142 cm\-3\) but remains too narrow, with a mean coverage of 0\.53\. Reducingwwto 0\.25 widens the ensemble, raising the coverage to 0\.63 and the spread–skill ratio to 0\.60, at the cost of a higher CRPS \(138 cm\-3\)\. Increasingwwto 2 narrows it, and the coverage drops to 0\.41\. The lowest CRPS of the sweep, 123 cm\-3, occurs atw=0\.75w=0\.75and is practically equal to that atw=0\.5w=0\.5\. Because the sweep was evaluated on the holdout, it was not used to choose the reported operating pointw=1w=1; selectingwwin an application would require separate validation data\.

Figure 14:Guidance scale trades CRPS against interval coverage\.Classifier\-free\-guidance sweep for FM\+CFG on the Hyytiälä temporal holdout\. \(a\) CRPS and correlation; \(b\) mean 80 % interval coverage and spread/skill ratio against the guidance scaleww\.w=1w=1\(dotted line\) is plain conditional sampling;w<1w<1blends towards the unconditional distribution\.

## Appendix EExtended spatial transfer

The paired tests in Sect\.[3\.4](https://arxiv.org/html/2610.33557#S3.SS4)score transferred and locally trained models on the same 60\-day holdout\. As a complementary test over a longer period, Table[9](https://arxiv.org/html/2610.33557#A5.T9)applies each model trained at one station, without retraining, to all valid windows at the other station\. Because these windows cover a different period from the local holdouts, the scores in this appendix should not be compared directly with those in Table[3](https://arxiv.org/html/2610.33557#S3.T3)\. Absolute errors are larger at Kumpula partly because its concentrations are higher\.

In this longer test, diffusion has the highest correlation in both directions \(0\.792 at Kumpula and 0\.689 at Hyytiälä\) and the lowest CRPS and RMSE in the Kumpula\-to\-Hyytiälä direction \(444\.4 cm\-3and 1188\.6 cm\-3, respectively\)\. In the Hyytiälä\-to\-Kumpula direction, the Drifting Model has the lowest CRPS \(1174\.9 cm\-3\) and FM\+CFG the lowest RMSE \(2494\.4 cm\-3\)\. Energy Matching occasionally produces extremely large samples under transfer\. Its RMSE, which squares these values, is therefore not reported in either direction, and in the Kumpula\-to\-Hyytiälä direction its CRPS is also numerically unreliable and is omitted\.

The physical constraints behave as in the local tests \(Table[10](https://arxiv.org/html/2610.33557#A5.T10)\)\. From Hyytiälä to Kumpula, GDE lowers the CRPS of FM from 1204\.7 cm\-3to 1188\.9 cm\-3, and number\+GDE gives the lowest value \(1174\.0 cm\-3\), at a lower correlation than GDE alone \(0\.723 against 0\.732\)\. In the reverse direction, GDE raises the correlation from 0\.660 to 0\.667 and lowers the CRPS from 830\.8 cm\-3to 732\.9 cm\-3; mass\+GDE gives a similar 732\.4 cm\-3\. These improvements are modest and do not remove the component discrepancies described in Appendix[G](https://arxiv.org/html/2610.33557#A7)\.

Figure 15:No Hyytiälä\-trained model reproduces the urban daytime small\-particle enhancement at Kumpula\.Zero\-shot application to SMEAR III for one day, compared with the observation \(top left\)\. Axes and color scale are as in Fig\.[2](https://arxiv.org/html/2610.33557#S2.F2)\.Table 9:Cross\-site \(zero\-shot\) full\-spectrum scores over all valid target\-site windows\. Correlation uses log\-transformed spectra; CRPS and RMSE are in cm\-3\.†\\daggerRMSE is not reported for Energy Matching in either direction because it is dominated by rare, extremely large samples\.‡\\ddaggerIn the Kumpula\-to\-Hyytiälä direction, Energy Matching’s CRPS is also omitted \(see text\)\. Best values are in bold\.\(a\) Hyytiälä\-trained→\\toKumpula

\(b\) Kumpula\-trained→\\toHyytiälä

Table 10:Physical\-constraint ablation under zero\-shot transfer, over all valid target\-site windows\. Columns follow Table[4](https://arxiv.org/html/2610.33557#S3.T4), omitting PSD RMSE; CRPS is in cm\-3and mass RMSE inμ\\mug m\-3\. These periods differ from the local temporal holdouts and the paired tests\. Best values are in bold\.\(a\) Hyytiälä\-trained→\\toKumpula

\(b\) Kumpula\-trained→\\toHyytiälä

## Appendix FDerived quantities and episode verification

The derived quantities of Sect\.[3\.3](https://arxiv.org/html/2610.33557#S3.SS3)are computed separately for every ensemble member before scores are calculated, so the predicted distribution of each quantity retains the dependence between bins\. The ultrafine particle number covers 2\.82–100 nm, and the condensation sink uses Eq\. \([12](https://arxiv.org/html/2610.33557#S2.E12)\) at the reference conditions of Appendix[C\.2](https://arxiv.org/html/2610.33557#A3.SS2)\. Table[F](https://arxiv.org/html/2610.33557#A6)gives the full scores for the four evaluation settings: the two local holdouts and the two transfer directions, all on the same target\-site windows\. The seasonal template and the historical ensemble are always built from the training data of the source station\. A positive CRPSS therefore means an improvement over what the source archive alone would provide, not over a locally trained model; Fig\.[16](https://arxiv.org/html/2610.33557#A6.F16)summarizes skill and coverage\.

Table 11:Derived\-quantity scores on paired target\-site temporal holdouts: observed\-range ultrafine particle \(UFP\) number and the resolved condensation sink\. CRPSS is relative to the source\-training historical ensemble \(“probabilistic template”\); bold marks the best CRPSS among the four learned models in each block\. Deterministic predictions have no interval coverage or width\. Unlike Tables[9](https://arxiv.org/html/2610.33557#A5.T9)and[10](https://arxiv.org/html/2610.33557#A5.T10), transfers here use the same target windows as the local models\.\(a\) Hyytiälä holdout

\(b\) Kumpula holdout

\(c\) Hyytiälä\-trained→\\toKumpula

\(d\) Kumpula\-trained→\\toHyytiälä

Figure 16:Local skill is positive for both quantities, but boreal\-to\-urban transfer loses all skill for ultrafine particle number\.Derived\-quantity skill and calibration on paired target\-site holdouts: ultrafine particle number \(upper group\) and the resolved condensation sink \(lower group\)\. Each group shows CRPSS relative to the source\-training historical ensemble above central\-80 % coverage, with methods on the vertical axis\. Dashed lines mark zero skill and nominal coverage of 0\.80\. Values are from Table[F](https://arxiv.org/html/2610.33557#A6)\.### F\.1Sustained high\-number episodes

An ensemble can also provide probabilities of events\. As an example, we define a sustained high\-number episode for a 24 h sequence by its largest 3 h running mean of ultrafine particle number,QQ, and call the window an event ifQQexceeds a thresholdqq\. The threshold is the 90th percentile ofQQin the training data of the source station and is kept unchanged under transfer\. The predicted event probability is the fraction of ensemble members whoseQQexceedsqq, and its accuracy is measured by the Brier score, the mean squared difference between predicted probability and observed outcome\([Brier, 1950](https://arxiv.org/html/2610.33557#bib.bib4)\)\. To avoid counting the same episode in many overlapping windows, we use one window per day, giving 60 daily cases per setting; consecutive days can still be correlated, and the one\-day bootstrap blocks used for the intervals do not account for longer dependence\. The event is a statistical definition, not a classification of new particle formation and not a health threshold\.

The holdout periods contain few events\. At Hyytiälä \(threshold 4745 cm\-3\) no day qualifies, and at Kumpula \(threshold 13 639 cm\-3\) only 3 of 60 days do\. Under transfer, the Hyytiälä threshold is exceeded on 29 of 60 days at Kumpula, whereas the Kumpula threshold is never exceeded at Hyytiälä, so a percentile threshold does not keep its rarity in a new environment\. Near\-zero Brier scores for samples without events do not show that events can be detected\. Kumpula’s local scores \(0\.031–0\.067\) rest on three events and cannot support a ranking, and in Hyytiälä\-to\-Kumpula transfer, scores of 0\.454–0\.484 are hardly better than always predicting no event \(Brier score29/60=0\.48329/60=0\.483\)\. Table[12](https://arxiv.org/html/2610.33557#A6.T12)and Fig\.[17](https://arxiv.org/html/2610.33557#A6.F17)document the available sample; reliable verification of event probabilities requires longer holdout periods\.

Table 12:Sustained high\-number episode verification, with evaluation settings in rows and methods in columns\. Cells give the Brier score above its 95 % block\-bootstrap interval \(block length 1 day\)\. Event counts give the observed days exceeding the source thresholdqq\(cm\-3\)\.Figure 17:The 60\-day holdouts contain too few high\-number episodes for reliable probability verification\.Top row: reliability diagram per evaluation setting \(mean predicted probability against observed event frequency; marker area scales with the number of days in each bin, and sparse bins are shown as they are\)\. Bottom row: two 24 h examples for the Hyytiälä FM\+GDE model, the highest observed episode \(left\) and the day with the largest mismatch between prediction and outcome, a failure case \(right\)\. The shaded band is the 10–90 % ensemble interval and the dashed line the thresholdqq\.
## Appendix GFitted size components

To see which part of the spectrum is lost under transfer, we describe each observed and predicted spectrum by three log\-normal components\([Hussein et al\., 2005](https://arxiv.org/html/2610.33557#bib.bib13)\),

n⁡\(Dp\)=∑k=13Nk2​π​log10​σk​exp⁡\[−\(log10⁡Dp−log10⁡Dg,k\)22​\(log10⁡σk\)2\],n\(D\_\{p\}\)=\\sum\_\{k=1\}^\{3\}\\frac\{N\_\{k\}\}\{\\sqrt\{2\\pi\}\\,\\log\_\{10\}\\sigma\_\{k\}\}\\exp\\\!\\left\[\-\\frac\{\\left\(\\log\_\{10\}D\_\{p\}\-\\log\_\{10\}D\_\{g,k\}\\right\)^\{2\}\}\{2\\left\(\\log\_\{10\}\\sigma\_\{k\}\\right\)^\{2\}\}\\right\],\(22\)whereNkN\_\{k\}is the number concentration,Dg,kD\_\{g,k\}the geometric mean diameter andσk\\sigma\_\{k\}the geometric standard deviation of componentkk\. The parameters are fitted by minimizing relative residuals, so that the smallest particles carry weight comparable to the peak of the spectrum, and the ranges allowed forDg,kD\_\{g,k\}do not overlap, so that each component keeps its identity across spectra\. We label the components C1–C3 in order of increasing diameter\. They are fitted components, not identified sources or the conventional Aitken, accumulation and coarse modes, and the 2\.82–794 nm grid does not resolve a coarse mode\. The fits describe the spectra well, with median log\-spaceR2R^\{2\}values of 0\.978–0\.991 \(Fig\.[19](https://arxiv.org/html/2610.33557#A7.F19)\), although a good fit does not guarantee that the component parameters are unique\.

For each component, the component numberNkN\_\{k\}is the integral of the fitted log\-normal, including its tails beyond the observed grid\. We summarize the amplitude bias of a model byRkR\_\{k\}, the median over windows of the ratio of predicted to observedNkN\_\{k\}, excluding windows with observedNkN\_\{k\}below 1 cm\-3\. The transfer statistic

Tk=RkcrossRksameT\_\{k\}=\\frac\{R\_\{k\}^\{\\rm cross\}\}\{R\_\{k\}^\{\\rm same\}\}\(23\)divides the ratio at the other station by the ratio at the training station;Tk=1T\_\{k\}=1means that transfer leaves the relative amplitude bias unchanged, not that the transferred prediction is unbiased\.

For Hyytiälä\-trained FM, the local ratiosRkR\_\{k\}are 1\.18, 0\.98 and 0\.99 for C1–C3, and the ratios at Kumpula are 0\.058, 0\.301 and 0\.829\. For Kumpula\-trained FM, the local ratios are 1\.07, 1\.30 and 0\.98, and the ratios at Hyytiälä are 17\.20, 3\.17 and 1\.15\. Across all FM variants, the normalized C1 ratios span 0\.035–0\.054 in the boreal\-to\-urban direction and 13\.8–16\.1 in the reverse direction, whereas C3 stays near unity and C2 lies in between \(Table[13](https://arxiv.org/html/2610.33557#A7.T13), Fig\.[18](https://arxiv.org/html/2610.33557#A7.F18)\)\. Transfer therefore affects the smallest particles most, which is consistent with the weak constraint that bulk mass places on small\-particle number; the fits do not attribute these differences to particular sources\. Locally, the median component ratios of FM are close to one \(Table[14](https://arxiv.org/html/2610.33557#A7.T14)\), but the interquartile range for C1 is broad \(0\.36–3\.89\), so agreement in the median does not imply accurate components in individual windows\.

Table 13:Normalized component transfer ratiosTkT\_\{k\}: cross\-site predicted/observed component number divided by the same model’s ratio at its training site\. C1–C3 are ordered fitted components\. Unity denotes an unchanged relative bias, not necessarily an unbiased target\-site prediction\.\(a\) Hyytiälä\-trained→\\toKumpula

\(b\) Kumpula\-trained→\\toHyytiälä

\\tophlineModelC1C2C3\\middlehlineFlow Matching16\.12\.441\.18FM \+ mass15\.32\.541\.23FM \+ GDE13\.82\.551\.03FM \+ mass \+ GDE15\.52\.491\.10FM \+ number14\.52\.721\.22FM \+ number \+ GDE15\.12\.561\.00\\bottomhline

Figure 18:Boreal\-to\-urban transfer changes the smallest fitted component most\.Component\-number ratios for Hyytiälä\-trained models, evaluated at Hyytiälä \(solid\) and Kumpula \(translucent\)\. Bars show the median predicted/observed number and whiskers the interquartile range\. C1–C3 are ordered fitted components on the 2\.82–794 nm grid; these raw ratios differ from the normalized ratios in Table[13](https://arxiv.org/html/2610.33557#A7.T13)\.Figure 19:Three log\-normal components describe the Hyytiälä holdout spectra well\.Left: distribution of log\-space fitR2R^\{2\}, with the median marked\. Right: an observed spectrum and its ordered fitted components C1–C3 on the 2\.82–794 nm grid\.Table 14:Median ratio of predicted to observed component number at Hyytiälä \(temporal holdout\), with interquartile range\. Unity indicates an unbiased median component amplitude; it does not imply accurate individual spectra\.
## Appendix HSOSAA\-FP and its comparison diagnostic

SOSAA is a one\-dimensional column model that simulates the evolution of gases and aerosol particles in the atmospheric boundary layer by combining meteorology, emissions, gas\-phase chemistry and aerosol dynamics\([Boy et al\., 2011](https://arxiv.org/html/2610.33557#bib.bib3);[Clusius et al\., 2026](https://arxiv.org/html/2610.33557#bib.bib5)\)\. In the SOSAA\-FP configuration of[Clusius et al\. \(2026\)](https://arxiv.org/html/2610.33557#bib.bib5), the column is simulated along backward air\-mass trajectories arriving at Hyytiälä, so that the simulated PSD reflects the emissions and processes encountered upstream of the station\. We use the archived simulations of[Clusius et al\. \(2026\)](https://arxiv.org/html/2610.33557#bib.bib5), whose 60\-bin grid \(1\.07–2000 nm\) contains the 50 diameters used here, so that observed, simulated and generated spectra can be compared bin by bin without interpolation\. The matched windows of Sect\.[3\.5](https://arxiv.org/html/2610.33557#S3.SS5)are those of the interpolation split \(Appendix[A](https://arxiv.org/html/2610.33557#A1)\) that fall within the simulated period\.

Hourly spectra emphasize short\-term differences\. To check whether the disagreement in Sect\.[3\.5](https://arxiv.org/html/2610.33557#S3.SS5)persists at longer time scales, we also apply the validation diagnostic of[Clusius et al\. \(2026\)](https://arxiv.org/html/2610.33557#bib.bib5): 12 d running means and seasonal medians of the number concentration in two size ranges \(Fig\.[20](https://arxiv.org/html/2610.33557#A8.F20)\)\. We use nominal ranges of 28\.2–100 and 100–1000 nm, chosen to match the 30–100 and 100–1000 nm ranges of that study as closely as our bin edges allow\. Both the observed and the simulated series are computed on the 50\-bin analysis grid, so the upper range is truncated at 794 nm in both\.

Figure 20:The validation diagnostic of[Clusius et al\. \(2026\)](https://arxiv.org/html/2610.33557#bib.bib5)applied to the matched Hyytiälä/SOSAA\-FP period\.Twelve\-day running means \(left\) and seasonal medians \(right\) for the nominal 28\.2–100 and 100–1000 nm ranges \(the latter truncated at 794 nm\)\. Solid lines and bars: observations; dashed lines and hatched bars: SOSAA\-FP\.
## Appendix ISampling cost

Iterative generative models trade accuracy against cost through the number of network evaluations \(NFE\) per sample\. For FM, the NFE is the number of Euler stepsKKin Eq\. \([15](https://arxiv.org/html/2610.33557#A2.E15)\); for diffusion, it is the number of denoising steps, which we vary using the deterministic DDIM sampler; the Drifting Model always uses a single evaluation\. Figure[21](https://arxiv.org/html/2610.33557#A9.F21)shows that FM needs about 50 steps for its CRPS to converge, whereas its correlation hardly depends on the step count\. Diffusion reaches its CRPS plateau within a few steps, but its correlation decreases as the number of steps grows\. The two methods therefore trade point accuracy against probabilistic skill along the sampling budget, and neither dominates\. The Drifting Model is not shown because its cost is fixed\.

Figure 21:Flow Matching needs about 50 integration steps, whereas diffusion reaches its CRPS plateau within a few steps at a cost in correlation\.\(a\) CRPS and \(b\) correlation of log\-transformed spectra on the Hyytiälä temporal holdout as a function of the number of function evaluations \(Euler steps for Flow Matching, DDIM denoising steps for diffusion\)\.\\authorcontribution

Z\.\-S\. Liu implemented the modelling, coding, and paper drafting; J\. Köykkä, W\. Peng, A\. Kheder, P\. Clusius and M\. Boy advised on atmospheric interpretation and paper revision\.

\\competinginterests

The authors declare that they have no conflict of interest\.

###### Acknowledgements\.

No funds are involved in this research work\.

## References

- Balcerak et al\. \(2025\)Balcerak, M\., Amiranashvili, T\., Terpin, A\., Shit, S\., Bogensperger, L\., Kaltenbach, S\., Koumoutsakos, P\., and Menze, B\.: Energy Matching: Unifying Flow Matching and Energy\-Based Models for Generative Modeling, in: Advances in Neural Information Processing Systems, 2025\.
- Barahona et al\. \(2026\)Barahona, D\., Breen, K\. H\., Block, K\., and Darmenov, A\.: Deep learning representation of the aerosol size distribution, Geoscientific Model Development, 19, 2437–2459,[10\.5194/gmd\-19\-2437\-2026](https://doi.org/10.5194/gmd-19-2437-2026), 2026\.
- Boy et al\. \(2011\)Boy, M\., Sogachev, A\., Lauros, J\., Zhou, L\., Guenther, A\., and Smolander, S\.: SOSA – a new model to simulate the concentrations of organic vapours and sulphuric acid inside the ABL – Part 1: Model description and initial evaluation, Atmospheric Chemistry and Physics, 11, 43–51,[10\.5194/acp\-11\-43\-2011](https://doi.org/10.5194/acp-11-43-2011), 2011\.
- Brier \(1950\)Brier, G\. W\.: Verification of Forecasts Expressed in Terms of Probability, Monthly Weather Review, 78, 1–3,[10\.1175/1520\-0493\(1950\)078<0001:VOFEIT\>2\.0\.CO;2](https://doi.org/10.1175/1520-0493(1950)078%3C0001:VOFEIT%3E2.0.CO;2), 1950\.
- Clusius et al\. \(2026\)Clusius, P\., Baykara, M\., Xavier, C\., Zhou, P\., Tyree, J\., Foreback, B\., Äijälä, M\., Graeffe, F\., Petäjä, T\., Kulmala, M\., Paasonen, P\., Palmer, P\. I\., and Boy, M\.: Modelling the impact of anthropogenic aerosols on CCN concentrations over a rural boreal forest environment, Atmospheric Chemistry and Physics, 26, 1967–1992,[10\.5194/acp\-26\-1967\-2026](https://doi.org/10.5194/acp-26-1967-2026), 2026\.
- Dal Maso et al\. \(2005\)Dal Maso, M\., Kulmala, M\., Riipinen, I\., Wagner, R\., Hussein, T\., Aalto, P\. P\., and Lehtinen, K\. E\. J\.: Formation and growth of fresh atmospheric aerosols: eight years of aerosol size distribution data from SMEAR II, Hyytiälä, Finland, Boreal Environment Research, 10, 323–336, 2005\.
- Deng et al\. \(2026\)Deng, M\., Li, H\., Li, T\., Du, Y\., and He, K\.: Generative Modeling via Drifting, arXiv preprint arXiv:2602\.04770,[10\.48550/arXiv\.2602\.04770](https://doi.org/10.48550/arXiv.2602.04770), 2026\.
- Gneiting and Raftery \(2007\)Gneiting, T\. and Raftery, A\. E\.: Strictly Proper Scoring Rules, Prediction, and Estimation, Journal of the American Statistical Association, 102, 359–378,[10\.1198/016214506000001437](https://doi.org/10.1198/016214506000001437), 2007\.
- Hari and Kulmala \(2005\)Hari, P\. and Kulmala, M\.: Station for Measuring Ecosystem–Atmosphere Relations \(SMEAR II\), Boreal Environment Research, 10, 315–322, 2005\.
- Hersbach et al\. \(2020\)Hersbach, H\., Bell, B\., Berrisford, P\., Hirahara, S\., Horányi, A\., Muñoz\-Sabater, J\., Nicolas, J\., Peubey, C\., Radu, R\., Schepers, D\., et al\.: The ERA5 global reanalysis, Quarterly Journal of the Royal Meteorological Society, 146, 1999–2049,[10\.1002/qj\.3803](https://doi.org/10.1002/qj.3803), 2020\.
- Ho and Salimans \(2021\)Ho, J\. and Salimans, T\.: Classifier\-Free Diffusion Guidance, in: NeurIPS 2021 Workshop on Deep Generative Models and Downstream Applications,[https://openreview\.net/forum?id=qw8AKxfYbI](https://openreview.net/forum?id=qw8AKxfYbI), 2021\.
- Ho et al\. \(2020\)Ho, J\., Jain, A\., and Abbeel, P\.: Denoising Diffusion Probabilistic Models, in: Advances in Neural Information Processing Systems, vol\. 33, 2020\.
- Hussein et al\. \(2005\)Hussein, T\., Dal Maso, M\., Petäjä, T\., Koponen, I\. K\., Paatero, P\., Aalto, P\. P\., Hämeri, K\., and Kulmala, M\.: Evaluation of an automatic algorithm for fitting the particle number size distributions, Boreal Environment Research, 10, 337–355, 2005\.
- Järvi et al\. \(2009\)Järvi, L\., Hannuniemi, H\., Hussein, T\., Junninen, H\., Aalto, P\. P\., Hillamo, R\., Mäkelä, T\., Keronen, P\., Siivola, E\., Vesala, T\., and Kulmala, M\.: The urban measurement station SMEAR III: Continuous monitoring of air pollution and surface–atmosphere interactions in Helsinki, Finland, Boreal Environment Research, 14, 86–109, 2009\.
- Kulmala et al\. \(2001\)Kulmala, M\., Dal Maso, M\., Mäkelä, J\. M\., Pirjola, L\., Väkevä, M\., Aalto, P\., Miikkulainen, P\., Hämeri, K\., and O’Dowd, C\. D\.: On the formation, growth and composition of nucleation mode particles, Tellus B, 53, 479–490, 2001\.
- Kulmala et al\. \(2012\)Kulmala, M\., Petäjä, T\., Nieminen, T\., Sipilä, M\., Manninen, H\. E\., Lehtipalo, K\., Dal Maso, M\., Aalto, P\. P\., Junninen, H\., Paasonen, P\., et al\.: Measurement of the nucleation of atmospheric aerosol particles, Nature Protocols, 7, 1651–1667,[10\.1038/nprot\.2012\.091](https://doi.org/10.1038/nprot.2012.091), 2012\.
- Lipman et al\. \(2023\)Lipman, Y\., Chen, R\. T\. Q\., Ben\-Hamu, H\., Nickel, M\., and Le, M\.: Flow Matching for Generative Modeling, in: International Conference on Learning Representations \(ICLR\), 2023\.
- Liu et al\. \(2023\)Liu, X\., Gong, C\., and Liu, Q\.: Flow Straight and Fast: Learning to Generate and Transfer Data with Rectified Flow, in: International Conference on Learning Representations \(ICLR\), 2023\.
- Liu et al\. \(2025\)Liu, Z\.\-S\., He, C\., Maier, R\., and Rupp, A\.: PUFM\+\+:PointCloudUpsamplingviaEnhancedFlow Matching, in: arXiv:2501\.15286, 2025\.
- Lovelace et al\. \(2023\)Lovelace, J\., Kishore, V\., Wan, C\., Shekhtman, E\., and Weinberger, K\.: Latent Diffusion for Language Generation, in: Advances in Neural Information Processing Systems, edited by Oh, A\., Naumann, T\., Globerson, A\., Saenko, K\., Hardt, M\., and Levine, S\., vol\. 36, pp\. 56 998–57 025,[10\.52202/075280\-2492](https://doi.org/10.52202/075280-2492), 2023\.
- Petäjä et al\. \(2009\)Petäjä, T\., Mauldin III, R\. L\., Kosciuch, E\., McGrath, J\., Nieminen, T\., Paasonen, P\., Boy, M\., Adamov, A\., Kotiaho, T\., and Kulmala, M\.: Sulfuric acid and OH concentrations in a boreal forest site, Atmospheric Chemistry and Physics, 9, 7435–7448,[10\.5194/acp\-9\-7435\-2009](https://doi.org/10.5194/acp-9-7435-2009), 2009\.
- Rombach et al\. \(2021\)Rombach, R\., Blattmann, A\., Lorenz, D\., Esser, P\., and Ommer, B\.: High\-Resolution Image Synthesis with Latent Diffusion Models, 2021\.
- Seinfeld and Pandis \(2016\)Seinfeld, J\. H\. and Pandis, S\. N\.: Atmospheric Chemistry and Physics: From Air Pollution to Climate Change, John Wiley & Sons, Hoboken, New Jersey, 3rd edn\., 2016\.
- Tong et al\. \(2024\)Tong, A\., Fatras, K\., Malkin, N\., Huguet, G\., Zhang, Y\., Rector\-Brooks, J\., Wolf, G\., and Bengio, Y\.: Improving and generalizing flow\-based generative models with minibatch optimal transport, Transactions on Machine Learning Research, 2024\.
- Vartiainen et al\. \(2025\)Vartiainen, A\., Mikkonen, S\., Leinonen, V\., Petäjä, T\., Wiedensohler, A\., Kühn, T\., and Miinalainen, T\.: Machine learning\-based downscaling of aerosol size distributions from a global climate model, Atmospheric Measurement Techniques, 18, 5763–5782,[10\.5194/amt\-18\-5763\-2025](https://doi.org/10.5194/amt-18-5763-2025), 2025\.
- Vaswani et al\. \(2017\)Vaswani, A\., Shazeer, N\., Parmar, N\., Uszkoreit, J\., Jones, L\., Gomez, A\. N\., Kaiser, Ł\., and Polosukhin, I\.: Attention Is All You Need, in: Advances in Neural Information Processing Systems, vol\. 30,[https://papers\.nips\.cc/paper/7181\-attention\-is\-all\-you\-need](https://papers.nips.cc/paper/7181-attention-is-all-you-need), 2017\.
- Wiedensohler et al\. \(2012\)Wiedensohler, A\., Birmili, W\., Nowak, A\., Sonntag, A\., Weinhold, K\., Merkel, M\., Wehner, B\., Tuch, T\., Pfeifer, S\., Fiebig, M\., et al\.: Mobility particle size spectrometers: harmonization of technical standards and data structure to facilitate high quality long\-term observations of atmospheric particle number size distributions, Atmospheric Measurement Techniques, 5, 657–685,[10\.5194/amt\-5\-657\-2012](https://doi.org/10.5194/amt-5-657-2012), 2012\.
- World Health Organization \(2021\)World Health Organization: WHO global air quality guidelines: particulate matter \(PM2\.5and PM10\), ozone, nitrogen dioxide, sulfur dioxide and carbon monoxide, World Health Organization, Geneva, 2021\.