# Physics-constrained deep learning postprocessing of temperature and humidity

###### Abstract

Weather forecasting centers currently rely on statistical postprocessing
methods to minimize forecast error. This improves skill but can lead to
predictions that violate physical principles or disregard dependencies between
variables, which can be problematic for downstream applications and for the
trustworthiness of postprocessing models, especially when they are based on
new machine learning approaches. Building on recent advances in physics-
informed machine learning, we propose to achieve physical consistency in deep
learning-based postprocessing models by integrating meteorological expertise
in the form of analytic equations. Applied to the post-processing of surface
weather in Switzerland, we find that constraining a neural network to enforce
thermodynamic state equations yields physically-consistent predictions of
temperature and humidity without compromising performance. Our approach is
especially advantageous when data is scarce, and our findings suggest that
incorporating domain expertise into postprocessing models allows the
optimization of weather forecast information while satisfying application-
specific requirements.

Postprocessing is a widely-used approach to reduce forecast error using
statistics, but it may lead to physical inconsistencies. This can be
problematic for trustworthiness and downstream applications. We present the
first machine learning-based postprocessing method intentionally designed to
strictly enforce physical laws. Our framework improves physical consistency
without sacrificing performance, and suggests that human expertise can be
incorporated into postprocessing models via analytic equations.

##  1 Introduction

The use of Weather forecasting centers heavily rely on statistical methods to
correct and refine numerical weather prediction (NWP) outputsis a long-
established effort. Over decades of developments, several kinds of models of
growing complexity have been proposed (vannitsem_statistical_2021), while ,
which improves skill at low computational cost (hemri_trends_2014;
vannitsem_statistical_2021). While the fundamental approach has remained the
same : for decades – statistically relating past NWP model outputs and other
additional data, such as topographic descriptors or seasonality, to
observations . Many weather forecasting centers currently rely on these
statistical methods, as they have consistently shown major skill improvements
(hemri_trends_2014; vannitsem_statistical_2021) for low computational costs
and are actively developing new post-processing tools to be integrated in
their operations– the traditional divide between physical and statistical
modeling is narrowing as increasingly more sophisticated models emerge to
harness the growing volume of available data (vannitsem_statistical_2021).

Current research focuses particularly on machine learning (ML) techniques,
with deep learning (DL) and artificial neural networks (ANNs) emerging as a
modern class of post-processing methods with the potential to outperform
traditional approaches in several aspects. For example, rasp_neural_2018 found
that simple feedforward ANNs can could significantly outperform traditional
regression based post-processing techniques, while being less computationally
demanding at inference time. The authors highlighted that ANNs can could
better incorporate non-linear relationships in a data-driven fashion, and
thanks to their flexibility are were more suited to handle the increasing
volumes of model and observation data thanks to their flexibility. ANNs have
also been combined with other statistical techniques such as Bernstein
polynomials (bremnes_ensemble_2020) for non-parametric probabilistic
predictions.

Furthermore, more sophisticated ANNs such as convolutional neural networks
(CNNs) have the ability to incorporate spatial and temporal data with
unprecedented flexibility. gronquist_deep_2021 used CNNs to improve forecasts
of global weather. hohlein_comparative_2020 and veldkamp_statistical_2021 used
CNNs for spatial downscaling of surface wind fieldfields. A process-specific
application was proposed by chapman_improving_2019;
chapman_probabilistic_2022, with the goal of improving the prediction of
atmospheric rivers, which are filaments of intense horizontal water vapor
transport (ralph_defining_2018). dai_spatially_2021 implemented a generative
adversarial network based on CNNs to produce physically realistic post-
processed forecasts of cloud cover. Thus, first attempts at using DL-based
approaches have shown promising improvements over traditional approaches, as
they better capture non-linear dependencies and often require less feature
engineering. Still, a number of challenges remain in applying ML approaches to
the post-processing world (vannitsem_statistical_2021; haupt_towards_2021),
and they cannot be considered a panacea for all problems. This stresses the
need to include more domain expertise into data-driven approaches in a hybrid
manner, which is facilitated by the availability of custom losses and
architectures in standard machine learning libraries (ebert-uphoff_cira_2021).

When traditional post-processing methods are applied, the ultimate goal is to
minimize the forecast error. This often leads to predictions that do not
exhibit the typical spatial and temporal correlation structure that emerges
from common patterns of atmospheric phenomena, or predictions that violate
physical principles and dependencies between variables. However, for various
applications, such as animated maps of meteorological parameters commonly
disseminated to the public, or in the context of hydrological forecasting
(cloke_ensemble_2009) and renewable energy (pinson_chapter_2018), it is
important to provide forecast scenarios that not only have a smaller error,
but also exhibit realistic spatio-temporal structures (e.g.
schefzik_ensemble_2017, for related work). Furthermore, consistency across
variables should be ensured in various applications. For hydrological
modeling, for example, temperature, radiation, and precipitation should be
consistent at all times.

On the one hand, The issue of consistency is particularly relevant in the
context of probabilistic postprocessing, where sampling from marginal
predictive distributions is an additional step that further breaks the
spatiotemporal and inter-variable consistency. Existing approaches try to
model dependencies from a statistical perspective, and not a physical one. We
believe the two are complementary, closely related yet different as noted by
mohrlen_chapter_2023. In the postprocessing field, limited research is
available in the postprocessing field on the issue of physical consistency, in
the sense of respecting physical principles or variable dependencies . On the
other handbased on analytic relationships. However, this question has recently
gained a lot of attention in the wider ML community, and some applications in
weather modeling are reviewed in kashinath_physics-informed_2021 or
willard_integrating_2022. In general, it has been shown that physical
consistency can be pursued by applying constraints to DL models in order to
prescribe specific physical processes. These constraints can take many forms.
The most widely used approach is to incorporate physics via soft
contraintsconstraints, by defining physics-based losses in addition to common
performance metrics such as mean absolute error (daw_physics-guided_2021).
Another popular approach is to design custom model architectures such that the
physical constraints are strictly enforced (e.g. beucler_enforcing_2021;
ling_reynolds_2016).

In this paper we explore ways to incorporate domain knowledge in DL-based
postprocessing models of temperature and humidity, and the related state
variables. Specifically, we evaluate the effect of imposing constraints based
on the ideal gas law and an empirical approximation of a thermodynamic state
equation, and we identify benefits and disadvantages of different approaches.
The goal of this paper is not to develop a highly optimized model for
operational use, but rather to provide some technical guidelines and insights
about incorporating meteorological expertise, in the form of analytic
equations, in postprocessing models of NWP. For this reason, we simplified the
problem under several aspects as explained in the next section. Most
importantly, we focus here on a deterministic setting, although we hope to
extend this framework to probabilistic predictions in the future.

##  2 Data and methods

###  2.1 Datasets

This studyis conducted to In this study, we make predictions at 128 131
weather station locations covering Switzerland, as shown in Fig. 1. We
consider forecast data from COSMO-E (klasa_evaluation_2018) for the
postprocessing task that spans January 2015-November 2020. 2017-December 2022.
COSMO-E is a limited area weather forecasting model used for the operational
weather forecasts for Switzerland by MeteoSwiss. It is operated as a 21-member
ensemble with a 2.2k̇m horizontal resolution. It runs two cycles per day (0000
and 1200 UTC) with a time horizon of five daysand covers the Alpine domain. We
interpolated . We interpolate the nearest model grid-cell to the selected
station locations, using the $k$-d tree method, and we only retained leadtimes
every 3 hours (. We retain model runs initialized at 0000 UTC and consider
leadtimes between +3 , and +6, +9, +12, etc. ) to reduce redundancy. 24 hours
with hourly steps. We use observational data from the SwissMetNet network
(meteoswiss_automatic_2022) to define our “ground truth”. The predictors and
predictands used in this study are summarised summarized in Table 1. They
represent instantaneous measurements with hourly granularity. We note that
while using both dew point temperature and dew point deficit is redundant,
this is necessary to guarantee a fair comparison between the different model
architectures in Sec. 3. In contrast to prior postprocessing studies, we train
a single model for all leadtimes and use the leadtime as a predictor, thereby
increasing the variance in the training set. This choice was motivated by the
ease of implementation, allowing us to focus on the main aspect of the
proposed methodology, namely the physical constraints.

###  2.2 Cross-validation and random seed

In order to consider the variability due to the data used during optimization,
each model configuration was trained on multiple cross-validation splits. When
dealing with time series, it is important to design the cross-validation split
strategy in a way that ensures independence between sets. At the same time, in
our setup it is desirable that samples of different sets are roughly equally
distributed throughout the year. To do so, we opted for a simple four-folds
cross-validation with a holdout set for testing. Specifically, of our five
years of data we used 20% for testing, 60% for training and 20% for
validation, and we removed values between sets in a gap of five days in order
to ensure independence. A total of four years was considered for training and
validation, meaning each of the four cross-validation folds considered a
different year in the dataset. The dataset partitioning is shown in the
Supplementary Material’s (SM) Fig. SM1. For each cross-validation split, we
applied a standard normalization to the inputs of each set based on the
training set’s mean and standard deviation. Moreover, to account for the
stochasticity of neural networks optimization (due to weights initialization
and gradient descent), each model configuration was trained with 3 different
random seeds. In total, for any given approach, 12 trials were conducted with
the trained models and subsequently all evaluated on the same holdout test
dataset.

Figure 1: The location of the 131 weather stations across Switzerland considered in our study. Stations are colored based on their elevation above sea level. We show the topography of Switzerland and its surroundings in the background. Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL. Table 1: List of predictors and predictands considered in this study. Predictor | Units | Symbol
---|---|---
Air temperature ensemble average | °C | $T$
Dew point temperature ensemble average | °C | $T_{d}$
Dew point deficit ensemble average | °C | $T_{\mathrm{def}}$
Relative humidity ensemble average | % | $RH$
Surface air pressure ensemble average | hPa | $P$
Water vapor mixing ratio ensemble average | g kg-1 | $r$
Forecast leadtime | h |
Sine component of day of the year |  |
Cosine component of day of the year |  |
Sine component of hour of the day |  |
Cosine component of hour of the day |  |
Predictand | Units | Symbol
2 m Air temperature | °C | $T$
2 m Dew point temperature | °C | $T_{d}$
2 m Surface air pressure | hPa | $P$
2 m Relative humidity | % | $RH$
2 m Water vapor mixing ratio | g kg-1 | $r$
|  |

###  2.3 Multi-task neural networks for postprocessing

To keep our DL framework general, the basic building block used for all models
in this study is a fully connected neural network (FCN, see Fig. 2) that takes
as inputs a vector containing the predictors in Table 1 and a vector with
station IDs. Station IDs are mapped into a n-dimensional, real-valued vector
$\textbf{z}\in\mathbb{R}^{4}$ $\textbf{z}\in\mathbb{R}^{6}$ via an embedding
layer and concatenated with the predictors. This approach, here referred to as
”unconstrained setting”, is the same proposed by rasp_neural_2018 and may be
regarded as state-of-the-art for the case of local postprocessing of surface
variables (schulz_machine_2022). It is convenient because it allows to train
for training a single model to do local postprocessing at many stations,
instead of training one model for each station, making its operational
implementation easier. We focus on deterministic forecasts and note that our
$\mathbb{R}^{6}$ embedding has more dimensions than the $\mathbb{R}^{2}$ found
in rasp_neural_2018: this is likely due to the fact that we target five
variables simultaneously and to the complex topography of the Alps (compared
to German stations). The forecasts are deterministic, and we train a model to
predict multiple target variables simultaneously. As such, we are dealing with
a case of multi-task learning (see crawshaw_multi-task_2020, for a survey on
the subject), where the individual errors of each task contribute to the
objective loss function. kendall_multi-task_2017 observed that the relative
weighting of each task’s loss had a strong influence on the overall
performance of such models. This is fairly intuitive in our application
because tasks have have different scales (due to different units) and
uncertainties. The authors proposed a weighting scheme based on the
homoscedastic uncertainty of each task, where the weights are learned during
optimization. To the best of our knowledge, this approach for multi-task
learning is new in the postprocessing field. If jointly postprocessing
multiple meteorological parameters simultaneously becomes more common, it will
be important to design optimal weighting schemes. We use the Mean Squared
Error (MSE) for the loss function $\mathcal{L}_{k}$ of each task $k$. For a
predicted value $\hat{y}_{k,i}$ in physical units (we will use the that
notation for predicted values throughout the rest of the paper) and an
observed value $y_{k,i}$, the task loss $\mathcal{L}_{k}$ is hence defined as:

$\mathcal{L}_{k}\overset{\mathrm{def}}{=}\sum_{i=1}^{N_{\mathrm{samples}}}\frac{\left(y_{k,i}-\hat{y}_{k,i}\right)^{2}}{N_{\mathrm{samples}}},$
(1)

where $N_{\mathrm{samples}}$ is the number of samples. For $p$ tasks, we
define the combined loss $\mathcal{L}$ as:

$\mathcal{L}\overset{\mathrm{def}}{=}\sum_{k=1}^{p}\left[\frac{1}{2}\frac{\mathcal{L}_{k}}{\sigma_{k}^{2}}+\log{\sigma_{k}}\right].$
(2)

Each task’s loss $\mathcal{L}_{k}$ is scaled by the homoscedastic uncertainty
represented by $\sigma^{2}_{k}$, and a regularising term $\log{\sigma_{k}}$ is
added to prevent degenerating towards a zero-weighted loss. In practice, for
improved numerical stability, we learn the log variance
$\eta_{k}\overset{\mathrm{def}}{=}\log{\sigma_{k}^{2}}$ for each task $k$ and
Eq. 2 becomes:

$\mathcal{L}=\sum_{k=1}^{p}\frac{\mathcal{L}_{k}\exp\left(-\eta_{k}\right)+\eta_{k}}{2},$
(3)

where the division by 2 can be ignored for optimization purposes as it does
not influence the minimization objective. To avoid having to learn large
biases, we initialize the bias vector in the model’s output layer using the
training set-averaged output vector, which facilitates optimization. The
optimization is insensitive to the initial values of this bias vector as long
as it has the same order of magnitude as the mean output.

###  2.4 Enforcing analytic constraints in neural networks

The methodology used here follows beucler_enforcing_2021 and was first applied
to neural networks emulating subgrid-scale parametrization for climate
modeling. Conservation laws are enforced during optimization via constraints
in the architecture or the loss function. In this study we aim to enforce
dependencies between variables using the ideal gas law and an an approximate
integral of the Clausius-Clapeyron equation used operationally at MeteoSwiss.
Specifically, we postprocess air temperature ($T$, °C), dew point temperature
($T_{d}$, °C), surface air pressure ($P$, hPa), relative humidity ($RH$, %)
and water vapor mixing ratio ($r$, g kg-1), and . We then aim to enforce the
following constraints:

$\displaystyle RH=f(T,T_{d})=\exp\left(\frac{a\cdot
T_{d}}{b+T_{d}}-\frac{a\cdot T}{b+T}\right),$ (4) $\displaystyle
r=g(P,T_{d})=\frac{622.0\cdot c\cdot\exp\left(\frac{a\cdot
T_{d}}{b+T_{d}}\right)}{P-c\cdot\exp\left(\frac{a\cdot
T_{d}}{b+T_{d}}\right)},$

where $a$, $b$ and $c$ are empirical coefficients, as explained below. The
system of interest includes five variables and two constraints, which leaves
us with three degrees of freedom. The constraints functions $f$ and $g$ are
derived from the following equations:

$e=c\cdot\exp\left(\frac{a\cdot T_{d}}{b+T_{d}}\right)\quad\text{ and }\quad
e_{s}=c\cdot\exp\left(\frac{a\cdot T}{b+T}\right),$ (5)

$RH=\frac{e}{e_{s}}\cdot 100,$ (6)

$r=1000\cdot\frac{0.622\cdot e}{p-e}.$ (7)

We note from Eq. 5 that the parameters of interest are linked by two
additional physical quantities: the water vapor pressure ($e$, hPa) and the
saturation water vapor pressure ($e_{s}$, hPa). Eq. 5 are is structurally
identical to the August–Roche–Magnus equation, an approximate integral of the
Clausius-Clapeyron relation accurate for standard weather conditions. We use
$a=17.368$, $b=238.83$ and $c=6.107$ hPa for $T\geq 0$; $a=17.856$,
$b=245.52$, $c=6.108$ hPa otherwise. We made this choice to ensure consistency
with MeteoSwiss’ internal processing of meteorological variables, but other
values can be found in the literature (e.g. lawrence_relationship_2005). Eq. 7
is a formula for water vapor mixing ratio derived from the ideal gas law for
dry air and water vapor, and can be found in many common textbooks (e.g.
emanuel_atmospheric_1994). We multiply by 1000 to express $r$ in g kg-1 rather
than g g-1.

Figure 2: Summary of the models used in this study. a, the basic building
block of all models, a fully connected network (FCN) preceded by an embedding
layer. b, the unconstrained setting used as a baseline, where all target
variables are predicted directly. c, the architecture constrained setting,
including a physical constraints layer that takes a subset of the target
variables as inputs and returns the complete prediction. d, the loss
constrained neural network, in which physical consistency is enforced by
adding a physics-based penalty $\mathcal{P}$ to the conventional loss
$\mathcal{L}$. e, an offline-constrained neural network, where constraints are
only applied after training using the constraints layer.

We proceed to implement and compare two approaches how to enforce physical
constrains in our networks:

#### architecture-constrained setting:

the constraints are enforced by using a layer that uses equations 4 to derive
$RH$ and $r$ from $T$, $T_{\mathrm{def}}$ and $P$. $T_{\mathrm{def}}$ is the
dew point deficit, to which we apply a Rectified Linear Unit (ReLU) activation
function to ensure positivity before computing dew point temperature as
$T_{d}\overset{\mathrm{def}}{=}T-T_{\mathrm{def}}$. This additional step
allows us to enforce that $T\geq T_{d}$ and $RH\in[0,100]$, two desirable
properties in this case. We show the constrained architecture in Figure 2c.
The trainable part of the model has the same number of layers and units as the
unconstrained architecture, but directly predicts only a subset of variables.
In our case, which has 5 outputs and 2 constraints, we directly predict 3
variables and derive the last 2 via a custom-defined layer encoding equations
4. An important point is that the choice of which variables are computed first
and which are derived analytically is arbitrary. Given $n=5$ variables and
$q=2$ constraints, the total number of possibilities in our case is
$\frac{n!}{q!(n-q)!}=10$. Nevertheless, there are differences in the actual
implementation that point in favor of some configurations. For example, one
may want the analytic constraints arranged in a way such that numerical
stability is not a concern, e.g., avoiding division by zero or asymptotes of
logarithmic functions.

#### loss-constrained setting:

the constraints are enforced by using an additional physics-based loss term
that includes a penalty term $\mathcal{P}$ based on residuals from our set of
analytic equations. As for the architecture-constrained approach, the
variables that we choose to calculate the residuals are arbitrary. Based on
Eq. 4, we define the following constraints:

$\displaystyle\left\\{\begin{array}[]{*{3}{r}l}\mathcal{P}_{RH}\overset{\mathrm{def}}{=}\hat{RH}-f(\hat{T},\hat{T_{d}})=0,\\\
\mathcal{P}_{r}\overset{\mathrm{def}}{=}\hat{r}-g(\hat{P},\hat{T_{d}})=0,\end{array}\right.$
(8)

where physical violations result in non-zero residuals. Using the L2-norm for
consistency with our MSE loss, we formulate the penalty term $\mathcal{P}$
used in the loss function as:

$\mathcal{P}\overset{\mathrm{def}}{=}\frac{(\mathcal{P}_{RH})^{2}}{\sigma^{2}_{RH}}+\frac{(\mathcal{P}_{r})^{2}}{\sigma^{2}_{r}},$
(9)

where we square the residuals $\mathcal{P}_{RH}$ and $\mathcal{P}_{r}$ to
penalize larger violations more, and then scale by the variance of the
observed values $\sigma^{2}_{RH}$ and $\sigma^{2}_{r}$ in order to normalize
the contribution of the two terms. Finally, the physical penalty term is added
to the conventional loss and our training objective becomes minimizing the
physically-constrained loss function $\mathcal{L}_{\mathcal{P}}$:

$\mathcal{L}_{\mathcal{P}}\overset{\mathrm{def}}{=}(1-\alpha)\mathcal{L}+\alpha\mathcal{P},$
(10)

where $\alpha\in\left[0,1\right]$ is a hyperparameter used to scale the
contribution of the physical penalty term. Note that in contrast with the hard
constraints in the architecture, with the soft constraints in the loss, we
have no guarantee that $\mathcal{P}=0$ because stochastic gradient descent
does not generally lead to a zero loss.

Finally, to assess whether enforcing constraints during optimization is
advantageous, we additionally introduce an “offline-constrained” setting (see
Figure 2e) in which constraints are enforced after training. More
specificallySpecifically, we train a model to minimize the MSE of $T$, $T_{d}$
(derived from $T$ and $T_{\mathrm{def}}$) and $P$. $RH$ and $r$ are then
calculated after training so as to exactly enforce our physical constraints.

###  2.5 Libraries, hyperparameters and training

We use the PyTorch deep learning library (paszke_pytorch_2019) to implement
our models, the Ray Tune library (liaw_tune_2018) for hyperparameter tuning,
and Snakemake (molder_sustainable_2021) to manage our workflow 111code at
www.github.com/frazane/pcpp-workflow. For training we use the Adam optimizer
(kingma_adam_2014) with the exponential decay rates set to $\beta_{1}=0.99$
and $\beta_{2}=0.999$ for the first and second moments (kingma_adam_2014),
implement an early-stopping rule based on a validation loss to avoid
overfitting. We use the aggregated normalized mean absolute error (NMAE)
aggregated over all 5 outputs as our validation loss:

$\mathrm{NMAE}\overset{\mathrm{def}}{=}\frac{1}{5}\sum_{k=1}^{5}\sum_{i=1}^{N_{\mathrm{samples}}}\frac{|y_{i,k}-\hat{y_{i,k}}|}{\sigma_{k}},$
(11)

and halt training after 5 epochs in the absence of improvements in the
validation loss. This metric was chosen for the validation and early stopping
because it proved to be more robust to sudden model changes during training,
compared to the training loss. We save the model state with the lowest
validation loss of all training epochs. The hyperparameters used to produce
the main results of this study are shown in Table 2. We chose them after
running a hyperparameter tuning algorithm for the unconstrained model that
considered the aggregated loss of all cross-validation splits. The best
performing hyperparameters configuration was then applied to all models. The
loss-constrained model also required a hyperparameter to scale the influence
of the physics-based penalty term. After testing different values, $\alpha$
was set to 0.995. We discuss this choice in the next section.

Table 2: Hyperparameters used to train the models, along with their optimal value and the search space used for tuning the unconstrained model, represented as the range or set of possible values followed by the sampling method. After selecting the five best configurations automatically, we chose the one with the lowest number of trainable parameters for parsimony. Parameter | Value | Search
---|---|---
Learning rate |  0.0008 0.0007 | [0.01-0.001] $\sim{}$ loguniform
Batch size |  1024 512 | {256, 512, 1024, 2048} $\sim{}$ choice
Units in layer 1 |  128 256 | {32, 64, 128, 256} $\sim{}$ choice
Units in layer 2 | 256 | {32, 64, 128, 256} $\sim{}$ choice
Embedding |  5 6 | {2, 3, 4, 5, 6} $\sim{}$ choice
Patience | 5 epochs |
$\alpha$ (loss-constrained) | 0.995 |
|  |

##  3 Results and discussion

In this section, we present the results of our models when evaluated on unseen
data. We will first compare the performance and physical consistency of
different architectures (Sec. 33.1) before discussing data efficiency (Sec.
33.2) and generalization ability (Sec. 33.3).

###  3.1 Predictive performance and physical consistency

We use two metrics to evaluate the overall performance of our models: the mean
absolute error (MAE) and the Mean Squared Skill Score (MSSS) calculated with
respect to the raw NWP forecast, defined as:

${\color[rgb]{1,0,0}\sout{MSSS}}{\color[rgb]{0,0,1}\uwave{\mathrm{MSSS}}}=1-{\color[rgb]{1,0,0}\sout{\frac{MSE_{PP}}{MSE_{NWP}}}}{\color[rgb]{0,0,1}\uwave{\frac{\mathrm{MSE_{PP}}}{\mathrm{MSE_{NWP}}}}},$
(12)

where $MSE_{PP}$ and $MSE_{NWP}$ $\mathrm{MSE_{PP}}$ and $\mathrm{MSE_{NWP}}$
represent the MSE for our postprocessed forecasts and the NWP forecast
respectively. The MSSS presented here was first computed on each station
individually and then averaged. The overall results are presented in Table 3
and Figure 3a for each variable. For both MAE and MSSS these results show that
the performance is comparable for all NN architectures. Unsurprisingly, the
offline-constrained approach performs poorly on the tasks for which it was not
optimized, that is $RH$ and $r$, and is overall the worst performing model,
with the exception of pressure. Other than that, there There are small
differences in performance of each setting on the different tasks, which we
hypothesize are related to the influence of the constraints coupled with the
multi-task loss weighting. Overall, we observe slightly better results for the
loss and architecture constrained approaches. Table 4(e) also shows the count
of best performing model at shows the results of the Diebold-Mariano
predictive performance tests for the MAE (we followed the implementation of
schulz_machine_2022) applied to each station and for each target variable,
based on the MAE averaged over all trials. These leadtime individually: the
reported values are the percentage of tests for which a method significantly
outperformed another method. Overall, these results appear to reflect those
from Table 3: we don’t observe a clear winner, but the offline constrained
approach appears to be generally worse than the others with the exception of
pressure. While we did not focus on obtaining the best possible performance
for our models, it is useful to know how other widely known postprocessing
approaches perform as reference: in the case of air temperature, we obtained a
MAE of 1.771.63 °C for an altitude-corrected NWP forecast (using a fixed lapse
rate of 6 °C km-1), and found a value of 1.181.01 °C for the MOSMIX product
provided by Deutscher Wetterdienst222https://www.dwd.de/, noting that this
product also includes real-time information from ground stations. In
comparison, our models ’ performance falls around a MAE of 1.5. exhibit MAEs
of approximately 1.35 °C. We note that the MSSS values are surprisingly high,
especially for pressure. These high MSSS values are due to the large errors in
the NWP model. For variables that are strongly tied to elevation, such as
pressure and temperature, the differences in the NWP model elevation and the
true station elevation result in consistently large biases. These elevation
differences can be larger than 100 meters. Taking the example of pressure, the
mean bias at certain stations is almost 90 hPa, which is reduced to almost 0
hPa by the postprocessing models, explaining the high MSSS values.

Table 3: Two performance metrics: The mean-absolute error (MAE, lower values are better) and mean squared skill score (MSSS, higher values are better) for each considered NN architecture (rows) and target variables (columns), averaged over the test set. | $T$ | $T_{d}$ | $P$ | $RH$ | $r$
---|---|---|---|---|---
| MAE | MSSS | MAE | MSSS | MAE | MSSS | MAE | MSSS | MAE | MSSS
Unconstrained |  1.516 1.355 |  0.554 0.351 |  1.657 1.514 |  0.452 0.318 |  1.006 0.529 |  0.998 0.618 |  9.330 8.737 |  0.343 0.290 |  0.589 0.550 |  0.341 0.252
Architecture constrained |  1.510 1.340 |  0.5570.357 |  1.653 1.517 |  0.449 0.305 |  1.005 0.527 |  0.998 0.624 |  9.3018.677 |  0.337 0.286 |  0.586 0.547 |  0.341 0.247
Loss constrained |  1.508 1.343 |  0.556 0.356 |  1.6481.501 |  0.450 0.326 |  1.009 0.530 |  0.998 0.623 |  9.382 8.722 |  0.337 0.297 |  0.5830.539 |  0.3480.274
Offline constrained |  1.517 1.356 |  0.554 0.341 |  1.658 1.530 |  0.451 0.295 |  1.004 0.533 |  0.998 0.607 |  9.516 8.911 |  0.313 0.261 |  0.587 0.551 |  0.340 0.238

Figure 3: a, MAE for each target variable and approach, where the boxplot
distribution represents the 9 trials using different cross-validation splits
and random seeds. Note that the ranges of these distributions are relatively
small compared to the absolute values of the error metric. b, Scatter plot
representing the distribution of physical violations $\mathcal{P}_{RH}$ in RH
units, as a function of $RH$, using all samples of all trials, where points
are color-coded by density. Physical violations are deviations from the zero
line. Table 4: Count For each target variable, the percentage of stations
tests for which an approach (on a row) produces significantly better forecasts
than another (on a column), according to Diebold-Mariano statistical tests
performed with the best MAE(we take . Reported is also on the average over
cross-validation splits and random seeds) percentage of ”wins” or ”losses” for
each approach. Tests are applied to each station and target variableleadtime
individually.

(a) Air temperature |  $T$ Unconstrained |  $T_{d}$ Architecture constrained |  $P$ Loss constrained |  $RH$ Offline constrained |  $r$Winning average
---|---|---|---|---|---
Unconstrained | 19 |  27 6.72 |  28 4.91 |  36 15.67 |  4 9.10
Architecture constrained |  44 32.42 | 20 |  21 16.61 |  55 34.13 |  31 27.72
Loss constrained |  56 26.75 |  79 12.69 | 21 |  38 33.33 |  86 24.26
Offline constrained |  12 14.32 |  5 6.83 |  61 8.98 | 2 |  10 10.04
Losing average | 24.50 | 8.75 | 10.17 | 27.71 |

(b) Dew point temperature | Unconstrained | Architecture constrained | Loss constrained | Offline constrained | Winning average
---|---|---|---|---|---
Unconstrained |  | 12.87 | 6.47 | 22.86 | 14.07
Architecture constrained | 6.18 |  | 3.82 | 17.92 | 9.31
Loss constrained | 26.32 | 33.62 |  | 49.07 | 36.34
Offline constrained | 2.29 | 1.85 | 1.27 |  | 1.81
Losing average | 11.60 | 16.12 | 3.85 | 29.95 |

(c) Surface air pressure | Unconstrained | Architecture constrained | Loss constrained | Offline constrained | Winning average
---|---|---|---|---|---
Unconstrained |  | 28.50 | 30.61 | 32.13 | 30.41
Architecture constrained | 32.82 |  | 36.28 | 38.75 | 35.95
Loss constrained | 32.82 | 28.10 |  | 36.17 | 32.36
Offline constrained | 22.90 | 22.43 | 28.28 |  | 24.54
Losing average | 29.52 | 26.34 | 31.72 | 35.68 |

(d) Relative humidity | Unconstrained | Architecture constrained | Loss constrained | Offline constrained | Winning average
---|---|---|---|---|---
Unconstrained |  | 9.63 | 11.89 | 38.10 | 19.87
Architecture constrained | 21.12 |  | 19.81 | 57.00 | 32.64
Loss constrained | 17.59 | 14.07 |  | 45.91 | 25.86
Offline constrained | 8.65 | 1.82 | 5.74 |  | 5.40
Losing average | 15.79 | 8.51 | 12.48 | 47.00 |

(e) Water vapor mixing ratio | Unconstrained | Architecture constrained | Loss constrained | Offline constrained | Winning average
---|---|---|---|---|---
Unconstrained |  | 7.74 | 1.45 | 10.58 | 6.59
Architecture constrained | 20.39 |  | 4.98 | 19.12 | 14.83
Loss constrained | 57.76 | 32.82 |  | 45.69 | 45.43
Offline constrained | 6.47 | 3.96 | 1.09 |  | 3.84
Losing average | 28.21 | 14.84 | 2.51 | 25.13 |

Figure 3b depicts the physical consistency of the predictions. On the vertical
axis is $\mathcal{P}_{RH}$, which is the difference between the predicted $RH$
and its physically-consistent counterpart derived from the constraint function
$f(T,T_{d})$ (Eq. 4), while we show predicted values on the horizontal axis.
Deviations from zero are therefore considered physical violations. , as well
as values that are not between 0 and 100. Compared to the unconstrained
approach, we observe a significant noticeable decrease of violations using the
loss-constrained approach, although large violations still occur at the ends
of the RH distribution, where values larger than 100% still occur. The
architecture-constrained models bring physical violations to zero to within
machine precision. These results are consistent with beucler_enforcing_2021.
As a side remark, we note that in the unconstrained approach larger violations
tend to occur more at the tails of the distribution, which could indicate that
it is more difficult to converge to a physically consistent solution if the
samples are scarce.

Figure 4: Effect of the hyperparameter $\alpha$ on both the overall physical
violation $\mathcal{P}$ (in red) and the NMAE (in blue) for the test dataset.

In order to choose an optimal value for the $\alpha$ hyperparameter, we tested
several values and compared both the NMAE and the physical consistency of the
predictions. The results are shown in Figure 4 for the following values:
$\alpha\in\\{0.,0.2,0.5,0.8,0.9,0.95,0.99,0.995,0.999,0.9999\\}$$\alpha\in\\{0.,0.2,0.5,0.8,0.9,0.95,0.99,0.995,0.999,0.9999,0.99999\\}$.
We observe that the trade-off between physical constraints and performance is
non-linear: up to $\alpha=0.995$, there is little to no drawback in terms of
performance. In contrast, for higher values of $\alpha$, the NMAE starts to
increase. We note two things: first, the choice of this hyperparameter
$\alpha$ can be chosen based on how much one wishes to prioritize physical
consistency over error reduction. Second, one should be aware that the choice
of the learning rate has a significant influence on $\alpha$’s impact (and
vice-versa) and thus on this trade-off, although this was not further
investigated in this study. We limit ourselves to observe that from a
practical standpoint, this relationship is inconvenient as it makes model
selection harder, and we consider it to be a drawback of the loss-constrained
approach.

###  3.2 Robustness to data scarcity

Figure 5: a, NMAE for each reduction and approach, where the boxplot
distribution represents the 9 trials using different cross-validation splits
and random seeds. As the size of the training dataset reduces, constrained
models perform relatively better. b, Scatter Box plot representing showing the
physical penalty term $\mathcal{P}$’s distribution as a function of physical
violations in training data size for the case of the 0.1%-factor
reductionunconstrained and loss constrained settings. See Figure 3b for
comparisonThe architecture- and offline-constrained approaches have zero
penalty by construction.

Among the potential advantages of constraining neural networks using physical
knowledge are improved robustness to data scarcity. The rationale is that
because we reduce the hypothesis space of the model to a subset of physically-
consistent solutions, we could expect the physics-constrained models to learn
with fewer training samples, or to require fewer parameters. We have therefore
designed an experiment in which we re-trained all models (with fewer
parameters in order to reduce the chance of overfitting, see Table SM2.) on
increasingly reduced training datasets, namely 10%, 1% and 0.120%, 5% and 1%
of the full dataset. The reduction is applied by station, in order to ensure
that all stations are still equally represented in the dataset. For instance,
with the 0.11% reduction we trained with roughly 150 200 samples for each
station. The results, shown in Figure 5a, indicate that enforcing physical
constraints (during model optimization) leads to better predictive performance
when seem to indicate a relatively smaller decrease in performance for the
architecture-constrained approach when the data is scarce. It is also
interesting to note that physical violations under data scarce conditions are
much larger than when using the full dataset, as shown in , although this
difference is rather small. Importantly, the added value of enforcing physical
consistency in data scarce situations is emphasized by Figure 5b. Our results
emphasize the added value of applying physical constraints in such casesFor
the unconstrained model in particular, the physical inconsistencies increase
as we reduce the number of training samples. Conversely, for the architecture
constrained approach, physical inconsistencies are always zero by
construction.

###  3.3 Generalization ability

A common finding of physics-informed ML is that physically-constraining models
could help them generalize to unseeen conditions (willard_integrating_2022).
In order to test the ability of our models to generalize to unseen weather
situations, we design an experiment in which models are trained on a dataset
that excludes the warm season (JJA), and then tested on the warm season only.
This choice was motivated by the increasing relevance of record-shattering
heat extremes in a warming climate. In such situations, the robustness of
postprocessing models is put to test as they have to process and predict
values never seen during training. We present our results in Figure 6, and
observe that physical constraints do not seem to impact the generalization
capabilities of the model to unseen temperature extremes. This result,
consistent with beucler_towards_2020, suggests that the constraints of
Equation 4 are insufficient to guarantee generalization capability for our
mapping.

Figure 6: NMAE for the test dataset containing samples from JJA, conditioned
on different quantiles of the univariate temperature distribution. As
expected, the error increases as temperature are more extremes, but the
relative performance of the considered architectures does not change
significantly.

##  4 Conclusion and outlook

In this study, we have adapted a physically-constrained, deep learning
framework to postprocess weather forecasts, which is new to our knowledge.
More generally, we demonstrated simple ways to integrate scientific expertise,
in the form of analytic equations, into a DL-based postprocessing model.
Compared to unconstrained or loss-constrained models, architecture-constrained
models enforce physical consistency to within machine precision without
degrading the performance of the considered variables. We have also shown that
physical constraints yield better predictions when data are scarce because of
the increased value of physical consistency. However, we did not observe a
significant advantage in terms of generalization capabilities. To interpret
these results, it is useful to distinguish the data efficiency and the
generalization experiment by their underlying challenge, that is,
interpolation and extrapolation, respectively. Physically constraining outputs
can help the model better interpolate data, but cannot mitigate the well-known
limitations of neural networks when it comes to out-of-distribution inputs
(extrapolation).

We believe that a significant value of the proposed methodology lies in its
simplicity: any kind of equation, as long as it is differentiable, can be
included in DL-based postprocessing models. Importantly, this extends beyond
the context of meteorology and physics-based constraints, as we could easily
imagine a similar methodology used to satisfy a diverse set of constraints
defined by the end users. For future research on this topic we foresee (i) an
extension to probabilistic forecasting, e.g. by adopting a generative approach
for the creation of physically consistent ensembles; and (ii) an extension to
a global postprocessing setup, where the model generalizes in space. Finally,
an open question is whether physical constraints have a stronger effect in
more challenging tasks, e.g. with higher-dimensional mappings or more marked
non-linearities.

###### Acknowledgements.

We thank the members of APPP team at MeteoSwiss for helpful comments and
feedback that significantly helped the project. FZ is supported by MeteoSwiss
and ETH Zürich, DN and ML are supported by MeteoSwiss, while TB is supported
by the Canton of Vaud in Switzerland. We also thank the Swiss National
Supercomputing Centre (CSCS) for computing infrastructure. The project’s
GitHub repository is accessible at https://github.com/frazane/pcpp-workflow.
The raw data used to train the models is free for research and education
purposes, and can be accessed via the IDAweb portal at
https://www.meteoswiss.admin.ch/services-and-publications/service/weather-and-
climate-products/data-portal-for-teaching-and-research.html.


# Rapidly Switchable X-ray Angular Momentum from a Free Electron Laser
Oscillator

Nanshun Huang  Zhangjiang Laboratory, Shanghai 201210, China  Haixiao Deng
<EMAIL_ADDRESS>Shanghai Advanced Research Institute, Chinese Academy of
Sciences, Shanghai 201210, China

###### Abstract

X-ray vortices carrying tunable Orbital Angular Momentum (OAM) are an emerging
tool for X-ray characterization technology. However, in contrast to the
generation of vortex beams in the visible wavelength region, the generation of
X-ray vortices in a controlled manner has proved challenging. Here, we
overcome this challenge using an X-ray free-electron laser oscillator (XFELO),
which can produce intense coherent X-rays with switchable OAM. Using the
pinhole mirror in an XFELO, this scheme adjusts only the kinetic energy of the
electron beam to produce vortex beams that can be programmed to dynamically
change between different OAM modes without the need for additional optical
elements. With the nominal parameters of currently constructing 1 MHz
repetition rate facility (i.e. SHINE), the active formation of the OAM modes
of $l=\pm 1$ and $l=\pm 2$ and the rapid switching between them by detuning
the electron beam energy of the XFELO are numerically illustrated. The real-
time switching can be achieved within 200 $\mu$s, while the output pulse
energy can reach the 100 $\mu$J level.

††preprint: APS/123-QED

Spin angular momentum (SAM) and orbital angular momentum (OAM) are two
distinct properties of light. The SAM of light take values only of $+1$ or
$-1$, referring to the left and right circular polarization states
respectively. Besides the SAM, photon beams with spiral wavefronts,
characterized by a transverse phase structure of $\exp(\pm il\phi)$, can carry
the OAM equivalent to $l\hbar$ per photon (where $\phi$ is the azimuthal
angle, and the integer $l$ is known as the vortex order or topological charge)
[1, 2]. Since the OAM can be much greater than the SAM by tuning the
topological charge, OAM light beams have been exploited in various
applications [3], including micromanipulation, imaging, high-capacity
communications, and quantum entanglement. A recent study demonstrated time-
dependent OAM beams in extreme ultraviolet light that result in a self-torque
[4].

Figure 1: General setup for XFELO optical vortex generation and OAM mode
switching. Controlling the electron beam energy together with the ultra-narrow
bandwidth, enables the manipulation of the OAM mode orders

However, the exploitation of X-ray OAM beams is a relatively novel domain. OAM
X-ray beams hold the potential to probe material properties in a way that is
inaccessible using traditional X-ray, such as probing a novel type of phase
dichroism in resonant diffraction from an artificial antiferromagne [5],
inducing charge current loops in fullerenes controlled by tuning the
topological charge [6], and enhancing microscopy [7] and X-ray imaging [8].
Very recently, helicity-dependent spectroscopic signal was experimentally
measured from a hard X-ray OAM beam at the iron (Fe) K edge (7.1 keV, 1.7 Å)
[9]. Overall, interest in optical vortices in the X-ray regime is increasing.

Most OAM-based proposals depend on the difference in the response of matter to
the topological charge $l$ of the vortex [8, 9, 10]. Rapid switching between
OAM modes is highly desirable to fully exploit the benefits of those OAM-based
methods. This requires a high-brightness X-ray source capable of producing
tunable OAM modes. In contrast to the generation of vortex beams at
ultraviolet wavelengths [11, 12], the generation of intense adjustable X-ray
vortices remains challenging and has been extensively investigated. On
timescales of tens of minutes, it was recently demonstrated that the OAM of
the X-ray beam from an artificial spin ice can be switched by increasing the
temperature from 270 K to 380 K or by applying a magnetic field [13]. However,
adjusting the X-ray OAM modes via a programed manner at high switching rates
is a challenge that has yet to be addressed. Considerable attention has been
given to their efficient creation and manipulation.

X-ray free-electron lasers (XFELs) are new generation accelerator-based light
sources with unprecedented brightness at the angstrom wavelength and
femtosecond time scales [14, 15, 16, 17, 18, 19, 20, 21, 22, 23]. With their
unique features, high-brightness FEL X-ray optical vortices can lead to new
research fields. Therefore, several FEL schemes have been studied to generate
vortex beams. In addition to direct conversion approaches [24, 25], these
schemes exploits a helical undulator with harmonic emission or harmonic
modulation [26, 27, 28]. However, to adjust the OAM in these schemes, the
operating harmonics of the helical undulator or transverse mode of the seed
laser must be changed mechanically. To date, no methods has been proposed for
implementing rapid and easy OAM switching using these schemes.

In this work, we propose a method for generating intense X-ray beams carrying
an OAM capable of rapid switching. It is based on an X-ray FEL oscillator
(XFELO) with a pinhole mirror, which consists of an optical cavity (to
circulate X-ray pulses) and an undulator (to provide gain). The XFELO is a
promising candidate for generating intense, fully coherent X-ray pulses at
high repetition rates [29, 30, 31, 32]. The scheme layout is shown in Fig. 1.
In the proposed method, an XFELO can directly generate X-ray vortices whose
topological charge can be adjusted by changing only the electron beam energy,
thereby allowing for very fast OAM switching and programming without changing
the mechanical setup. The underlying principle of modifying the OAM order
through electron beam energy adjustment is rooted in the FEL nature, where the
gain spectrum for various transverse modes experiences a shift due to mode-
dependent Gouy phase detuning, which was initially introduced for the direct
generation of OAM beams [33].

The impact of the Gouy phase on FEL amplification was previously studied [34,
35, 14, 33]. The presence of the Gouy phase causes additional phase slippage
when the radiation field passes through the focus. This additional phase
slippage increases the phase speed of light, and a faster electron velocity is
required to maintain resonance. Therefore, the optimal detuning point is
shifted to higher beam energies [14]. Thus, for certain FEL parameters, it is
possible to detune such that the OAM modes have maximum gain.

This effect can be derived from the change in the FEL resonance condition. The
paraxial beam $E(r,\phi,z)$ can be expanded in a series of Laguerre–Gauss (LG)
modes as follows:

$\mathrm{LG}_{p}^{l}(r,\phi,z)=A_{p}^{|l|}(r,\phi,z)e^{il\phi}e^{i\varphi(z)}$
(1)

with

$\varphi(z)=(1+2p+|l|)\arctan(\frac{z}{z_{\mathrm{R}}}),$

where $\varphi$ is the Gouy phase, $A_{p}^{l}$ is the complex amplitude, $l$
denotes the azimuthal mode number, and $p$ is the radial mode number.
$z_{\mathrm{R}}={\pi w^{2}_{0}}/{\lambda}$ denotes the Rayleigh range with the
wavelength $\lambda$ and beam waist $w_{0}$.

Then, in the basic FEL equations, the relative electron phase associated with
the radiation and the undulator fields is given by $\theta=(k+k_{u})z-\omega
t-\varphi(z)+l\phi$. Effective energy exchange between the electron and
radiation ocurrs if $d\theta/dt=0$. Then, the resonance energy can be obtained
as [36]

$\frac{\gamma_{R}^{2}}{\gamma_{R0}^{2}}=\frac{1}{1-\frac{1}{k_{u}}\frac{d\varphi}{dz}},$
(2)

where $\gamma_{R0}=\sqrt{\lambda_{u}(1+K^{2}/2)/2\lambda}$ denotes the
original FEL resonance energy with the radiation wavelength $\lambda$,
undulator parameter $K$, and undulator period $\lambda_{u}$.
$k_{u}=2\pi/\lambda_{u}$ denotes undulator wavenumber. Owing to the presence
of $z$-dependent phase $\varphi(z)$, the resonance energy for the individual
$l$ modes is different. And the change in resonance energy
$\Delta\gamma_{R}=\gamma_{R}-\gamma_{R0}$ increases with $N=1+2p+|l|$. As a
result, energy detuning can be controlled to maximize the gain of a specific
mode. Through this approach, the OAM modes can be generated and amplified in
an XFELO.

Figure 2: Characterization of the XFELO in the OAM modes. (a), (b), (c), and
(d) show the X-ray OAM beam with $l=1$, $l=-1$, $l=2$, and $l=-2$,
respectively. The numbers indicate the gain curve, amplitude, phase and
spectrum of the OAM beam respectively.

It should be noted that the $l$ and $-l$ modes are not degenerate in XFELO.
This occurs because the interaction between the light and electrons leads to a
micro-bunching structure in a FEL, which is imprinted by the wavefront. Owing
to orthogonality ($\int\exp[i(l-k)\phi]\propto\delta_{l,k}$), the micro-
bunching structures formed by different $l$ modes are not degenerate. This
leads to the robust generation of a single OAM mode. However, it’s worth
mentioning that degenerate states might also appear when the micro-bunching is
not strong enough, specifically during linear gain regions. Moreover, a
detailed analysis can be found in Supplementaries.

It should also be remarked that the XFELO cannot distinguish between a change
in $l$ and a change in $p$ if $N$ doesn’t change. For example, if
$N=1+2p+|l|=3$, the three modes exhibit optimal detuning: $p=1,l=0$ and
$p=0,l=\pm 2$. In particular, undesirable $p>0$ modes are amplified. To
spatially suppress the $p>0$ modes, a pinhole mirror can be used, which
introduces additional loss into these modes. Therefore, by employing a
pinhole, when changing the detuning, only $l$ changes, without exciting
higher-order modes with $p>0$.

As analyzed above, lasing and switching in the OAM modes are realized by
maximizing the gain via careful tuning of the electron beam energy and the use
of a pinhole mirror. Such a small energy change is easily realized by
programming the low-level RF system of the last accelerating module, as
designed in the European XFEL [21]. In addition, the pinhole diamond mirror
has already been proposed as an X-ray regenerative amplifier free-electron
laser [31].

We explored OAM mode switching with three-dimensional numerical simulations
within the context of the Shanghai High-Repetition-Rate XFEL and Extreme Light
Facility (SHINE) [37], which is the first hard X-ray FEL facility in China
(currently under construction). The primary parameters are listed in Table 1.
To demonstrate the robustness of the OAM XFELO, a shot-to-shot jitter was
applied, where the angular fluctuation of the mirrors was set at 15 nrad
(RMS).

Table 1: Main Parameters  Beam Energy | 8 GeV
---|---
Relative Energy Spread | 0.01%
Normalized Emittance | 0.45 mm$\cdot$mrad
Peak Current | 750 A
Undulator Period Length | 26 mm
Undulator Segment Length | 4 m
Photon Energy | 6.96 keV
Mirror Material | Diamond
Bragg Mirror Reflectivity | 85%
Hole Diameter |  $70~{}\mu$m
Bunch Transverse Size |  $\sim 37~{}\mu$m
CRL Focal Length |  $54.7$ m
Figure 3: Top panel: (a) mode decomposition and programed electron beam
energies, (b) power evolution. Right panel: (c) snapshot images of the
transverse intensity and phase distributions, showing OAM state switching
between $|l|=1$ and $|l|=2$. In this case, the switching time is 200 $\mu$s.

## Direct lasing of different OAM modes

To investigate the performance of the XFELO in different OAM modes, electron
beams with different energies were applied (while the XFELO configuration was
fixed). Two beam energies were used. The first beam energy was $\gamma=15680$,
corresponding to $\frac{\delta\gamma}{\gamma_{R0}}=0.055\%$, and was used to
generate the $|l|=1$ mode. The second one was $\gamma=15687$, corresponding to
$\frac{\delta\gamma}{\gamma_{R0}}=0.1\%$, and was used to generate the $|l|=2$
mode. The results are presented in Fig. 2.

Fig. 2 (a) and (b) depicts the cavity output and mode decomposition for the
$|l|=1$ charged vortex case. The pulse energy output from the optical cavity
coupling reaches 100 $\mu$J. Typical transverse intensity and phase
distributions at the exit of the undulator are shown in Fig. 2 (a2, b2) and
Fig. 2 (a3, b3), respectively. The hollow intensity distribution, which is a
characteristic of an optical vortex, can be found. The phase distribution
represents the $l=\pm 1,p=0$ LG mode, confirming the unit charge of the vortex
beam.

Fig. 2 (c) and (d) depicts the XFELO performance in the OAM mode of $|l|=2$.
The pulse energy coupling out from the resonant cavity is approximately 120
$\mu$J at saturation. Because of the relatively low gain, the radiation in the
$|l|=2$ modes requires more round trips to reach saturation. The lower gain is
due to the weaker transverse coupling and the larger emission angle of larger
$|l|$ modes. In addition, the lower single-pass gain in the $|l|=2$ modes
leads to more severe mode competition. Thus, it takes longer to establish the
dominant mode, compared to the $|l|=1$ case. Fig. 2 (c2, d2) and Fig. 2 (c3,
d3) show the The transverse intensity and phase distributions, respectively,
clearly proving that the vortex beam carries a topological charge of $|l|=2$.
In both cases, high-purity optical vortices are obtained. In addition,
distinctive narrow bandwidths are obtained from the these cases, as shown in
right column of Fig. 2.

## Dynamic switching of OAM modes

In addition to the direct lasing of the different vortex beams for macro-pulse
facilities, real-time switching between the OAM modes can be achieved for
continuous-wave facilities. Fig. 3 shows the cases for electron $\gamma$
values from 15676 to 15688, corresponding to OAM mode switching between
$|l|=1$ and $|l|=2$.

The left panel in Fig. 3 shows the (a) mode purity and associated electron
beam energy, (b) power evolution, and (c) snapshots of intensity and phase
distribution. The initial beam energy is $\gamma=15680$, and the XFELO
operates in a vortex state with the $l=1$ mode. As the electron beam energy is
varied to $\gamma=15688$, the OAM mode shifts from $l=1$ to $l=2$, shown in
Fig. 3 (a). Subsequently, the system is switched to the emitted vortex state
of $l=1$ by tuning the electron beam energy to $\gamma=15676$. The snapshots
of intensity and phase distribution in plot (c) of Fig. 3, collected at
different times, clearly show the OAM state switching between $l=1$ and $l=2$.
In this case, the switching time is determined as 200 $\mu$s. These results
also suggest that X-ray OAM can be programmed by the electron beam energy.

It is interesting to observe that during the switch from the $l=1$ mode, the
field is shifted, then causing two vortices with a topological charge of 1
each to merge into a single phase singularity with $l=2$. Conversely, when
switching from $l=2$ to $l=1$, the single phase singularity splits into two
vortices, each with a topological charge of 1. This phenomenon occurs in the
interaction among modes with $|l|>1$ due to the symmetry in phase around the
axis. The dominant mode moves towards the center, causing the fusion and
fragmentation of phase singularities and their relocation onto the axis.
Therefore, the switching of $l$ always commences from the adjacent mode,
implying that the switch could follow paths such as$l=2\rightarrow$ or
$l=-1\rightarrow-2$, and so on. It is crucial to note that during these
switch, the sign does not alter.

In conclusion, we presented a scheme for efficiently generating and
manipulating X-ray vortices in an XFELO. This scheme adjusts only the electron
energy, thereby allowing a fast switching between different OAM modes. To the
best of our knowledge, this is the first study to propose a method for
realizing fast OAM state switching in the X-ray region. Therefore, our
switching rate ( kHz level in the context of SHINE) is expected to be two to
three orders of magnitude faster than those of previous optical approaches.
These findings represent a essential step toward the realization of
programmable X-ray OAM beams. More extensively, the proposed approach holds
promise for at-source solutions for a wide range of applications involving
vortex beams in different areas of fundamental and applied physics, ranging
from super-resolution imaging to material characterization, and from quantum
optics to optical manipulation.

###### Acknowledgements.

We thank Zihan Zhu for providing the start-to-end simulated electron bunch.
This work was supported by the CAS Project for Young Scientists in Basic
Research (YSBR-042), the National Natural Science Foundation of China
(12125508, 11935020),Program of Shanghai Academic/Technology Research Leader
(21XD1404100), and Shanghai Pilot Program for Basic Research - Chinese Academy
of Sciences, Shanghai Branch (JCYJ-SHFY-2021-010). The authors declare no
competing interests. The code used in this study is available from the
corresponding authors upon request.

## References

  * Allen  _et al._ [1992] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, “Orbital angular momentum of light and the transformation of laguerre-gaussian laser modes,” Phys. Rev. A 45, 8185–8189 (1992).
  * Rubinsztein-Dunlop  _et al._ [2016] H. Rubinsztein-Dunlop, A. Forbes, M. V. Berry, M. R. Dennis, D. L. Andrews, M. Mansuripur, C. Denz, C. Alpmann, P. Banzer, T. Bauer, E. Karimi, L. Marrucci, M. Padgett, M. Ritsch-Marte, N. M. Litchinitser, N. P. Bigelow, C. Rosales-Guzmán, A. Belmonte, J. P. Torres, T. W. Neely, M. Baker, R. Gordon, A. B. Stilgoe, J. Romero, A. G. White, R. Fickler, A. E. Willner, G. Xie, B. McMorran, and A. M. Weiner, “Roadmap on structured light,” J. Opt. 19, 013001 (2016).
  * Shen  _et al._ [2019] Y. Shen, X. Wang, Z. Xie, C. Min, X. Fu, Q. Liu, M. Gong, and X. Yuan, “Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities,” Light Sci. Appl. 8, 90 (2019).
  * Rego  _et al._ [2019] L. Rego, K. M. Dorney, N. J. Brooks, Q. L. Nguyen, C.-T. Liao, J. San Román, D. E. Couch, A. Liu, E. Pisanty, M. Lewenstein, L. Plaja, H. C. Kapteyn, M. M. Murnane, and C. Hernández-García, “Generation of extreme-ultraviolet beams with time-varying orbital angular momentum,” Science 364, eaaw9486 (2019).
  * McCarter  _et al._ [2022] M. R. McCarter, A. I. U. Saleheen, A. Singh, R. Tumbleson, J. S. Woods, A. S. Tremsin, A. Scholl, L. E. De Long, J. T. Hastings, S. A. Morley, and S. Roy, “Antiferromagnetic real-space configuration probed by x-ray orbital angular momentum phase dichroism,”  (2022), arXiv:2205.03475 [cond-mat, physics:physics].
  * Wätzel  _et al._ [2016] J. Wätzel, Y. Pavlyukh, A. Schäffer, and J. Berakdar, “Optical vortex driven charge current loop and optomagnetism in fullerenes,” Carbon 99, 439–443 (2016).
  * Kohmura  _et al._ [2020] Y. Kohmura, K. Sawada, M. Mizumaki, K. Ohwada, T. Watanuki, and T. Ishikawa, “X-ray microscope for imaging topological charge and orbital angular momentum distribution formed by chirality,” Opt. Express 28, 24115–24122 (2020).
  * Pancaldi  _et al._ [2023] M. Pancaldi, F. Guzzi, C. S. Bevis, M. Manfredda, J. Barolak, S. Bonetti, I. Bykova, D. De Angelis, G. De Ninno, M. Fanciulli, L. Novinec, E. Pedersoli, A. Ravindran, B. Rösner, C. David, T. Ruchon, A. Simoncig, M. Zangrando, D. E. Adams, P. Vavassori, M. Sacchi, G. Kourousias, G. F. Mancini, and F. Capotondi, “High-resolution ptychographic imaging at a seeded free-electron laser source using OAM beams,”  (2023).
  * Rouxel  _et al._ [2022] J. R. Rouxel, B. Rösner, D. Karpov, C. Bacellar, G. F. Mancini, F. Zinna, D. Kinschel, O. Cannelli, M. Oppermann, C. Svetina, A. Diaz, J. Lacour, C. David, and M. Chergui, “Hard X-ray helical dichroism of disordered molecular media,” Nature Photonics 16, 570–574 (2022).
  * Ye  _et al._ [2019] L. Ye, J. R. Rouxel, S. Asban, B. Rösner, and S. Mukamel, “Probing Molecular Chirality by Orbital-Angular-Momentum-Carrying X-ray Pulses,” J. Chem. Theory Comput. 15, 4180–4186 (2019).
  * Gauthier  _et al._ [2017] D. Gauthier, P. R. Ribič, G. Adhikary, A. Camper, C. Chappuis, R. Cucini, L. F. DiMauro, G. Dovillaire, F. Frassetto, R. Géneaux, P. Miotti, L. Poletto, B. Ressel, C. Spezzani, M. Stupar, T. Ruchon, and G. De Ninno, “Tunable orbital angular momentum in high-harmonic generation,” Nat. Commun. 8, 14971 (2017).
  * Heras  _et al._ [2022] A. d. l. Heras, A. d. l. Heras, A. K. Pandey, A. K. Pandey, J. S. Román, J. Serrano, E. Baynard, G. Dovillaire, M. Pittman, C. G. Durfee, L. Plaja, S. Kazamias, O. Guilbaud, and C. Hernández-García, “Extreme-ultraviolet vector-vortex beams from high harmonic generation,” Optica 9, 71–79 (2022).
  * Woods  _et al._ [2021] J. S. Woods, X. M. Chen, R. V. Chopdekar, B. Farmer, C. Mazzoli, R. Koch, A. S. Tremsin, W. Hu, A. Scholl, S. Kevan, S. Wilkins, W.-K. Kwok, L. E. De Long, S. Roy, and J. T. Hastings, “Switchable X-Ray Orbital Angular Momentum from an Artificial Spin Ice,” Phys. Rev. Lett. 126, 117201 (2021).
  * Pellegrini, Marinelli, and Reiche [2016] C. Pellegrini, A. Marinelli, and S. Reiche, “The physics of x-ray free-electron lasers,” Rev. Mod. Phys. 88, 015006 (2016).
  * Huang  _et al._ [2021] N. Huang, H. Deng, B. Liu, D. Wang, and Z. Zhao, “Features and futures of X-ray free-electron lasers,” The Innovation 2, 100097 (2021).
  * Ackermann  _et al._ [2007] W. Ackermann, G. Asova, V. Ayvazyan, _et al._ , “Operation of a free-electron laser from the extreme ultraviolet to the water window,” Nat. Photonics 1, 336–342 (2007).
  * Emma  _et al._ [2010] P. Emma, R. Akre, J. Arthur, _et al._ , “First lasing and operation of an ångstrom-wavelength free-electron laser,” Nat. Photonics 4, 641–647 (2010).
  * Ishikawa  _et al._ [2012] T. Ishikawa, H. Aoyagi, T. Asaka, _et al._ , “A compact X-ray free-electron laser emitting in the sub-ångström region,” Nat. Photonics 6, 540–544 (2012).
  * Kang  _et al._ [2017] H.-S. Kang, C.-K. Min, H. Heo, _et al._ , “Hard X-ray free-electron laser with femtosecond-scale timing jitter,” Nat. Photonics 11, 708–713 (2017).
  * Prat  _et al._ [2020] E. Prat, R. Abela, M. Aiba, _et al._ , “A compact and cost-effective hard X-ray free-electron laser driven by a high-brightness and low-energy electron beam,” Nat. Photonics 14, 748–754 (2020).
  * Decking  _et al._ [2020] W. Decking, S. Abeghyan, P. Abramian, _et al._ , “A MHz-repetition-rate hard X-ray free-electron laser driven by a superconducting linear accelerator,” Nat. Photonics 14, 391–397 (2020).
  * Allaria  _et al._ [2012] E. Allaria, R. Appio, L. Badano, _et al._ , “Highly coherent and stable pulses from the FERMI seeded free-electron laser in the extreme ultraviolet,” Nat. Photonics 6, 699–704 (2012).
  * Feng  _et al._ [2022] C. Feng, T. Liu, S. Chen, K. Zhou, K. Zhang, Z. Qi, D. Gu, Z. Wang, Z. Jiang, X. Li, B. Wang, X. Wang, W. Zhang, L. Feng, C. Li, T. Lan, B. Li, M. Zhang, H. Deng, D. Xiang, B. Liu, and Z. Zhao, “Coherent and ultrashort soft x-ray pulses from echo-enabled harmonic cascade free-electron lasers,” Optica 9, 785–791 (2022).
  * Rebernik Ribič  _et al._ [2017] P. Rebernik Ribič, B. Rösner, D. Gauthier, E. Allaria, F. Döring, L. Foglia, L. Giannessi, N. Mahne, M. Manfredda, C. Masciovecchio, R. Mincigrucci, N. Mirian, E. Principi, E. Roussel, A. Simoncig, S. Spampinati, C. David, and G. De Ninno, “Extreme-Ultraviolet Vortices from a Free-Electron Laser,” Phys. Rev. X 7, 031036 (2017).
  * Loetgering  _et al._ [2020] L. Loetgering, M. Baluktsian, K. Keskinbora, R. Horstmeyer, T. Wilhein, G. Schütz, K. S. E. Eikema, and S. Witte, “Generation and characterization of focused helical x-ray beams,” Sci. Adv. 6, eaax8836 (2020).
  * Hemsing  _et al._ [2009] E. Hemsing, P. Musumeci, S. Reiche, R. Tikhoplav, A. Marinelli, J. B. Rosenzweig, and A. Gover, “Helical Electron-Beam Microbunching by Harmonic Coupling in a Helical Undulator,” Phys. Rev. Lett. 102, 174801 (2009).
  * Hemsing, Marinelli, and Rosenzweig [2011] E. Hemsing, A. Marinelli, and J. B. Rosenzweig, “Generating Optical Orbital Angular Momentum in a High-Gain Free-Electron Laser at the First Harmonic,” Phys. Rev. Lett. 106, 164803 (2011).
  * Hemsing and Marinelli [2012] E. Hemsing and A. Marinelli, “Echo-Enabled X-Ray Vortex Generation,” Phys. Rev. Lett. 109, 224801 (2012).
  * Kim, Shvyd’ko, and Reiche [2008] K.-J. Kim, Y. Shvyd’ko, and S. Reiche, “A Proposal for an X-Ray Free-Electron Laser Oscillator with an Energy-Recovery Linac,” Phys. Rev. Lett. 100, 244802 (2008).
  * Dai, Deng, and Dai [2012] J. Dai, H. Deng, and Z. Dai, “Proposal for an X-Ray Free Electron Laser Oscillator with Intermediate Energy Electron Beam,” Phys. Rev. Lett. 108, 034802 (2012).
  * Marcus  _et al._ [2020] G. Marcus, A. Halavanau, Z. Huang, J. Krzywinski, J. MacArthur, R. Margraf, T. Raubenheimer, and D. Zhu, “Refractive guide switching a regenerative amplifier free-electron laser for high peak and average power hard x rays,” Phys. Rev. Lett. 125, 254801 (2020).
  * Rauer  _et al._ [2023] P. Rauer, W. Decking, D. Lipka, D. Thoden, T. Wohlenberg, I. Bahns, U. Brueggmann, S. Casalbuoni, M. Di Felice, M. Dommach, J. Grünert, S. Karabekyan, A. Koch, D. La Civita, B. Rio, L. Samoylova, H. Sinn, M. Vannoni, C. Youngman, W. Hillert, and J. Rossbach, “Cavity based x-ray free electron laser demonstrator at the European X-ray Free Electron Laser facility,” Physical Review Accelerators and Beams 26, 020701 (2023).
  * Huang and Deng [2021] N. Huang and H. Deng, “Generating X-rays with orbital angular momentum in a free-electron laser oscillator,” Optica 8, 1020 (2021).
  * Saldin, Schneidmiller, and Yurkov [1993] E. L. Saldin, E. A. Schneidmiller, and M. V. Yurkov, “On a linear theory of an FEL amplifier with an axisymmetric electron beam,” Optics Communications 97, 272–290 (1993).
  * Hemsing  _et al._ [2008] E. Hemsing, A. Marinelli, S. Reiche, and J. Rosenzweig, “Longitudinal dispersion of orbital angular momentum modes in high-gain free-electron lasers,” Phys. Rev. ST Accel. Beams 11, 070704 (2008).
  * Kim, Huang, and Lindberg [2017] K.-J. Kim, Z. Huang, and R. Lindberg,  _Synchrotron Radiation and Free-Electron Lasers: Principles of Coherent X-Ray Generation_ (Cambridge University Press, Cambridge, 2017).
  * Huang  _et al._ [2023] N.-S. Huang, Z.-P. Liu, B.-J. Deng, Z.-H. Zhu, S.-H. Li, T. Liu, Z. Qi, J.-W. Yan, W. Zhang, S.-W. Xiang, Y.-Y. Lei, Y. Zhu, Y.-Z. He, Q.-B. Yuan, F. Gao, R.-B. Deng, S. Sun, Z.-D. Lei, Z.-Q. Jiang, M.-Q. Duan, Y. Zhuan, X.-F. Huang, P.-C. Dong, Z.-L. Li, S.-Y. Si, L. Xue, S. Chen, Y.-F. Liu, Y.-J. Tong, H.-X. Deng, and Z.-T. Zhao, “The MING proposal at SHINE: megahertz cavity enhanced X-ray generation,” Nuclear Science and Techniques 34, 6 (2023).