# Neural Multi-Quantile Forecasting for Optimal Inventory Management

Federico Garza Ramírez
ITAM
<EMAIL_ADDRESS>

###### Abstract

In this work we propose the use of quantile regression and dilated recurrent
neural networks with temporal scaling (MQ-DRNN-s) and apply it to the
inventory management task. This model showed a better performance of up to
3.2% over a statistical benchmark (the quantile autoregressive model with
exogenous variables, QAR-X), being better than the MQ-DRNN without temporal
scaling by 6%. The above on a set of 10,000 time series of sales of El Globo
over a 53-week horizon using rolling windows of 7-day ahead each week.

_K_ eywords Deep Learning $\cdot$ Time Series $\cdot$ Quantile Regression

##  1 Introduction

With the increasing generation of data over time and ever-increasing storage
capacity, companies face the opportunity to make better decisions. Time series
forecasting has proven particularly useful to leverage this information in
inventory management. Inventory management benefits from improved forecasting
and the generation of different scenarios, which translates into less
transportation (reducing costs and CO2 emissions) and better production
planning (less waste), as noted by [34].

In this work, inspired by the recent success of neural networks in the
datasets of the third and fourth Makridakis competitions for time series
forecasting (M3, [28] and M4, [29]) and the Tourism dataset [3], and the
necessity to improve the quality of forecasts for the inventory management
task, we propose a model based on _Dilated Recurrent Neural Networks_ capable
of generating multiple quantiles with adaptive scaling (MQ-DRNN-s) that proved
to be up to 3.1% better than statistical models. These results were obtained
on 10,000 products from one of the largest retail companies in Mexico for 53
weeks.

The work developed by [9] is the main inspiration for the model described in
this document. They propose the _Dilated Recurrent Neural Network_ as a faster
and more accurate alternative to the _Recurrent Neural Network_ , which is
consistent with the results found in this paper. The MQ-DRNN-s explored here
also includes data preprocessing within the model, a fundamental feature for
the proper functioning of machine learning models. This data preprocessing
handled by the model facilitates its deployment to production since less data
engineering is required to obtain reasonably good results.

In addition to this introduction, this document’s structure includes 5
sections as follows. Section 5.2 contains a review on literature on optimal
inventory management and deep learning, Section 8 introduces notation and
describes the MQ-DRNN-s model. Section 8 explores the performance of the model
on the task of peak demand prediction. Finally Section 13 comments conclusions
about this work.

The contribution of this work is framed in the literature of time series
forecasting using machine learning and contributes to three specific aspects.

  1. 1.

Non-parametric quantile estimation: A model is proposed to non-parametrically
forecast multiple quantiles simultaneously from a set of time series. This is
achieved through the quantile loss.

  2. 2.

Temporal scaling technique: The model itself handles the necessary processing
of the data to obtain better results than statistical models. This
characteristic facilitates its use in production environments.

  3. 3.

Forecasting comparative analysis: Finally, we provide a comparison with
statistical models and models similar to ours without the inclusion of our
proposed preprocessing. Our proposed model outperformed classical statistical
benchmarks by 3.2%. And the naively implemented MQ-DRNN without scaling by 6%.

##  2 Notation

In this section, we briefly introduce the notation used throughout the work.
We denote vectors in bold and scalars in non-bold. A time-series vector is
denoted by ${\bm{y}}$, while a particular value for the observation at time
$t$ is denoted by $y_{t}$. The subscript $T$ denotes the last observation in
the time series; that is, $y_{T}$. $H$ gives the prediction horizon;
meanwhile, we denote by $\hat{y}_{T+h}$ the forecasts associated with this
horizon, where $1\leq h\leq H$ and by $\hat{y}^{(q)}_{T+h}$ the forecast of
the quantile $q$ at $T+h$; the actual values are denoted by $y_{T+h}$ and
$y^{(q)}_{T+h}$ respectively. We use the subscript $i$ when it is necessary to
denote different time series, i.e., ${\bm{y}}_{i}$. Finally, we denote the set
of natural numbers less than or equal to $d\in\mathbb{N}$ as $[d]$.

##  3 Optimal Inventory Managment

In time series forecasting, it is standard to estimate the conditional
expectation of the time series in the future. Although it is a widespread
practice, there are scenarios where it is convenient to generate forecasts
above or below the mean value. Quantile forecasting addresses the generation
of such scenarios. Theoretically, the problem of the optimal inventory
allocation has been explored in the _News Vendor Problem_ from operations
research literature [37]. Its solution justifies the choice of a quantile; and
states that maximizing News Vendor profits is given by choosing a quantile.

As described by [37], the newsvendor problem is an excellent framework to work
on inventory-related problems. At its most basic setting, the problem involves
three economic agents: the supplier, the buyer (the newsvendor), and the
consumers, interacting within a single period.

This problem focuses on determining an optimal stock quantity $y^{s}$ of a
single good (newspapers) for the newsvendor to buy from the supplier at a
given cost $v>0$. From the consumers’ side, their demanded quantity is given
by the random variable $Y$ (with density and distribution functions $f_{Y}$
and $F_{Y}$, respectively), at an exogenous price $p>v$. The supplier side
assumes that there are no capacity restrictions and zero transactional costs.
During the single period, the quantity $y^{s}$, a priori chosen by the
newsvendor, might exceed or fall behind the realization of $Y$.

Finally, if the stocked quantity $y^{s}$ exceeds the quantity $y^{d}$ demanded
by the consumers, the remaining units have a unitary recovery cost of $g<v$;
if the opposite happens, then each unit of demand not met has a cost $B>0$ for
the buyer.

Thus, the profits function that reflects the problem the buyer faces is

$\displaystyle\begin{split}\Pi(y^{s},y^{d})&=\left(py^{d}-vy^{s}+g(y^{s}-y^{d})\right)\mathbbm{1}_{\\{y^{d}\leq
y^{s}\\}}\\\
&+\left(py^{s}-vy^{s}-B(y^{d}-y^{s})\right)\mathbbm{1}_{\\{y^{d}>y^{s}\\}}.\end{split}$
(1)

Given that the quantity $y^{s}$ is chosen by the buyer prior to observing
$y^{d}$, the buyer maximizes

$\displaystyle\begin{split}\mathbb{E}[\Pi(y^{s})]=&\int_{0}^{y^{s}}\left(py^{d}-vy^{s}+g(y^{s}-y^{d})\right)f_{Y}(y^{d})dy^{d}\\\
&+\int_{y^{s}}^{\infty}\left(py^{s}-vy^{s}-B(y^{d}-y^{s})\right)f_{Y}(y^{d})dy^{d}.\end{split}$
(2)

Given $\mathbb{E}_{Y}=\mu$, the previous equation reduces to

$\displaystyle\begin{split}\mathbb{E}[\Pi(y^{s})]=&(p-g)\mu-(v-g)y^{s}\\\
&-(p-g+B)\int_{y^{s}}^{\infty}(y^{d}-y^{s})f_{Y}(y^{d})dy^{d}.\end{split}$ (3)

Therefore, the first order conditions needed for an optimal allocation satisfy

$\displaystyle\frac{d}{dy^{s}}\mathbb{E}[\Pi(y^{s})]$
$\displaystyle=(p-v+B)-(p-g+B)F_{Y}(y^{s}),$ (4)
$\displaystyle\frac{d^{2}}{d{y^{s}}^{2}}\mathbb{E}[\Pi(y^{s})]$
$\displaystyle=-(p-g+B)f_{Y}(y^{s})<0.$ (5)

Thus, $\mathbb{E}[\Pi(y^{s})]$ reaches its global maximum at
${y^{s}}^{*}=F_{Y}^{-1}\left(\frac{p-v+B}{p-g+B}\right)$ (Theorem 1, Target
quantile minimizes the quantile loss). The above shows the existence of a
unique optimal quantile $y_{q}^{*}={y^{s}}^{*}$ from the consumers’ demand
density function that maximizes the buyer’s expected profits. Such quantile
$q$ accumulates $\frac{p-v+B}{p-g+B}$ from $F_{Y}$, which represents the ratio
between the costs the buyer faces when exceeding or not the consumers’ demand.

Therefore, companies can search for the quantile that maximizes their profits.
Although there are methods to approximate the value of the quantile based on
unit costs of overestimation and underestimation [21], in this work, we will
assume that this quantile is known. Thus the _News Vendor Problem_ justifies
the existence of an optimal quantile that will serve as input for the time
series forecast for that quantile.

##  4 Deep Learning Models

###  4.1 Universal Approximation Theorem

Probabilistic time series forecasting consists of finding a conditional
probability on covariates for a definite horizon. In practice, finding this
probability function is not feasible. For this purpose, it is necessary to
have a sufficiently flexible family of function approximators with sufficient
accuracy. In particular, neural networks have good properties that allow
estimating this function.

Aside from excellent empirical results, neural networks happen to also have
interesting theoretical guarantees that states the universality of its
approximation. The Universal Approximation Theorem described in the Appendix
(Theorem 2) states that any Lipschitz function (Definition 1) can be
approximated by a 3-layer Neural Network (Definition 3). Although we present
this result, there are versions of the theorem that prove a general version.
See for example [24] and [35].

In particular for _Optimal Inventory Management_ , sufficiently flexible and
powerful approximators will allow to have excellent forecasting accuracy. We
complete this theoretical guarantees with the empirical validation from
Section 8.

###  4.2 Neural Forecasting

Beyond the fully connected neural networks that were proven universal
approximators, in recent years, the use of deep learning for time series
forecasting has become widespread. Since the victory of the hybrid ESRNN model
in the M4 competition [40], fundamental advances in the area have emerged.
This section began by recapping the relevant work starting with Recurrent
Neural Networks and advances made in probabilistic modeling.

_Recurrent Neural Networks_ use hidden states propagated in time (RNN; [15,
47]). The hidden states work as encoders of the previous observations, which
allows the model to be able to memorize relevant information from the series;
thanks to this feature, this model has been used in different contexts of
_Sequence-to-Sequence models_ (Seq2Seq) such as applied to natural language
processing [17, 18] and machine translation [43, 4]. Its widespread use has
allowed an important development in the variants that have emerged; such is
the case of new training techniques and designs that have proven to be
functional as _Long Short Term Memory_ (LSTM; [16]) and the _Gated Recurrent
Units_ (GRU; [10]).

One of the major milestones and precursors of recent advances in sequential
modeling has been developing techniques to learn deeper representations
through longer memory. This advance has been made possible by the use of
convolutions and skip-connections within the recurrent structures. Such is the
case of _Temporal Convolutional Networks_ (TCN; [5]) as well as _Dilated RNN_
(DRNN; [9]), and _WaveNet_ for audio generation and machine translation [44,
12, 22].

As pointed out in [6], the research done in fields such as machine translation
and audio generation, framed in the sequence-to-sequence models literature,
has been transferred to time series modeling with great success. This
achievement is due to the ability of these models to learn complex temporal
dependencies. Within the time series literature, there are models such as the
_Multi-Quantile Convolutional Neural Network_ (MQCNN; [45]), the _Exponential
Smoothing Recurrent Neural Network_ (ESRNN; [40]), or the _Neural Basis
Expansion Analysis_ (NBEATS; [33]) and its exogenous version (NBEATS-x; [32]).

An area of particular interest in time series forecasting with machine
learning techniques has been probabilistic output. In this field, the aim is
to know the conditional distribution of the future values of the time series
and not only its expected value [6]. There are two approaches to the subject.
The first one assumes a parametric form of the output distribution such as the
_Autoregressive Recurrent Networks_ (DeepAR; [39]) or the _Associative and
Recurrent Mixture Density Networks_ (ARMDN; [30]). On the other hand, there is
literature studying nonparametric models; these models produce quantiles
directly by minimizing the quantile loss [23].

Although academia has resisted the widespread adoption of neural networks for
time series forecasting [41], these models have become state-of-the-art in
many areas. For example, in weather forecasting [31], energy markets [25, 13]
and large online retailers [39, 45]. This performance has made its use more
and more constant by the research community [26, 6].

##  5 Quantile Regression

###  5.1 Quantile Loss

The newsvendor problem states that profit maximization is reached at an
optimal quantile of the distribution of the demanded quantity. In business
scenarios where this distribution is unknown, the quantile $q$ can be
heuristically defined. The theory developed by [23] on quantile regression
helps us to generate forecasts for such quantiles. Thus, following the
notation by [45], models can be built to predict the conditional quantile

$y^{(q)}_{T+h}|y_{T},$ (6)

of the target distribution given by

$\mathbb{P}(y_{T+h}\leq y^{(q)}_{T+h}|y_{T})=q,$ (7)

seeking to minimize the quantile loss, given by

$\begin{split}QL_{q}\left(y_{T+h},\hat{y}_{T+h}^{(q)}\right)=&q\max\left(0,y_{T+h}-\hat{y}^{(q)}_{T+h}\right)+\\\
&(1-q)\max\left(0,\hat{y}^{(q)}_{T+h}-y_{T+h}\right),\end{split}$ (8)

where $q$ is the quantile to be forecasted.

As shown in Theorem 1 (Target quantile minimizes the quantile loss), the
_quantile loss_ reaches its minimum at exactly the quantile $q$. Moreover, if
$q=0.5$, the Mean Absolute Error is recovered. Figure 1 visually recovers this
metric; the slope of the plot reflects the desired imbalance by overestimating
or underestimating. The _quantile loss_ over the entire forecast horizon $H$
is calculated arithmetically, averaging the quantile losses of the individual
forecasts.

Figure 1: The Figure shows a plot of the _quantile loss_. As can be seen, this
loss allows to estimate a particular quantile of interest by selecting the
parameter $q$. This quantile gives the slope of the loss.

Finally, it is useful to calculate the sampling probability of the quantile
forecasts as a measure of calibration. That is, the percentage of times that
the actual value is below the quantile forecast. This metric is given by

$CL_{q}({\bm{y}},\hat{{\bm{y}}}^{(q)})=\frac{1}{H}\sum^{H}_{h=1}1_{\\{y_{T+h}\leq\hat{y}_{T+h}^{(q)}\\}}.$
(9)

The forecasts are expected to be well-calibrated, i.e., the difference between
the percentage of times that the actual value is below the quantile forecast
and the quantile is small. Formally, the equation

$|CL_{q}({\bm{y}},\hat{{\bm{y}}}^{(q)})-q|<\varepsilon,$ (10)

is verified for some $\varepsilon>0$ and small. However, as observed by [21]
strategies can be defined where perfect calibration is approximated
asymptotically without information of the process $y$.

###  5.2 Quantile Autoregression

The idea of robust regressions is ancient, as pointed out in [42] potentially
as early as 1760. Subsequently, the framework of quantile regression was
developed primarily in [23]; in their work, they extend the use of least
squares to linearly estimate the conditional mean of $Y$ given $X$ for the
quantile case. Instead of minimizing the mean of the quadratic errors, the
authors develop the idea of minimizing the mean of the quantile loss. [48]
reviews statistical characteristics of the model.

In time series, _quantile autoregression with covariates_ (_QAR-X_) seeks to
model the quantiles through a linear regression where the exogenous variables
are specific lags of the time series. The QAR-X is a local model; it is
defined by

$y_{t}^{(q)}=\beta_{0}+\sum_{l\in L}\beta_{l}y_{t-l}+\mathbb{\alpha}{\bm{x}},$
(11)

for each quantile $q$ and for each time series, where $L$ is a set of lags,
and ${\bm{x}}$ is a vector of exogenous variables, e.g., trend or day of the
week. The forecast can be generated by a _Seasonal Naive_ (C.2) of the fitted
values or in a recursive manner111Using one-step ahead forecasts iteratively
to forecast a horizon greater than one is known as teacher forcing.. It
consists of forecasting the next value and using this recursively to generate
the next one. This method is characterized by accumulating errors [6].
Finally, minimizing the quantile loss finds the best quantile regression
model.

This model has been widely used for problems in different areas. In finance,
it has been used mainly for portfolio construction to minimize the risk,
encoded in a selected quantile $q$ [49, 2, 27]. Similarly, it has been applied
in supply chain logistics for demand forecasting [8] and in electricity-load
to avoid blackouts [14].

##  6 Dilated Recurrent Neural Network

_Recurrent Neural Networks_ (RNN) were first proposed and popularized in the
work of [38, 15, 46]. The fundamental idea of this architecture consists of
connecting the hidden units of the neural network with themselves in a
previous period. In each time step $t+1$, the network receives two inputs:
some external value corresponding to that time step $t+1$ and the output of
the hidden units in the previous time step $t$ [15]. Figure 2 describes the
architecture of this model. All timesteps share the network weights.

Figure 2: The Figure shows the architecture of a Recurrent Neural Network
(RNN). The network receives an external input $x_{t}$ and the output of the
previous hidden unit $h_{t-1}$ to produce the output $y_{t}$.

Although RNNs have shown great success in different applications, one of their
main problems is handling very long sequences [9]. In particular, RNNs present
problems handling very long term dependencies in conjunction with short and
medium-term dependencies; also, the training of this architecture presents
vanishing and exploiting gradient and back-propagation through-time problems;
finally, this back-propagation algorithm increases the training time [36].

We use the _Dilated Recurrent Neural Network_ (DRNN; [9]) model in our method,
shown in Figure 3. In this architecture the hidden unit at $t+1$ no longer
depends on the hidden unit at $t$ but on the hidden unit at $t+1-l$ where
$l>1$ is the dilation size. This dilation size can be different for each of
the hidden layers222Observe that if $l=1$ for each of the hidden layers, the
RNN is recovered.. This architecture alleviates gradient problems, improves
computational time, and better learns long-term dependencies.

Figure 3: The Figure shows the architecture of the _Dilated Recurrent Neural
Network_ (DRNN). For each hidden layer, the size of the dilation changes. The
hidden units of each layer at $t+1$ depend on the hidden units of the same
layer at $t+1-l$ where $l$ is the dilation size.

##  7 Temporal Scaling

One of the main inspirations for the model we propose is the _Exponential
Smoothing-Recurrent Neural Networks_ model (ESRNN; [40]). This model is
hybrid: it scales the time series locally (for each time series) using
_Exponential Smoothing_ (ES; [20]) (C.3), decomposing it into levels and
seasonality, to model the shared trends of all series using the RNN. This
model turned out to be the winner of the fourth version of the Makridakis
competitions (M4; [29]).

Scaling the time series before using the RNN is a crucial part of the ESRNN
algorithm. In our proposal, we take up this idea. We process the time series
locally (for each time series itself) in our model using the median. We then
model the median residuals globally (the whole residuals of the time series
together) with a DRNN. Finally, the model constructs the forecasts by adding
the predicted residuals to the median signal of the time series.

Formally, the time series is given by $y_{t}$. We get the level of each time
series trough its median and it is denoted by $l_{t}$ (12). The difference
between the level and the series is given by $z_{t}$ (13). We use a DRNN to
model the residuals $z_{t}$ jointly. The DRNN receives as arguments the
residuals $z_{t}$, temporal exogenous variables denoted by $x_{t}$, and static
exogenous variables denote by $s_{t}$ (14). The model predicts the level using
the Naive method (C.1), that is, with the last observed value of the same
level (15). Finally, the model computes the forecasts by adding the predicted
level of the series with the outputs of the DRNN (16).

Level: $\displaystyle\quad l_{t}=\text{median}(y_{t})$ (12) Residual:
$\displaystyle\quad z_{t}=y_{t}-l_{t}$ (13) NN Forecast:
$\displaystyle\quad\hat{z}_{t}=\text{DRNN}(z_{t},x_{t},s_{t})$ (14) Level
Forecast: $\displaystyle\quad\hat{l}_{t}=\text{Naive}(l_{t})$ (15) Forecast:
$\displaystyle\quad\hat{y}_{T+h}=\hat{l}_{T+h}+\hat{z}_{T+h}$ (16)

##  8 Multi-Quantile Loss

Machine learning models such as the DRNN we use in our proposed model can
generate multiple outputs at the same time [50]. The above allows us to
generate multiple quantiles for the same model. This ability is achieved
through the use of _multi-quantile loss_

$MQL_{Q}\left({\bm{y}},\left(\hat{{\bm{y}}}^{(q)}\right)_{q\in
Q}\right)=\frac{1}{H|Q|}\sum_{q\in
Q}\sum^{H}_{h=1}QL_{q}\left(y_{T+h},\hat{y}_{T+h}^{(q)}\right),$ (17)

where $H$ is the forecast horizon, and $Q$ is the set of quantiles [45]. To
calculate the _multi-quantile loss_ , the quantile loss (Equation (8)) is
arithmetically averaged for each forecast window and each time series for
multiple horizons and multiple quantiles.

As can be observed, the quantile loss $QL_{q}$ is not comparable between time
series because the error is calculated locally, i.e., for each of the time
series; in particular, this problem arises with product demand data where each
of them may exhibit different levels from each other. Finally, the described
metric is calculated for a set of time series by simply calculating the
arithmetic average of the loss of each time series.

##  9 El Globo Demand Dataset

In order to evaluate our proposed model, we selected a sample of products from
one of the largest companies in Mexico, El Globo. This sample consists of
daily sales of 10,000 products selected for having the highest number of sales
during 2019—the time series range from 2018-04-05 to 2020-02-26. As data
preprocessing, missing values (days without sales) were replaced with zeros to
complete the time series (for a total of 693 observations per product).

Also, the proposed model is capable of receiving exogenous variables and
static variables (for each time series). The day of the week was chosen as the
exogenous variable, while the static variables correspond to the sales center
and the product family. The sales center corresponds to the establishment
where the purchase was made. The product family corresponds to the product
category.

##  10 Training Methodology

We divide the data into training, validation, and test sets, as shown in
Figure 4, to train the proposed model on the dataset described above. The
training set spans from April 5, 2018, to November 28, 2018; the validation
set from November 29, 2018, to February 20, 2019; finally, the test set from
February 21, 2019, to February 26, 2020 (corresponding to 53 weeks). The $.3$,
$.5$, $.7$ and $.9$ quantiles (correspondingly the $30$, $50$, $70$ and $90$
percentiles) were chosen as they were the ones requested by the operation
owning the dataset 333The operations managers were interested in a symmetric
interval around the median covering the 40% of cases and an upper case, the 90
percentile..

The optimal hyperparameters of the model were chosen by observing those that
minimized the _multi-quantile loss_ ($MQL$) in the validation set. The
validation set forecasts were generated with a rolling window each week with a
7-day horizon, considering the task horizon. This technique is known as time
series cross-validation [19]. To carry out the hyperparamter optimization, we
used a Bayesian iterative algorithm (HyperOpt; [7]) using python over the
validation set. Table 1 lists the hyperparameter space over which the
optimization was performed while Table 5 lists the optimal set of
hyperparameters.

Figure 4: Representation of train, validation and test sets. The validation
loss provides the signal for hyperparameter optimization. Validation forecasts
were calculated using 7-day-ahead rolling windows. Early stopping sets were
selected at random (highlighted windows) during the test phase; it consists of
three windows of four weeks each. Forecasts for test set were constructed
using 7-day-ahead rolling windows.

The test set consists of forecasting an entire year, divided into 53 weeks.
Forecasts were generated each week with a 7-day horizon to reflect the
operational performance that many companies follow for forecast generation.
For each week of the test set, the model was retrained with the previous data
and the optimal hyperparameters found using the validation set. We used early
stopping using three randomly selected windows of four weeks each as
regularization. Figure 4 recovers this procedure visually. Finally, the final
forecasts were obtained by ensembling with the median ten instances of the
trained model (parameters were randomly initialized in each instance
controlled by a seed).

Hyperparameter | Description | Considered Values
---|---|---
Architecture Parameters
input_size_multiplier | Input size is a multiple of the output horizon. | [4]
output_size | Output size is the forecast horizon for day ahead forecasting. | [7]
add_nl_layer | Insert a tanh layer between the RNN stack and the output layer. | [True, False]
cell_type | Type of RNN cell. | [GRU, RNN, LSTM, ResLSTM]
dilations | Each list represents one chunk of Dilated LSTMS. | [[[1, 2], [4, 8]], [[1, 2]]]
state_hsize | Dimension of hidden state of the RNN. | range(10, 100)
Optimization and Regularization parameters
learning_rate | Initial learning rate for regression problem. | range(5e-4,1e-3)
lr_decay | The decay constant allows large initial lr to escape local minima. | [0.5]
lr_scheduler_step_size | Number of times the learning rate is halved during train. | [2]
per_series_lr_multip | Multiplier for per-series parameters. | [1]
rnn_weight_decay | Parameter of L2/Tikhonov regularization of the RNN parameters. | [0]
batch_size | The number of samples for each gradient step. | [32]
n_iterations | Maximum number of iterations of gradient descent. | [1,000]
early_stopping | Consecutive iterations without validation loss improvement before stop. | [10]
gradient_clipping_threshold | Max norm of gradient vector. | range(10, 100)
noise_std | Standard deviation of white noise added to input during training. | [0.001]
hyperopt_iters | Number of iterations of hyperparameter search. | [50]
Data Parameters
idx_to_sample_freq | Step size between each window. | [7]
val_to_sample_freq | Step size between each window for validation. | [7]
window_sampling_limit | Number of time windows included in the full dataset. | [all available]

Table 1: Hyperparameter search space for the MQ-DRNN and MQ-DRNN-s models. We
select the hyperparameter set with the best performance in the validation set
using the HyperOpt library.

##  11 Model pipeline

Hyperparameter optimization was performed using a worker with 4 CPUs, 1 GPU
and 15 GB of RAM. Forecast generation was performed in parallel using dask
[11] in python on a cluster of 20 workers with the above configuration raised
with coiled444Coiled is a cloud service developed to automate and scale tasks.
https://coiled.io/. To perform parallel processing, we generated 53
corresponding partitions and uploaded the data to Amazon Web Services.
Parallel processing reduced the computational time by more than 20 times. The
full pipeline is described in Figure (5).

Figure 5: The Figure describes the pipeline for obtaining the forecasts by the
MQ-DRNN-s model. First, the set of 10,000 products $Z$ for which the forecast
will be made is obtained. Then the objective variable $y$, the exogenous
variable $x$ containing the day of the week, the static variables $s$
consisting of the sales center and the product family are obtained.
Subsequently, hyperparameter optimization is performed on the validation set.
With the optimal hyperparameters, the forecasts for each of the 53 weeks are
generated in a rolling fashion.

##  12 Main Results

Table 2 shows the comparison between the _QAR-X_ , _MQ-DRNN_ , and _MQ-DRNN-s_
models on the test dataset. The first column shows the analyzed percentile.
The QAR-X column shows the results associated with the autoregressive quantile
model, while the MQ-DRNN column shows the results associated with the base
model, and finally MQ-DRNN-s reports the results of our proposed model.

As can be seen, the proposed model is better for the evaluation metrics for
the 30, 50, 70, and 90 percentiles (P30, P50, P70, P90). In the case of the 90
percentile, the models have comparable performance, with the QAR-X model being
better by an insignificant difference. The proposed model obtained a
percentage improvement of 3.2% compared to the QAR-X model and a percentage
improvemente of 6% compared to the MQ-DRNN model. Also, the models achieve
adequate calibration for the percentiles analyzed as shown in Table 3.

To test whether the results are statistically significant, we performed a 2-
sided paired t-test following the practice in [33]. The statistical test was
performed for each QAR-X and MQ-DRNN model against our proposed model, MQ-
DRNN-s. Table 4 shows the results. For each of the percentiles and each model,
the difference is statistically significant at 99%.

Forecast Evaluation
---
Percentile | QAR-X | MQ-DRNN | MQ-DRNN-s
P30 | 1.15 | 1.17 | 1.11
P50 | 1.35 | 1.37 | 1.28
P70 | 1.21 | 1.24 | 1.16
P90 | 0.64 | 0.70 | 0.66
Table 2: Quantile forecast accuracy measures for 7-day-ahead demand during 53
weeks for 10,000 time series for models QAR-X, MQ-DRNN and MQ-DRNN-s. The
reported metric is the _quantile loss_ (QL) at P30, P50, P70 and P90
percentile levels. Best results are highlighted in bold. Forecast Calibration
---
Percentile | QAR-X | MQ-DRNN | MQ-DRNN-s
P30 | 30.28 | 29.17 | 30.53
P50 | 49.86 | 49.80 | 49.72
P70 | 69.44 | 71.24 | 70.16
P90 | 89.22 | 91.47 | 90.42
Table 3: Calibration Metric (CL, 9) for 7-day-ahead demand during 53 weeks for
10,000 time series for models QAR-X, MQ-DRNN and MQ-DRNN-s. The reported
metric is the _calibration_ (CL) at P30, P50, P70 and P90 percentile levels.
Qualitative analysis to verify that the forecasts are close to the predicted
quantile. Statistical performance comparison
---
Percentile | QAR-X | MQ-DRNN
P30 | -0.032 | -0.051
(0.00162) | (0.00161)
P50 | -0.069 | -0.091
(0.00173) | (0.00182)
P70 | -0.045 | -0.079
(0.00172) | (0.00193)
P90 | 0.016 | -0.047
(0.00152) | (0.00191)
Table 4: Average quantile loss difference against MQ-DRNN-s, computed as MQ-
DRNN-s - QAR-X and MQ-DRNN-s - MQ-DRNN for P30, P50, P70 and P90 percentile
levels. Standard error of the mean displayed in parenthesis. Bold entries are
significant at the 99% level (2-sided paired t-test).

Figures 6, 7 and 8 show the graphical results of the proposed model’s
forecasts for the whole year. It includes three time series in which we
observe the difference in level for each of the forecasted percentiles and
their ascending order. As an additional feature, we observe that the forecasts
of all quantiles adapt to the behavior of the time series.

Figure 6: The Figure shows the quantile forecasts of the MQ-DRNN-s model for
percentiles P30, P50, P70, and true demand $y$ for 3 sample products from the
dataset. Finally the plot shows 53 sets of forecasts calculated each week.
Figure 7: The Figure shows the quantile forecasts of the MQ-DRNN-s model for
percentiles P30, P50, P70, and true demand $y$ for 3 sample products from the
dataset. Finally the plot shows 53 sets of forecasts calculated each week.
Figure 8: The Figure shows the quantile forecasts of the MQ-DRNN-s model for
percentiles P30, P50, P70, and true demand $y$ for 3 sample products from the
dataset. Finally the plot shows 53 sets of forecasts calculated each week.

##  13 Conclusion

In this work we explore temporal scaling for multi-quantile time series
forecasting using Recurrent Neural Networks (MQ-DRNN-s). This model showed a
better performance of up to 3.1% over an statistical benchmark (the quantile
autoregressive model with exogenous variables, QAR-X) being this one better
than the MQRNN without temporal scaling by 6%. This results were found on a
set of 10,000 time series of sales of one of the most important companies in
Mexico over a 53-week horizon in a rolling fashion of 7-day ahead each week.

For companies seeking to make better decisions based on estimates of the
future, the model proposed here offers three significant benefits: i) on the
one hand, the data engineering required to put it into production is
significantly reduced since an essential part of the data preprocessing occurs
within the model, ii) the generation of multiple quantile scenarios is
satisfied by training a single model; that is, the model does not have to be
trained for each quantile, and finally, iii) the accuracy of the inferences
generated over statistical baselines.

The results are consistent with the existing literature, demonstrating the
remarkable capacity of machine learning models to forecast time series. These
results also have significant business implications as this model’s use would
bring benefits in terms of accuracy and the generation of different scenarios.

##  Appendix A Quantile Loss Optimality

######  Theorem 1.

Let $q\in(0,1)$ the target quantile, $y\in\mathbb{R}$ the target variable and
$\hat{y}\in\mathbb{R}$ the forecast. The quantile loss is defined by

$L_{q}(y,\hat{y})=q\max(y-\hat{y},0)+(1-q)\max(\hat{y}-y,0).$

Then the minimum of the expected quantile loss,
$\hat{y}^{*}=\operatorname*{arg\,min}_{\hat{y}}\mathop{\mathbb{E}_{Y}}[L_{q}(y,\hat{y})]$
is the $q$ quantile, i.e.

$F_{Y}(\hat{y}^{*})=q.$

###### Proof.

We can rewrite the quantile loss as

$L_{q}(y,\hat{y})=q(y-\hat{y})\mathbbm{1}_{y\geq\hat{y}}+(1-q)(\hat{y}-y)\mathbbm{1}_{y<\hat{y}}.$

Then we have

$\mathop{\mathbb{E}_{Y}}[L_{q}(y,\hat{y})]=q\int_{\hat{y}}^{\infty}(y-\hat{y})dF_{Y}(y)+(1-q)\int_{-\infty}^{\hat{y}}(\hat{y}-y)dF_{Y}(y).$

Using the Leibniz integral rule we have,

$\displaystyle\frac{\partial\mathop{\mathbb{E}_{Y}}[L_{q}(y,\hat{y})]}{\partial\hat{y}}$
$\displaystyle=q\left(\int_{\hat{y}}^{\infty}\frac{\partial(y-\hat{y})}{\partial\hat{y}}dF_{Y}(y)-(y-\hat{y})\biggr{\rvert}_{y=\hat{y}}\frac{\partial\hat{y}}{\partial\hat{y}}\right)$
$\displaystyle+(1-q)\left(\int_{-\infty}^{\hat{y}}\frac{\partial(\hat{y}-y)}{\partial\hat{y}}dF_{Y}(y)+(\hat{y}-y)\biggr{\rvert}_{y=\hat{y}}\frac{\partial\hat{y}}{\partial\hat{y}}\right)$
$\displaystyle=-q\int_{\hat{y}}^{\infty}dF_{Y}(y)+(1-q)\int_{-\infty}^{\hat{y}}dF_{Y}(y)$
$\displaystyle=-q(1-F_{Y}(\hat{y}))+(1-q)F_{Y}(\hat{y})$
$\displaystyle=-q+qF_{Y}(\hat{y})+F_{Y}(\hat{y})-qF_{Y}(\hat{y})$
$\displaystyle=F_{Y}(\hat{y})-q.$

Since the quantile loss is a convex function of $q$ the minimum is given by

$\frac{\partial\mathop{\mathbb{E}_{Y}}[L_{q}(y,\hat{y})]}{\partial\hat{y}}=0,$

and therefore,

$F_{Y}(\hat{y}^{*})=q.$

∎

##  Appendix B Optimal Hyperparameters

Description | Hyperparameter | MQRNN | MQRNN-s
---|---|---|---
Architecture | input_size_multiplier | 4 | 4
output_size | 7 | 7
add_nl_layer | False | False
cell_type | ResLSTM | LSTM
dilations | ((1, 2), (4, 8)) | ((1, 2), (4, 8))
state_hsize | 65.0 | 96.0
Optimization | learning_rate | 0.000621 | 0.00085
lr_decay | 0.5 | 0.5
lr_scheduler_step_size | 2 | 2
per_series_lr_multip | 1 | 1
rnn_weight_decay | 0 | 0
batch_size | 32 | 32
n_iterations | 2000 | 2000
early_stopping | True | True
gradient_clipping_threshold | 46.0 | 73.0
noise_std | 0.001 | 0.001
Data | idx_to_sample_freq | 7 | 7
val_idx_to_sample_freq | 7 | 7
window_sampling_limit | 100000 | 100000
Table 5: The Table shows the optimal hyperparameters obtained with the signal
from the validation set. The hyperparameters have been grouped under three
main headings. The model architecture, which corresponds to the specific
instance of the proposed model. The optimization hyperparameters, which is
constituted by the characteristics associated with the optimizer. Finally, the
data sampling parameters. The set of optimal hyperparameters for each model is
shown.

##  Appendix C Basic Models

###  C.1 Naive

The _Naive_ model is the simplest model to forecast time series. It simply
assumes that the futures values will be the last one observed in the time
series. Formally, the forecast for $y_{T+h}$ is given by

$\hat{y}_{T+h|T}=y_{T}\enspace\forall h\in\\{1,...,H\\}.$ (18)

###  C.2 Seasonal Naive

For time series seasonality is a repeated pattern of a fixed size. For example
the same month (for yearly data), the same week of the year (for weekly data)
or the same day (for daily data). The _Sesonal Naive_ model assumes that the
next season value will be equal to the last values observed for that season.
For example, with daily data considering a weekly seasonality the forecast for
the next Monday will be the same as the previous Monday. Formally, using $m$
as the seasonal period the forecast for $y_{T+h}$ is given by

$\hat{y}_{T+h|T}=y_{T+h-m}\enspace\forall h\in\\{1,...,H\\}.$ (19)

###  C.3 Exponential Smoothing

The Exponential Smoothing family of models seeks to describe the trend and
seasonality of the data using weighted averages of past observations. In its
additive version it is described as follows

$\displaystyle\hat{y}_{t+h|t}$
$\displaystyle=l_{t}+\phi_{h}b_{t}+s_{t+h-m(k-1)},$ (20) $\displaystyle l_{t}$
$\displaystyle=\alpha(y_{t}-s_{t-m})+(1-\alpha)(l_{t-1}+\phi b_{t-1}),$ (21)
$\displaystyle b_{t}$ $\displaystyle=\beta(l_{t}-l_{t-1})+(1-\beta)\phi
b_{t-1},$ (22) $\displaystyle s_{t}$ $\displaystyle=\gamma(y_{t}-l_{t-1}-\phi
b_{t-1})+(1-\gamma)s_{t-m},$ (23)

where $l_{t}$ denotes the series at time $t$, $b_{t}$ denotes the slope at
time $t$, $s_{t}$ denotes the seasonal component of the series at time $t$ and
$m$ denotes the season; $\alpha$, $\beta$, $\gamma$ and $\phi$ are smoothing
parameters, $\phi_{h}=\sum_{i=1}^{h}\phi^{i}$ and $k$ is the integer part of
$(h-1)/m$.

The rest of the models in the family are constructed by changing the
assumptions on trend, allowing it to be non-existent, additive or additive
damped; and also on seasonality, allowing it to be non-existent, additive or
multiplicative. For a more detailed explanation see [19].

##  Appendix D Universal Approximation Theorem

######  Definition 1.

Let $f:[0,1]^{d}\rightarrow\mathbb{R}$ with $d\in\mathbb{N}$. We say $f$ is
K-Lipschitz if there exists a constant $K>0$ such that for each
${\bm{x}}_{1},{\bm{x}}_{2}\in[0,1]^{d}$

$|f({\bm{x}}_{1})-f({\bm{x}}_{2})|\leq K||{\bm{x}}_{1}-{\bm{x}}_{2}||_{2}.$

######  Definition 2.

Let $d\in\mathbb{N}$. We say $C\subset\mathbb{R}^{d}$ is a cell if there
exists $(l_{i})_{i\in[d]},(r_{i})_{i\in[d]}\in\mathbb{R}^{d}$ such that

$C=\times_{i=1}^{d}[l_{i},r_{i}].$

######  Lemma 1.

Let $d\in\mathbb{N}$ and $f:[0,1]^{d}\rightarrow\mathbb{R}$ a K-Lipschitz
function. Let $\delta>0$ and $P_{\delta}=\\{C_{1},\dots,C_{N}\\}$ a partition
of $[0,1]^{d}$ where each $C_{i}$ is a cell with side length at most $\delta$.
Let $\\{{\bm{x}}_{1},\dots,{\bm{x}}_{N}\\}$ where ${\bm{x}}_{i}\in C_{i}$ for
each $i\in[N]$ then

  1. 1.

$\sup_{i\in[N]}\sup_{{\bm{x}}^{\prime}\in
C_{i}}|f({\bm{x}}^{\prime})-f({\bm{x}}_{i})|\leq K\delta$.

  2. 2.

The function $g:\mathbb{R}^{d}\rightarrow\mathbb{R}$ given by
$g({\bm{x}})=\sum_{i=1}^{N}1_{\\{{\bm{x}}\in C_{i}\\}}f({\bm{x}}_{i})$
satisfies

$\sup_{{\bm{x}}\in[0,1]^{d}}|f({\bm{x}})-g({\bm{x}})|\leq K\delta.$

###### Proof.

1. Let $i\in[N]$ and ${\bm{x}}^{\prime}\in C_{i}$. Then as ${\bm{x}}^{\prime}\in C_{i}$ we have $||{\bm{x}}^{\prime}-{\bm{x}}_{i}||_{2}\leq\delta$. By Lipschitzness of $f$

$\displaystyle|f({\bm{x}}^{\prime})-f({\bm{x}}_{i})|$ $\displaystyle\leq
K||{\bm{x}}^{\prime}-{\bm{x}}_{i}||_{2}$ $\displaystyle\leq K\delta.$

2. Let ${\bm{x}}\in[0,1]^{d}$. Since $P_{\delta}$ is a partition there exists $i\in[N]$ such that ${\bm{x}}\in C_{i}$. Therefore $1_{\\{{\bm{x}}\in C_{i}\\}}=1$ and $1_{\\{{\bm{x}}\in C_{j}\\}}=0$ for each $i\neq j$. Thus $g({\bm{x}})=f({\bm{x}}_{i})$ and then by 1

$|f({\bm{x}})-g({\bm{x}})|=|f({\bm{x}})-f({\bm{x}}_{i})|\leq K\delta.$

∎

######  Lemma 2.

Let $d\in\mathbb{N}$, $C=\times_{i=1}^{d}[l_{i},r_{i}]\subset\mathbb{R}^{d}$ a
cell and ${\bm{x}}\in\mathbb{R}^{d}$. Then ${\bm{x}}\in C$ if and only if for
each $\gamma\in(0,1)$ it follows that

$\sum_{i=1}^{d}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq r_{i}\\}}\geq
2d-1+\gamma.$

###### Proof.

Necessity. Let ${\bm{x}}\in C$ and $\gamma\in(0,1)$. Then since $C$ is a cell
$l_{i}\leq x_{i}\leq r_{i}$ for each $i\in[d]$. Therefore $1_{\\{x_{i}\geq
l_{i}\\}}=1$ and $1_{\\{x_{i}\leq r_{i}\\}}=1$ for each $i\in[d]$. Thus

$\sum_{i=1}^{d}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq
r_{i}\\}}=\sum_{i=1}^{d}2=2d.$

Since $\gamma\in(0,1)$ it is true that $\gamma-1<0$ and therefore $2d\geq
2d-1+\gamma$. Thus

$\sum_{i=1}^{d}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq r_{i}\\}}=2d\geq
2d-1+\gamma.$

Sufficiency. Let ${\bm{x}}\in\mathbb{R}^{d}$ and suppose that for each
$\gamma\in(0,1)$ it is true that

$\sum_{i=1}^{d}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq r_{i}\\}}\geq
2d-1+\gamma.$

Towards a contradiction suppose ${\bm{x}}\notin C$. Then there exists
$i^{*}\in[d]$ such that $x_{i^{*}}\leq l_{i}$ or $x_{i^{*}}\geq r_{i}$; if
$x_{i^{*}}\leq l_{i^{*}}$ then $1_{\\{x_{i^{*}}\geq l_{i^{*}}\\}}=0$ and
$1_{\\{x_{i^{*}}\leq r_{i^{*}}\\}}=1$. On the other hand if $x_{i^{*}}\geq
r_{i^{*}}$ then $1_{\\{x_{i^{*}}\geq l_{i^{*}}\\}}=1$ and $1_{\\{x_{i^{*}}\leq
r_{i^{*}}\\}}=0$. Therefore $1_{\\{x_{i^{*}}\geq
l_{i^{*}}\\}}+1_{\\{x_{i^{*}}\leq r_{i^{*}}\\}}=1$ and $1_{\\{x_{i}\geq
l_{i}\\}}+1_{\\{x_{i}\leq r_{i}\\}}\leq 2$ for each $i\neq i^{*}$. Thus

$\displaystyle\sum_{i=1}^{d}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq
r_{i}\\}}$ $\displaystyle=1_{\\{x_{i^{*}}\geq
l_{i^{*}}\\}}+1_{\\{x_{i^{*}}\leq r_{i^{*}}\\}}+\sum_{i\neq
i^{*}}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq r_{i}\\}}$ $\displaystyle\leq
1+2(d-1)$ $\displaystyle=2d-1.$

Which is a contradiction since $\gamma\in(0,1)$.

∎

######  Lemma 3.

Let $\tau\geq 0$ and $x\in\mathbb{R}$. Let
$\sigma:\mathbb{R}\rightarrow\mathbb{R}$ given by $\sigma(x)=\max(0,x)$ and
$g_{\tau}:\mathbb{R}\rightarrow\mathbb{R}$ given by $g_{\tau}(x)=|1_{\\{x\geq
0\\}}-(\sigma(\tau x)-\sigma(\tau x-1))|$. Then

$g_{\tau}(x)=\begin{cases}\leq 1&x\in\left[0,\frac{1}{\tau}\right]\\\
0&oc.\end{cases}$

###### Proof.

Let’s use cases.

Case 1. Let $x\in\mathbb{R}$ such that $x<0$. Then $1_{\\{x\geq 0\\}}=0$,
$\sigma(\tau x)=0$ and $\sigma(\tau x-1)=0$. Therefore, $g_{\tau}(x)=0$.

Case 2. Let $x\in\mathbb{R}$ such that $x>\frac{1}{\tau}$. Then $1_{\\{x\geq
0\\}}=1$, $\tau x>1$ and thus $\sigma(\tau x)=\tau x$ and $\sigma(\tau
x-1)=\tau x-1$. Therefore

$g_{\tau}(x)=|1-(\tau x-\tau x+1)|=0.$

Case 3. Let $x\in\mathbb{R}$ such that $x\in\left[0,\frac{1}{\tau}\right]$.
Then $1_{\\{x\geq 0\\}}=1$, $\tau x\in[0,1]$ and thus $\sigma(\tau x)=\tau x$
and $\sigma(\tau x-1)=0$. Therefore

$g_{\tau}(x)=|1-\tau x|=1-\tau x.$

∎

######  Definition 3.

Let $d\in\mathbb{N}$, $L,N\in\mathbb{N}$,
$\theta^{(l)}_{j,i},\beta^{(l)}_{j}\in\mathbb{R}$ for each $i,j\in[N]$,
$l\in[L]$, and $\sigma:\mathbb{R}\rightarrow\mathbb{R}$ given by
$\sigma(x)=\max(0,x)$. A feed-forward Neural Network with $L+1$ layers and $N$
units or neurons per layer is a function
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ given by

$\displaystyle f^{(1)}_{j}$ $\displaystyle=x_{j},$ $\displaystyle f^{(l)}_{j}$
$\displaystyle=\sigma\left(\beta^{(l-1)}_{j}+\sum_{i=1}^{N}\theta^{(l)}_{j,i}f^{(l-1)}_{i}\right),$
$\displaystyle f({\bm{x}})$
$\displaystyle=\sigma\left(\beta^{(L+1)}+\sum_{i=1}^{N}\theta^{(L+1)}_{i}f^{(L)}_{i}\right).$

For $l\in[L]$, $l>1$.

######  Lemma 4.

Let $d\in\mathbb{N}$, $C=\times_{i=1}^{d}[l_{i},r_{i}]\subset\mathbb{R}^{d}$ a
cell. Then there exists a 2-layer Neural Network $\hat{h}^{(2)}_{\tau}$
parametrized by $\tau>0$ such that

$\lim_{\tau\rightarrow\infty}\int_{[0,1]^{d}}|\hat{h}^{(2)}({\bm{x}})-1_{\\{{\bm{x}}\in
C\\}}|d{\bm{x}}=0.$

###### Proof.

Intuitively the proof follows the following steps. At first we define a
function of $\tau$ (a $\sigma$ differential as defined in the Lemma 3) that
approximates the indicator function when $\tau\rightarrow\infty$. Based on
this function we define the neural network to approximate the indicator
function of $C$. Formally, Let $\gamma\in(0,1)$ and
$h:\mathbb{R}^{d}\rightarrow\mathbb{R}$ given by

$h_{C}({\bm{x}})=1_{\\{\sum_{i=1}^{d}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq
r_{i}\\}}\geq 2d-1+\gamma\\}}.$

Then by Lemma 2 we have that $1_{\\{{\bm{x}}\in C\\}}=h_{C}({\bm{x}})$ for all
${\bm{x}}\in C$. Then for $\tau>0$ we define
$\hat{1}^{\tau}:\mathbb{R}\rightarrow\mathbb{R}$ given by

$\hat{1}^{\tau}_{\\{x>=0\\}}=\sigma(\tau x)-\sigma(\tau x-1).$

For $x\in\mathbb{R}$. Observe that for $x,y\in\mathbb{R}$ it follows that
$\hat{1}^{\tau}_{\\{x\geq y\\}}=\hat{1}^{\tau}_{\\{x-y\geq 0\\}}$ and
$\hat{1}_{\\{x\leq y\\}}=\hat{1}_{\\{y-x\geq 0\\}}$. Then we construct
$\hat{h}^{(1)}_{\tau}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ given by

$\hat{h}^{(1)}_{\tau}({\bm{x}})=1_{\\{\sum_{i=1}^{d}\hat{1}^{\tau}_{\\{x_{i}\geq
l_{i}\\}}+\hat{1}^{\tau}_{\\{x_{i}\leq r_{i}\\}}\geq 2d-1+\gamma\\}}.$

And we construct $\hat{h}^{(2)}_{\tau}:\mathbb{R}^{d}\rightarrow\mathbb{R}$
given by

$\hat{h}^{(2)}_{\tau}({\bm{x}})=\hat{1}^{\tau}_{\\{\sum_{i=1}^{d}\hat{1}^{\tau}_{\\{x_{i}\geq
l_{i}\\}}+\hat{1}^{\tau}_{\\{x_{i}\leq r_{i}\\}}\geq 2d-1+\gamma\\}}.$

The function $\hat{h}^{(2)}_{\tau}=\hat{h}^{(2)}$ is our 2-layer Neural
Network. Then by triangle inequality

$\displaystyle|\hat{h}^{(2)}({\bm{x}})-1_{\\{{\bm{x}}\in C\\}}|$
$\displaystyle=|\hat{h}^{(2)}_{\tau}({\bm{x}})-h_{C}({\bm{x}})|$
$\displaystyle\leq|\hat{h}^{(2)}_{\tau}({\bm{x}})-\hat{h}^{(1)}_{\tau}({\bm{x}})|+|\hat{h}^{(1)}_{\tau}({\bm{x}})-h_{C}({\bm{x}})|.$

Let’s define $L_{i}=\left[l_{i},l_{i}+\frac{1}{\tau}\right]$ and
$L=\times_{i=1}^{d}L_{i}$; also
$R_{i}=\left[r_{i}-\frac{1}{\tau},r_{i}\right]$ and $R=\times_{i=1}^{d}R_{i}$.
Then we have that
$|\hat{h}^{(2)}_{\tau}({\bm{x}})-\hat{h}^{(1)}_{\tau}({\bm{x}})|\neq 0$ if and
only if ${\bm{x}}\in[0,1]^{d}\cap(L\cup R)$. Also
$|\hat{h}^{(2)}_{\tau}({\bm{x}})-\hat{h}^{(1)}_{\tau}({\bm{x}})|\leq 1$ and

$\displaystyle\int_{[0,1]\cap L_{i}}dx_{i}$ $\displaystyle\leq\frac{1}{\tau},$
$\displaystyle\int_{[0,1]\cap R_{i}}dx_{i}$ $\displaystyle\leq\frac{1}{\tau}.$

Therefore we have that,

$\displaystyle\int_{[0,1]^{d}}|\hat{h}^{(2)}_{\tau}({\bm{x}})-\hat{h}^{(1)}_{\tau}({\bm{x}})|d{\bm{x}}$
$\displaystyle=\int_{[0,1]^{d}\cap(L\cup
R)}|\hat{h}^{(2)}_{\tau}({\bm{x}})-\hat{h}^{(1)}_{\tau}({\bm{x}})|d{\bm{x}}$
$\displaystyle\leq\int_{[0,1]^{d}\cap(L\cup R)}d{\bm{x}}$
$\displaystyle=\int_{[0,1]^{d}\cap L}d{\bm{x}}+\int_{[0,1]^{d}\cap
R}d{\bm{x}}$ $\displaystyle=\prod_{i=1}^{d}\int_{[0,1]\cap
L_{i}}dx_{i}+\prod_{i=1}^{d}\int_{[0,1]\cap R_{i}}dx_{i}$
$\displaystyle\leq\left(\frac{1}{\tau}\right)^{d}+\left(\frac{1}{\tau}\right)^{d}$
$\displaystyle=2\left(\frac{1}{\tau}\right)^{d}.$

By the other hand, we have that
$|\hat{h}^{(1)}_{\tau}({\bm{x}})-h_{C}({\bm{x}})|\neq 0$ if and only if

$\sum_{i=1}^{d}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq
r_{i}\\}}\neq\sum_{i=1}^{d}\hat{1}^{\tau}_{\\{x_{i}\geq
l_{i}\\}}+\hat{1}^{\tau}_{\\{x_{i}\leq r_{i}\\}}.$

Then we have that,

$\displaystyle 1_{\\{\sum_{i=1}^{d}1_{\\{x_{i}\geq l_{i}\\}}+1_{\\{x_{i}\leq
r_{i}\\}}\neq\sum_{i=1}^{d}\hat{1}^{\tau}_{\\{x_{i}\geq
l_{i}\\}}+\hat{1}^{\tau}_{\\{x_{i}\leq r_{i}\\}}\\}}$
$\displaystyle\leq\sum_{i=1}^{d}1_{\\{1_{\\{x_{i}\geq
l_{i}\\}}\neq\hat{1}^{\tau}_{\\{x_{i}\geq l_{i}\\}}\\}}$
$\displaystyle+\sum_{i=1}^{d}1_{\\{1_{\\{x_{i}\geq
l_{i}\\}}\neq\hat{1}^{\tau}_{\\{x_{i}\geq l_{i}\\}}\\}}$ $\displaystyle\leq
d+d$ $\displaystyle=2d.$

Observe that $1_{\\{x_{i}\geq l_{i}\\}}\neq\hat{1}^{\tau}_{\\{x_{i}\geq
l_{i}\\}}$ for all $x_{i}\in[0,1]\cap L_{i}$ and $1_{\\{x_{i}\leq
r_{i}\\}}\neq\hat{1}^{\tau}_{\\{x_{i}\leq r_{i}\\}}$ for all
$x_{i}\in[0,1]\cap R_{i}$. Therefore

$\displaystyle\int_{[0,1]^{d}}|\hat{h}^{(1)}_{\tau}({\bm{x}})-h_{C}({\bm{x}})|d{\bm{x}}$
$\displaystyle\leq\int_{[0,1]^{d}\cap(L\cup R)}(2d)d{\bm{x}}$
$\displaystyle=2d\left(\int_{[0,1]^{d}\cap L}d{\bm{x}}+\int_{[0,1]^{d}\cap
R}d{\bm{x}}\right)$ $\displaystyle=2d\left(\prod_{i=1}^{d}\int_{[0,1]\cap
L_{i}}dx_{i}+\prod_{i=1}^{d}\int_{[0,1]\cap R_{i}}dx_{i}\right)$
$\displaystyle\leq
2d\left(\left(\frac{1}{\tau}\right)^{d}+\left(\frac{1}{\tau}\right)^{d}\right)$
$\displaystyle=4d\left(\frac{1}{\tau}\right)^{d}.$

Finally, we have that

$\displaystyle\int_{[0,1]^{d}}|\hat{h}^{(2)}({\bm{x}})-1_{\\{{\bm{x}}\in
C\\}}|d{\bm{x}}$
$\displaystyle\leq\int_{[0,1]^{d}}|\hat{h}^{(2)}_{\tau}({\bm{x}})-\hat{h}^{(1)}_{\tau}({\bm{x}})|d{\bm{x}}$
$\displaystyle+\int_{[0,1]^{d}}|\hat{h}^{(1)}_{\tau}({\bm{x}})-h_{C}({\bm{x}})|d{\bm{x}}$
$\displaystyle\leq
2\left(\frac{1}{\tau}\right)^{d}+4d\left(\frac{1}{\tau}\right)^{d}$
$\displaystyle=(2+4d)\left(\frac{1}{\tau}\right)^{d}.$

Therefore

$\lim_{\tau\rightarrow\infty}\int_{[0,1]^{d}}|\hat{h}^{(2)}({\bm{x}})-1_{\\{{\bm{x}}\in
C\\}}|d{\bm{x}}=0.$

∎

######  Theorem 2.

Let $d\in\mathbb{N}$ and $f:[0,1]^{d}\rightarrow\mathbb{R}$ a K-Lipschitz
function. Then for all $\varepsilon>0$ there exists a 3-layer Neural Network
$\hat{f}$ with $O\left(d\right(\frac{L}{\varepsilon}\left){}^{d}\right)$
neurons and ReLU activations such that

$\int_{[0,1]^{d}}|f({\bm{x}})-\hat{f}({\bm{x}})|d{\bm{x}}\leq\varepsilon.$

###### Proof.

Let $\varepsilon>0$, $\delta=\frac{\varepsilon}{K}$ and
$P_{\delta}=\\{C_{1},\dots,C_{N}\\}$ a partition of $[0,1]^{d}$ where each
$C_{i}$ is a cell with side length at most $\delta$. By Lemma 4 there exists a
2-layer Neural Network $\hat{h}^{(2)}$ such that

$\int_{[0,1]^{d}}|\hat{h}^{(2)}({\bm{x}})-1_{\\{{\bm{x}}\in
C\\}}|d{\bm{x}}\rightarrow 0.$

Let’s define $\hat{f}:\mathbb{R}^{d}\rightarrow\mathbb{R}$ given by

$\hat{f}({\bm{x}})=\sum_{i=1}^{N}\hat{h}^{(2)}({\bm{x}})f({\bm{x}}).$

This is our 3-layer Neural Network. Observe that

$\hat{f}({\bm{x}})\rightarrow\sum_{i=1}^{N}1_{\\{{\bm{x}}\in
C\\}}f({\bm{x}}).$

Since $\hat{h}^{(2)}$ can approximate arbitrarily well $1_{\\{{\bm{x}}\in
C\\}}$ we have by Lemma 1

$\sup_{{\bm{x}}\in[0,1]^{d}}|f({\bm{x}})-\hat{f}({\bm{x}})|\leq K\delta.$

Hence

$\displaystyle\int_{[0,1]^{d}}|f({\bm{x}})-\hat{f}({\bm{x}})|d{\bm{x}}$
$\displaystyle\leq\int_{[0,1]^{d}}K\delta d{\bm{x}}$ $\displaystyle=K\delta$
$\displaystyle=K\left(\frac{\varepsilon}{K}\right)$
$\displaystyle=\varepsilon.$

∎

## References

  * [1]
  * Allen et al. [2011] Allen, D., Powell, R. and Singh, A. [2011], ‘Quantile regression as a tool for portfolio investment decisions during times of financial distress’, Annals of Financial Economics 06, 1150003.
  * Athanasopoulos et al. [2011] Athanasopoulos, G., Hyndman, R. J., Song, H. and Wu, D. C. [2011], ‘The tourism forecasting competition’, International Journal of Forecasting 27(3), 822–844.  Special Section 1: Forecasting with Artificial Neural Networks and Computational Intelligence Special Section 2: Tourism Forecasting.
https://www.sciencedirect.com/science/article/pii/S016920701000107X

  * Bahdanau et al. [2015] Bahdanau, D., Cho, K. and Bengio, Y. [2015], ‘Neural machine translation by jointly learning to align and translate’, CoRR abs/1409.0473.
  * Bai et al. [2018] Bai, S., Kolter, J. Z. and Koltun, V. [2018], ‘An empirical evaluation of generic convolutional and recurrent networks for sequence modeling’, CoRR abs/1803.01271.
http://arxiv.org/abs/1803.01271

  * Benidis et al. [2020] Benidis, K., Rangapuram, S. S., Flunkert, V., Wang, B., Maddix, D. C., Türkmen, A. C., Gasthaus, J., Bohlke-Schneider, M., Salinas, D., Stella, L., Callot, L. and Januschowski, T. [2020], ‘Neural forecasting: Introduction and literature overview’, CoRR abs/2004.10240.
https://arxiv.org/abs/2004.10240

  * Bergstra et al. [2013] Bergstra, J., Yamins, D. and Cox, D. [2013], Making a science of model search: Hyperparameter optimization in hundreds of dimensions for vision architectures, in S. Dasgupta and D. McAllester, eds, ‘Proceedings of the 30th International Conference on Machine Learning’, Proceedings of Machine Learning Research, PMLR, Atlanta, Georgia, USA, pp. 115–123.
http://proceedings.mlr.press/v28/bergstra13.html

  * Bruzda [2020] Bruzda, J. [2020], ‘Multistep quantile forecasts for supply chain and logistics operations: bootstrapping, the GARCH model and quantile regression based approaches’, Central European Journal of Operations Research 28(1), 309–336.
https://doi.org/10.1007/s10100-018-0591-2

  * Chang et al. [2017] Chang, S., Zhang, Y., Han, W., Yu, M., Guo, X., Tan, W., Cui, X., Witbrock, M. J., Hasegawa-Johnson, M. and Huang, T. S. [2017], ‘Dilated recurrent neural networks’, CoRR abs/1710.02224.
http://arxiv.org/abs/1710.02224

  * Chung et al. [2014] Chung, J., Gülçehre, Ç., Cho, K. and Bengio, Y. [2014], ‘Empirical evaluation of gated recurrent neural networks on sequence modeling’, CoRR abs/1412.3555.
http://arxiv.org/abs/1412.3555

  * Dask Development Team [2016] Dask Development Team [2016], Dask: Library for dynamic task scheduling.
https://dask.org

  * Dauphin et al. [2016] Dauphin, Y. N., Fan, A., Auli, M. and Grangier, D. [2016], ‘Language modeling with gated convolutional networks’, CoRR abs/1612.08083.
http://arxiv.org/abs/1612.08083

  * Dudek et al. [2021] Dudek, G., Pełka, P. and Smyl, S. [2021], ‘A hybrid residual dilated lstm and exponential smoothing model for midterm electric load forecasting’, IEEE Transactions on Neural Networks and Learning Systems pp. 1–13.
  * Elamin and Fukushige [2018] Elamin, N. and Fukushige, M. [2018], ‘Quantile regression model for peak load demand forecasting with approximation by triangular distribution to avoid blackouts’, International Journal of Energy Economics and Policy 8, 119–124.
  * Elman [1990] Elman, J. L. [1990], ‘Finding structure in time’, COGNITIVE SCIENCE 14(2), 179–211.
https://crl.ucsd.edu/~elman/Papers/fsit.pdf

  * Gers et al. [1999] Gers, F. A., Cummins, F. and Schmidhuber, J. [1999], ‘Learning to forget: continual prediction with lstm’, IET Conference Proceedings pp. 850–855(5).
  * Graves [2013] Graves, A. [2013], ‘Generating sequences with recurrent neural networks’, CoRR abs/1308.0850.
http://arxiv.org/abs/1308.0850

  * Hermans and Schrauwen [2013] Hermans, M. and Schrauwen, B. [2013], Training and analysing deep recurrent neural networks, in C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani and K. Q. Weinberger, eds, ‘Advances in Neural Information Processing Systems 26’, Curran Associates, Inc., pp. 190–198.
http://papers.nips.cc/paper/5166-training-and-analysing-deep-recurrent-neural-
networks.pdf

  * Hyndman and Athanasopoulos [2018] Hyndman, R. and Athanasopoulos, G. [2018], Forecasting: Principles and Practice, 2nd edn, OTexts, Australia.
  * Hyndman et al. [2002] Hyndman, R. J., Koehler, A. B., Snyder, R. D. and Grose, S. [2002], ‘A state space framework for automatic forecasting using exponential smoothing methods’, International Journal of Forecasting 18(3), 439–454.
  * Jose and Winkler [2009] Jose, V. R. R. and Winkler, R. L. [2009], ‘Evaluating quantile assessments’, Operations Research 57(5), 1287–1297.
http://www.jstor.org/stable/25614838

  * Kalchbrenner et al. [2016] Kalchbrenner, N., Espeholt, L., Simonyan, K., van den Oord, A., Graves, A. and Kavukcuoglu, K. [2016], ‘Neural machine translation in linear time’, CoRR abs/1610.10099.
http://arxiv.org/abs/1610.10099

  * Koenker and Bassett [1978] Koenker, R. and Bassett, G. [1978], ‘Regression quantiles’, Econometrica 46(1), 33–50.
http://www.jstor.org/stable/1913643

  * Kratsios and Bilokopytov [2020] Kratsios, A. and Bilokopytov, I. [2020], ‘Non-euclidean universal approximation’, CoRR abs/2006.02341.
https://arxiv.org/abs/2006.02341

  * Lago et al. [2018] Lago, J., De Ridder, F. and De Schutter, B. [2018], ‘Forecasting spot electricity prices: Deep learning approaches and empirical comparison of traditional algorithms’, Applied Energy 221, 386 – 405.
  * Längkvist et al. [2014] Längkvist, M., Karlsson, L. and Loutfi, A. [2014], ‘A review of unsupervised feature learning and deep learning for time-series modeling’, Pattern Recognition Letters 42, 11–24.
  * Ma and Pohlman [2008] Ma, L. and Pohlman, L. [2008], ‘Return forecasts and optimal portfolio construction: A quantile regression approach’, European Journal of Finance 14, 409–425.
  * Makridakis and Hibon [2000] Makridakis, S. and Hibon, M. [2000], ‘The m3-competition: results, conclusions and implications’, International Journal of Forecasting 16(4), 451–476.  The M3- Competition.
https://www.sciencedirect.com/science/article/pii/S0169207000000571

  * Makridakis et al. [2018] Makridakis, S., Spiliotis, E. and Assimakopoulos, V. [2018], ‘The m4 competition: Results, findings, conclusion and way forward’, International Journal of Forecasting 34(4), 802 – 808.
  * Mukherjee et al. [2018] Mukherjee, S., Shankar, D., Ghosh, A., Tathawadekar, N., Kompalli, P., Sarawagi, S. and Chaudhury, K. [2018], ‘ARMDN: associative and recurrent mixture density networks for eretail demand forecasting’, CoRR abs/1803.03800.
http://arxiv.org/abs/1803.03800

  * Nascimento et al. [2019] Nascimento, R. C., Souto, Y. M., Ogasawara, E. S., Porto, F. and Bezerra, E. [2019], ‘Stconvs2s: Spatiotemporal convolutional sequence to sequence network for weather forecasting’, CoRR abs/1912.00134.
http://arxiv.org/abs/1912.00134

  * Olivares et al. [2021] Olivares, K. G., Challu, C., Marcjasz, G., Weron, R. and Dubrawski, A. [2021], ‘Neural basis expansion analysis with exogenous variables: Forecasting electricity prices with nbeatsx’, International Journal of Forecasting, submitted Working Paper version available at arXiv:2104.05522.
  * Oreshkin et al. [2019] Oreshkin, B. N., Carpov, D., Chapados, N. and Bengio, Y. [2019], ‘N-BEATS: neural basis expansion analysis for interpretable time series forecasting’, CoRR abs/1905.10437.
http://arxiv.org/abs/1905.10437

  * Oreshkin et al. [2020] Oreshkin, B. N., Carpov, D., Chapados, N. and Bengio, Y. [2020], ‘Meta-learning framework with applications to zero-shot time-series forecasting’, CoRR abs/2002.02887.
https://arxiv.org/abs/2002.02887

  * Park et al. [2020] Park, S., Yun, C., Lee, J. and Shin, J. [2020], ‘Minimum width for universal approximation’, CoRR abs/2006.08859.
https://arxiv.org/abs/2006.08859

  * Pascanu et al. [2012] Pascanu, R., Mikolov, T. and Bengio, Y. [2012], ‘Understanding the exploding gradient problem’, CoRR abs/1211.5063.
http://arxiv.org/abs/1211.5063

  * Qin et al. [2011] Qin, Y., Wang, R., Vakharia, A., Chen, Y. and Seref, M. [2011], ‘The newsvendor problem: review and directions for future research. european journal of operational research, 213, 361-374’, European Journal of Operational Research 213, 361–374.
  * Rumelhart et al. [1986] Rumelhart, D. E., Hinton, G. E. and Williams, R. J. [1986], ‘Learning Representations by Back-propagating Errors’, Nature 323(6088), 533–536.
http://www.nature.com/articles/323533a0

  * Salinas et al. [2020] Salinas, D., Flunkert, V., Gasthaus, J. and Januschowski, T. [2020], ‘Deepar: Probabilistic forecasting with autoregressive recurrent networks’, International Journal of Forecasting 36(3), 1181–1191.
  * Smyl [2019] Smyl, S. [2019], ‘A hybrid method of exponential smoothing and recurrent neural networks for time series forecasting’, International Journal of Forecasting .
  * Spyros et al. [2018] Spyros, M., Spiliotis, E. and Assimakopoulos, V. [2018], ‘Statistical and machine learning forecasting methods: Concerns and ways forward’, PLoS One .
  * Stigler [1984] Stigler, S. M. [1984], ‘Studies in the history of probability and statistics xl boscovich, simpson and a 1760 manuscript note on fitting a linear relation’, Biometrika 71(3), 615–620.
http://www.jstor.org/stable/2336572

  * Sutskever et al. [2014] Sutskever, I., Vinyals, O. and Le, Q. V. [2014], ‘Sequence to sequence learning with neural networks’, CoRR abs/1409.3215.
http://arxiv.org/abs/1409.3215

  * Van Den Oord et al. [2016] Van Den Oord, A., Dieleman, S., Zen, H., Simonyan, K., Vinyals, O., Graves, A., Kalchbrenner, N., Senior, A. W. and Kavukcuoglu, K. [2016], ‘Wavenet: A generative model for raw audio’, CoRR abs/1609.03499.
http://arxiv.org/abs/1609.03499

  * Wen et al. [2017] Wen, R., Torkkola, K., Narayanaswamy, B. and Madeka, D. [2017], ‘A Multi-Horizon Quantile Recurrent Forecaster’, arXiv e-prints p. arXiv:1711.11053.
  * Werbos [1988] Werbos, P. J. [1988], ‘Generalization of backpropagation with application to a recurrent gas market model’, Neural Networks 1(4), 339–356.
  * Werbos [1990] Werbos, P. J. [1990], ‘Backpropagation through time: what it does and how to do it’, Proceedings of the IEEE 78(10), 1550–1560.
  * Xiao [2012] Xiao, Z. [2012], 9 - time series quantile regressions, in T. Subba Rao, S. Subba Rao and C. Rao, eds, ‘Time Series Analysis: Methods and Applications’, Vol. 30 of Handbook of Statistics, Elsevier, pp. 213–257.
  * Xiao and Koenker [2009] Xiao, Z. and Koenker, R. [2009], ‘Conditional quantile estimation for generalized autoregressive conditional heteroscedasticity models’, Journal of the American Statistical Association 104(488), 1696–1712.
http://www.jstor.org/stable/40592373

  * Xu et al. [2020] Xu, D., Shi, Y., Tsang, I. W., Ong, Y.-S., Gong, C. and Shen, X. [2020], ‘Survey on multi-output learning’, IEEE Transactions on Neural Networks and Learning Systems 31(7), 2409–2429.


# Foldy-Wouthuysen transformation and multiwave states of a graphene electron
in external fields and free (2+1)-space

Alexander J. Silenko1,2,3<EMAIL_ADDRESS>1Bogoliubov Laboratory of
Theoretical Physics, Joint Institute for Nuclear Research, Dubna 141980,
Russia  2Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou
730000, China  3Research Institute for Nuclear Problems, Belarusian State
University, Minsk 220030, Belarus

###### Abstract

The relativistic Foldy-Wouthuysen transformation is used for an advanced
description of planar graphene electrons in external fields and free
(2+1)-space. It is shown that the initial Dirac equation should contain usual
Dirac matrices but not the Pauli ones. The spin of graphene electrons is not
the isotopic spin and takes the values $\displaystyle\pm 1/2$. The exact
Foldy-Wouthuysen Hamiltonian of a graphene electron in uniform and nonuniform
magnetic fields is derived. The exact energy spectrum agreeing with
experimental data and exact Foldy-Wouthuysen wave eigenfunctions are obtained.
These eigenfunctions describe multiwave (structured) states in (2+1)-space. It
is proven that the Hermite-Gauss beams exist even in the free space. In the
multiwave Hermite-Gauss states, graphene electrons acquire nonzero effective
masses dependent on a quantum number and move with group velocities which are
less than the Fermi velocity. Graphene electrons in a static electric field
also can exist in the multiwave Hermite-Gauss states defining non-spreading
coherent beams. These beams can be accelerated and decelerated.

Dirac equation in (2+1)-space; relativistic Foldy-Wouthuysen transformation;
Landau levels for graphene electrons; multiwave states in graphene

##  I Introduction

The graphene is an allotrope of carbon consisting of a single layer of atoms
arranged in a hexagonal lattice nanostructure GeimNov . It possesses many
extraordinary properties. In particular, current carriers in the graphene
(electrons) are massless. Since the graphene has a 2-dimensional structure,
graphene electrons exist in (2+1)-space. Certainly, they are described by the
Dirac equation in this space. A solution of this equation defines all
fundamental properties of the graphene electron. In the present study, we
consider the planar graphene electron.

Unlike _all_ previous studies of the graphene electron, we use the
relativistic Foldy-Wouthuysen (FW) transformation of the initial Dirac
equation developed in Refs. JMP ; PRA2015 . As a result, we fulfill the exact
description of the graphene electron in a magnetic field and the high-
precision description in an electric field. The relativistic FW transformation
can also be used for a particle in nonstationary fields PRAnonstat . The FW
representation gives one the Schrödinger form of the relativistic quantum
mechanics (QM). The solution for the magnetic field is similar to the Landau
result Landau ; LL3  and gives one the exact FW wave eigenfunctions and the
exact energy spectrum agreeing with experimental data GrapheneMiller ;
GraphenePCheng .

An analysis of paraxial multiwave states of graphene electrons in the free
(2+1)-space and a static uniform electric field demonstrates the existence of
Hermite-Gauss (HG) states and the possibility to accelerate and decelerate the
electrons without a lost of coherence. We also describe extraordinary
properties of graphene electrons in the HG states, namely, acquiring nonzero
effective masses dependent on a quantum number and moving with group
velocities which are less than the carrier velocity. The latter velocity is
equal to the Fermi one $\displaystyle v_{F}\approx 1\times 10^{6}$ m/s.

The classical description of spin in (2+1)-space fulfilled in Sec. II shows
that the spin is a real and important physical parameter. We consider
relativistic FW transformations in Sec. III. The exact solutions for a
graphene electron in a magnetic field are obtained in Sec. IV. In Secs. V and
VI, we analyze multiwave (structured) HG states and extraordinary properties
of graphene electrons in the free space and a static electric field. Section
VII summarizes the results obtained.

We assume that $\displaystyle\hbar=1,~{}c=1$ but include $\displaystyle\hbar$
and $\displaystyle c$ into some formulas when this inclusion clarifies the
problem. In the considered case, the role of $\displaystyle c$ plays the Fermi
velocity $\displaystyle v_{F}$. While the graphene electrons are massless, we
take into account the existence of massive chiral quasiparticles in the
bilayer graphene (see Ref. GrapheneGeim  and references therein) and keep the
mass in many formulas.

##  II Classical description of spin in (2+1)-space

There is a wonderful similarity between classical and quantum-mechanical
Hamiltonians and equations of spin motion in (3+1)-space. This similarity
becomes evident when the FW representation FW  is used (see Ref. PRAFW  and
references therein). It is very natural to suppose that such a similarity
remains valid in (2+1)-space. Therefore, a classical analysis of spin in
(2+1)-space plays an important role.

In this case, the antisymmetric tensor of electromagnetic field takes the form
($\displaystyle\mu,\nu=0,1,2$)

$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}=\left(\begin{array}[]{ccc}0&E_{1}&E_{2}\\\
-E_{1}&0&-B\\\ -E_{2}&B&0\end{array}\right),$ (1)

where $\displaystyle A_{\mu}=(A_{0},A_{1},A_{2})=(A_{0},-\bm{A})$ is the field
potential and $\displaystyle\bm{E}=(E_{1},E_{2})$ and $\displaystyle B$ are
the electric and magnetic fields. The difference between upper and lower
indices is only for the field potential ($\displaystyle
A^{1}=-A_{1},\,A^{2}=-A_{2}$).

Vector products of two components of any vectors correspond to third
components in the (3+1)-dimensional space:
$\displaystyle\bm{C}\times\bm{D}=C_{1}D_{2}-C_{2}D_{1}$. Certainly, this
_scalar_ quantity can have any sign. To avoid any misleading, we will often
write such a vector product in the form
$\displaystyle(\bm{C}\times\bm{D})_{3}$ demonstrating that this is a scalar.

The (2+1)-OAM tensor is given by

${\cal L}^{\mu\nu}=\sum{(x^{\mu}p^{\nu}-x^{\nu}p^{\mu})}.$ (2)

The spatial components of this tensor read

${\cal L}^{12}=-{\cal L}^{21}=L,$ (3)

where $\displaystyle L$ is the OAM which can have any sign or be equal to
zero. The summation takes place for an ensemble of particles.

In classical physics and QM, the _conventional_ spin is defined in the
particle rest frame. The total angular momentum is the sum of the orbital
angular momentum (OAM) and the rest-frame spin. While the total angular
momentum can also be defined as the sum of the OAM and the laboratory-frame
spin, this definition needs the use of noncommutative geometry (see Ref. PRAFW
and references therein).

The conventional classical (rest-frame) spin is defined as the internal OAM of
a compound particle in its rest frame. The classical spin depends on a hidden
internal motion in the compound particle. The OAM relative to any axis is
conditioned by the distance to this axis and the momentum of compound particle
as a whole. For the compound particle, there is no summation in Eq. (2).

Since the OAM can have any sign or be equal to zero, the spin also has this
property. As a result, the spin in QM of a Dirac particle can be equal to
$\displaystyle\pm 1/2$. Therefore, the spin coupling with the magnetic field,
$\displaystyle sB$, can be nonzero and can has any sign even for a fixed sign
of $\displaystyle B$. QM in the FW representation always agrees with classical
physics and this agreement has rigorous substantiation PRAFW . Classical and
quantum-mechanical equations of spin motion are similar. The spin is an
additional degree of freedom. The main distinguishing feature of spin in QM is
its quantization.

The analysis fulfilled shows that the spin in (2+1)-space is a real and
important physical parameter which cannot be reduced to the isotopic spin.

##  III Relativistic Foldy-Wouthuysen transformation in (2+1)-space

The relativistic FW transformations in (2+1)- and (3+1)-spaces are identical.
Therefore, we use in the present study the relativistic method first developed
in Ref. JMP  and presented in the final form in Ref. PRA2015 . In the general
case, this method uses an expansion of the final result in a series in
increasing powers of $\displaystyle\hbar$ and gives one _exact_ expressions
for leading terms in the FW Hamiltonian proportional to the zero and first
powers of the Planck constant and for such terms proportional to
$\displaystyle\hbar^{2}$ which describe contact interactions.

A great advantage of the FW representation FW  is the simple form of operators
corresponding to classical observables. In this representation, the
Hamiltonian and all operators are even, i.e., block-diagonal (diagonal in two
spinors). The passage to the classical limit usually reduces to a replacement
of operators in quantum-mechanical Hamiltonians and equations of motion with
the corresponding classical quantities. The possibility of such a replacement,
explicitly or implicitly used in practically all works devoted to the FW
transformation, has been rigorously proved for the stationary case in Ref.
JINRLett12 .

The initial Hamiltonian operator can be split into even and odd operators
commuting and anticommuting with the Dirac operator $\displaystyle\beta$,
respectively:

${\cal H}=\beta{\cal M}+{\cal E}+{\cal O},~{}~{}~{}\beta{\cal M}={\cal
M}\beta,~{}~{}~{}\beta{\cal E}={\cal E}\beta,~{}~{}~{}\beta{\cal O}=-{\cal
O}\beta.$ (4)

The even operators $\displaystyle{\cal M}$ and $\displaystyle{\cal E}$ and the
odd operator $\displaystyle{\cal O}$ are diagonal and off-diagonal in two
spinors, respectively. For a Dirac particle, the $\displaystyle{\cal M}$
operator usually reduces to the particle rest energy $\displaystyle mc^{2}$:

${\cal H}_{D}=\beta mc^{2}+{\cal E}+{\cal O}.$ (5)

In the general case, the relativistic FW transformation of the initial
Hamiltonian results in PRA2015 ; PRA2008

$\begin{array}[]{c}{\cal H}_{FW}=\beta\epsilon+{\cal
E}+\frac{1}{4}\Biggl{\\{}\frac{1}{2\epsilon^{2}+\\{\epsilon,{\cal
M}\\}},\biggl{(}\beta\left[{\cal O},[{\cal O},{\cal M}]\right]\\\ -[{\cal
O},[{\cal O},{\cal F}]]\biggr{)}\Biggr{\\}},\quad\epsilon=\sqrt{{\cal
M}^{2}+{\cal O}^{2}},\quad{\cal F}={\cal E}-i\hbar\frac{\partial}{\partial
t}.\end{array}$ (6)

In this equation, external fields can be nonstationary.

The (2+1) Dirac equation has the same form as the usual one but the space is
reduced:

$i\gamma^{\mu}(\partial_{\mu}+ieA_{\mu})\psi=m\psi,\qquad\mu=0,1,2.$ (7)

For graphene electron, $\displaystyle m=0$. All vectors have two components,
e.g., $\displaystyle\bm{p}=-i\hbar\left(\frac{\partial}{\partial
x^{1}},\frac{\partial}{\partial x^{2}}\right)$, and $\displaystyle p_{0}={\cal
H}$. The Dirac equation (7) takes the form

$[\gamma^{0}(p_{0}-eA_{0})-\bm{\gamma}\cdot\bm{\pi}+m]\psi=0,\qquad\bm{\pi}=\bm{p}-e\bm{A},$
(8)

where $\displaystyle\bm{\pi}$ and $\displaystyle\bm{p}$ are the kinetic and
canonical (generalized) momenta. As a result, the initial equation for the
Dirac Hamiltonian reads

${\cal H}_{D}\psi=(\bm{\alpha}\cdot\bm{\pi}+\beta
m+eA_{0})\psi,\quad\bm{\alpha}=\gamma^{0}\bm{\gamma},\quad\beta\equiv\gamma^{0}.$
(9)

For this equation,

$\displaystyle{\cal E}=eA_{0},\qquad{\cal O}=\bm{\alpha}\cdot\bm{\pi}.$

For electrons, $\displaystyle e=-|e|$.

There are different forms of the Dirac matrices in (2+1)-space. The
conventional approach used, e.g., in gravity (see Refs. Unal ; Koke ) and
electromagnetism RJackiw  consists in the reduction of the matrix space:

$\beta=\gamma^{0}=\sigma_{3},\qquad\gamma^{1}=i\sigma_{1},\qquad\gamma^{2}=i\sigma_{2},$
(10)

where $\displaystyle\sigma_{i}~{}(i=1,2,3)$ are the $\displaystyle 2\times 2$
Pauli matrices. However, the following analysis shows that this reduction
results in misleading conclusions. It is sufficient to consider the planar
graphene electron in a magnetic field. For the reduced spin matrices (10), the
relativistic FW Hamiltonian takes the form

$\begin{array}[]{c}{\cal
H}_{FW}=\sigma_{3}\sqrt{m^{2}+\bm{\pi}^{2}-e\sigma_{3}B}.\end{array}$ (11)

Since this Hamiltonian is diagonal, states with positive and negative total
energies are separated. Lower and upper components of the two-component wave
function vanish for positive-energy and negative-energy states, respectively.
As a result, Eq. (11) can be presented as follows:

$\begin{array}[]{c}{\cal
H}_{FW}=\left(\begin{array}[]{cc}\sqrt{m^{2}+\bm{\pi}^{2}-eB}&0\\\
0&-\sqrt{m^{2}+\bm{\pi}^{2}+eB}\end{array}\right).\end{array}$ (12)

Equations (11) and (12) clearly show that the Dirac equation with the reduced
matrices (10) does not describe the true spin. It follows from these equations
that, in particular, a spin-1/2 particle in positive-energy states has _only
one_ spin value. The spin characterized by Eqs. (11) and (12) is not an
additional degree of freedom. These properties sharply contrast with the
results obtained in the precedent section. Therefore, we can conclude that the
conventional approach consisting in the reduction of the matrix space is
mistaken.

The reduction of the matrix space is not necessary and one can use the
original $\displaystyle 4\times 4$ Dirac matrices Lozovik . We apply this
approach in the present study. When the original Dirac matrices are utilized,
the relativistic FW Hamiltonian has the form

$\begin{array}[]{c}{\cal
H}_{FW}=\beta\epsilon+eA_{0}+\frac{e}{8}\biggl{\\{}\frac{1}{\epsilon(\epsilon+m)},\Sigma_{3}(\bm{\pi}\times\bm{E}\\\
-\bm{E}\times\bm{\pi})_{3}\biggr{\\}},\qquad\epsilon=\sqrt{m^{2}+\bm{\pi}^{2}-e\Sigma_{3}B}.\end{array}$
(13)

In this equation, only the field potentials and their first derivatives are
taken into account. The applied method gives _exact_ expressions for such
terms PRA2015 . Here and below, we use the standard form of the Dirac matrices
BLP . The total energy can be positive and negative.

Equation (13) explicitly shows that the rudimentary spin in (2+1)-space is an
important characteristic of the particle polarization but not an isotopic
spin. The eigenvalues $\displaystyle\mathfrak{S}_{3}=\pm 1$ of the operator
$\displaystyle\Sigma_{3}$ are scalar values but not projections of any vector.
The spin is an additional degree of freedom and its value can be positive or
negative. Thus, the use of the original $\displaystyle 4\times 4$ Dirac
matrices conserves the agreement between QM in the FW representation and
classical physics. This fundamental conclusion can be extended on Dirac
particles in (2+1)-gravity.

We should also mention that spin interactions of anyons with an
electromagnetic field have been considered in Ref. Horvathy .

##  IV Exact solutions for a graphene electron in a magnetic field

When the original $\displaystyle 4\times 4$ Dirac matrices are used, the FW
Hamiltonian for a graphene electron in a static magnetic field is exact and
reads

$\begin{array}[]{c}{\cal
H}_{FW}=\beta\sqrt{m^{2}+\bm{\pi}^{2}-\Sigma_{3}B}.\end{array}$ (14)

We underline that the magnetic field can be nonuniform. A similar situation
takes place for a Dirac particle JMP ; Case ; Energy1 ; Energy3  and a spin-1
particle with the normal ($\displaystyle g=2$) magnetic moment spin1  in
(3+1)-space. The corresponding FW Hamiltonians are also exact. The similarity
of the FW Hamiltonians in (2+1)- and (3+1)-spaces substantially simplifies
finding the exact energy spectrum and exact wave eigenfunctions for the
Hamiltonian (14).

Fortunately, exact solutions for a graphene electron ($\displaystyle m=0$) in
a uniform magnetic field can be obtained from the Landau solutions for
(3+1)-space. We suppose that $\displaystyle B>0$. There are two appropriate
gauges LL3 . For the symmetric gauge $\displaystyle A_{\phi}=Br/2,\,A_{r}=0$,
the squared Eq. (14) takes the form

$\begin{array}[]{c}-\nabla^{2}+ieB\frac{\partial}{\partial\phi}+\frac{e^{2}B^{2}r^{2}}{4}-2es_{z}B=\mathfrak{E}^{2},\\\
\nabla^{2}=\frac{\partial^{2}}{\partial
r^{2}}+\frac{1}{r}\frac{\partial}{\partial
r}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\phi^{2}},\end{array}$ (15)

where $\displaystyle\mathfrak{E}$ is the energy.

The energy spectrum is defined by the well-known eigenvalues of the operator
$\displaystyle\bm{\pi}^{2}$ LL3  and is gived by ($\displaystyle e=-|e|$)

$\begin{array}[]{c}\mathfrak{E}=\pm\sqrt{\left(2n+1+|\ell|+\ell+2s_{z}\right)|e|B}.\end{array}$
(16)

Here $\displaystyle s_{z}=\pm 1/2$, $\displaystyle n=0,1,2,\dots$, and
$\displaystyle\ell$ (more precisely, $\displaystyle\hbar\ell$) is the OAM
which is a scalar (but not a pseudovector) in (2+1)-space. $\displaystyle\ell$
is integer. The positive and negative signs define states with positive and
negative total energies (graphene electrons and holes, respectively). The
corresponding wave eigenfunctions are also defined by the Landau solution for
the symmetric gauge LL3  (see also Refs. Energy3 ; paraxialLandau ). They are
based on the Laguerre polynomials and are given by

$\begin{array}[]{c}\Phi_{FW}={\cal
A}\exp{(i\ell\phi)},\qquad\int{\Phi_{FW}^{\dagger}\Phi_{FW}rdrd\phi}=1,\\\
{\cal
A}=\frac{C_{n\ell}}{w_{m}}\left(\frac{\sqrt{2}r}{w_{m}}\right)^{|\ell|}L_{n}^{|\ell|}\left(\frac{2r^{2}}{w^{2}_{m}}\right)\exp{\left(-\frac{r^{2}}{w^{2}_{m}}\right)}\eta,\\\
C_{n\ell}=\sqrt{\frac{2n!}{\pi(n+|\ell|)!}},\qquad
w_{m}=\frac{2}{\sqrt{|e|B}},\end{array}$ (17)

where the real function $\displaystyle{\cal A}$ defines the amplitude of the
beam, and $\displaystyle L_{n}^{|\ell|}$ is the generalized Laguerre
polynomial. The wave function $\displaystyle\Phi_{FW}$ is the upper spinor for
the FW wave function in positive-energy states:
$\displaystyle\Psi_{FW}=\left(\begin{array}[]{c}\Phi_{FW}\\\
0\end{array}\right)$. If negative-energy states are disregarded, the spin
function $\displaystyle\eta$ is an eigenfunction of the Dirac operator
$\displaystyle\Sigma_{z}$ (cf. Ref. Energy3 ):

$\displaystyle\Sigma_{z}\eta^{\pm}=\pm\eta^{\pm},\quad\eta^{+}=\left(\begin{array}[]{c}1\\\
0\end{array}\right),\quad\eta^{-}=\left(\begin{array}[]{c}0\\\
1\end{array}\right).$

Here 0 and 1 are the zero and unit $\displaystyle 2\times 2$ matrices.

Another gauge is also possible Landau ; LL3 : $\displaystyle
A_{x}=-By,\,A_{y}=0$. For this gauge, the energy eigenvalues can be obtained
from the following equation:

$\begin{array}[]{c}-\nabla^{2}+2ieBy\frac{\partial}{\partial
x}+e^{2}B^{2}y^{2}-2es_{z}B=\mathfrak{E}^{2},\\\
\nabla^{2}=\frac{\partial^{2}}{\partial x^{2}}+\frac{\partial^{2}}{\partial
y^{2}}.\end{array}$ (18)

The energy spectrum is given by

$\begin{array}[]{c}\mathfrak{E}=\pm\sqrt{\left(2N+1+2s_{z}\right)|e|B},\end{array}$
(19)

where $\displaystyle N=0,1,2,\dots$. Evidently, Eqs. (16) and (19) define the
same set of energy levels. We underline the importance of the spin. The
eigenfunctions can be simply taken from Refs. Landau ; LL3  (with the
rejection of $\displaystyle z$ coordinate and the addition of the spin wave
function):

$\begin{array}[]{c}\Phi_{FW}=\exp{(ip_{x}x)}\psi(y)\eta,\qquad\int{\psi^{\dagger}\psi
dy}=1,\\\
\psi(y)=C_{N}\exp{\left[-\frac{2(y-y_{0})^{2}}{w_{m}^{2}}\right]}H_{N}\left[\frac{2(y-y_{0})}{w_{m}}\right],\\\
C_{N}=\frac{1}{\pi^{1/4}w_{m}^{1/2}\sqrt{2^{N-1}N!}},\qquad
y_{0}=-\frac{p_{x}}{eB},\end{array}$ (20)

where $\displaystyle H_{N}$ is the Hermite polynomial. We note that the exact
solutions (17) and (20) define multiwave states formed by infinite continuums
of partial de Broglie waves but not by single de Broglie waves having definite
momenta.

Equations (16) and (19) perfectly agree with experimental data GrapheneMiller
; GraphenePCheng  and do not contradict to earlier theoretical derivations
RJackiw ; Haldane . However, we take into account the spin contribution. In
particular, the existence of the zero energy level
$\displaystyle\mathfrak{E}=0$ can be explained only with allowance for this
contribution. In fact, all energy levels are shifted by the value of the Dirac
energy $\displaystyle E_{D}$ GraphenePCheng .

##  V Multiwave Hermite-Gauss states and extraordinary properties of graphene
electrons in the free space

Fortunately, multiwave HG states of graphene electrons in the free (2+1)-space
can also be deduced from the corresponding HG states of particles in the free
(3+1)-space. The latter states are well known. They are defined by the wave
functions localized in $\displaystyle x$ and $\displaystyle y$ dimensions (see
Refs. Kogelnik ; Siegman ; Pampaloni ):
$\displaystyle\Phi_{FW}=\psi_{n}(x,z)\psi_{m}(y,z)$. These functions are
derived in the paraxial approximation
($\displaystyle|p_{z}|>>|\bm{p}_{\bot}|$) and are the solutions of the
paraxial equation. In the relativistic QM in the FW representation, this
equation has the form photonPRA

$\begin{array}[]{c}\left(\nabla^{2}_{\bot}+2ik\frac{\partial}{\partial
z}\right)\\!\Phi_{FW}=0,\quad\nabla^{2}_{\bot}=\frac{\partial^{2}}{\partial
x^{2}}+\frac{\partial^{2}}{\partial y^{2}},\end{array}$ (21)

where $\displaystyle\hbar k=p=\sqrt{\bm{p}_{\bot}^{2}+p_{z}^{2}}$ and
$\displaystyle\hbar$ is omitted. In Eq. (21), an appropriate shift of the
squared particle momentum has been made paraxialLandau ; arXiv . It is similar
to a shift of the zero energy level in Schrödinger QM.

Wave functions in optics are similar to FW ones. The functions
$\displaystyle\psi_{n}(x,z)$ and $\displaystyle\psi_{m}(y,z)$ are normalized
to unit:

$\displaystyle\int{\psi_{n}^{\dagger}(x,z)\psi_{n}(x,z)dx}=1,\quad\int{\psi_{m}^{\dagger}(y,z)\psi_{m}(y,z)dy}=1.$

However, a HG wave in the free (3+1)-space can be localized in $\displaystyle
x$ direction and unlocalized in $\displaystyle y$ one. Such a beam can also
exist in the free (2+1)-space. In this space, the paraxial equation takes the
form

$\begin{array}[]{c}\left(\frac{\partial^{2}}{\partial
x^{2}}+2ik\frac{\partial}{\partial y}\right)\Phi_{FW}=0.\end{array}$ (22)

The corresponding wave function reads

$\begin{array}[]{c}\Phi_{FW}=\exp{(ip_{y}y)}\psi_{N}(x,y)\eta,\\\
\int{\psi_{N}^{\dagger}(x,y)\psi_{N}(x,y)dx}=1,\\\
\psi_{N}(x,y)=C_{N}H_{N}\left[\frac{\sqrt{2}x}{w(y)}\right]\exp{\biggl{[}-\frac{x^{2}}{w^{2}(y)}-ik\frac{x^{2}}{2R(y)}}\\\
{+i\left(N+\frac{1}{2}\right)\Phi_{G}(y)\biggr{]}},\quad
C_{N}=\frac{2^{1/4}}{\pi^{1/4}w^{1/2}(y)\sqrt{2^{N}N!}},\\\
w(y)=w_{0}\sqrt{1+\frac{4y^{2}}{k^{2}w_{0}^{4}}},\quad
R(y)=y+\frac{k^{2}w_{0}^{4}}{4y},\\\
\Phi_{G}(y)=\arctan{\left[\frac{kw^{2}(y)}{2R(y)}\right]}=\arctan{\left(\frac{2y}{kw_{0}^{2}}\right)},\end{array}$
(23)

where $\displaystyle k$ is the wave number of the beam, $\displaystyle w_{0}$
is the minimum beam width, $\displaystyle R(y)$ is the radius of curvature of
the wave front, and $\displaystyle\Phi_{G}(y)$ is the Gouy phase.

Equation (23) describes multiwave (structured) HG beams being infinite
continuums of partial de Broglie waves in the free (2+1)-space. These beams do
not spread. As well as the corresponding beams in the free (3+1)-space
photonPRA , such beams possess extraordinary properties. First of all, current
carriers (electrons and holes) in multiwave states move slower than the same
carriers in structureless ones. QM gives the simple explanation of the similar
effect in the free (3+1)-space photonPRA .

All beam parameters are defined by expectation values or eigenvalues of
related operators. In the considered case, the group velocity operator depends
on a hidden motion along the $\displaystyle x$ axis. As follows from Eq. (14)
at $\displaystyle B=0$,

$\begin{array}[]{c}\bm{v}=\frac{\partial{\cal
H}_{FW}}{\partial\bm{p}}=\frac{v_{F}\bm{p}}{p},\quad v=v_{F}.\end{array}$ (24)

We use the term “hidden motion” for a motion which does not contribute to
expectation values of operators defining some components of the velocity and
momentum but affects expectation values of squares of both these operators and
eigenvalues of the energy operator. In the considered case,
$\displaystyle\langle v_{x}\rangle=0,~{}\langle v_{y}\rangle\neq 0$. However,
the absolute value of the total momentum, $\displaystyle p$, is equal to
$\displaystyle\mathfrak{E}/v_{F}$ and is the same for all partial waves. The
paraxial approximation leads to $\displaystyle
p_{y}=\sqrt{p^{2}-p_{x}^{2}}\approx p\left(1-\frac{p_{x}^{2}}{2p^{2}}\right)$.
Thus, the beam moves as a whole with the averaged group velocity (cf. Ref.
photonPRA )

$\begin{array}[]{c}\langle\bm{v}\rangle=v^{2}_{F}\frac{\langle
p_{y}\rangle\bm{e}_{y}}{\mathfrak{E}}\approx
v_{F}\bm{e}_{y}\left(1-\frac{v^{2}_{F}\langle
p_{x}^{2}\rangle}{2\mathfrak{E}^{2}}\right).\end{array}$ (25)

This effect is similar to the subluminality of structured light (see Ref.
photonPRA  and references therein).

Graphene electrons in the HG states also acquire nonzero effective masses
dependent on a quantum number. In the laboratory frame, a graphene electron
can be considered as a centroid with the constant laboratory frame energy
$\displaystyle\mathfrak{E}$ and the velocity
$\displaystyle|\langle\bm{v}\rangle|<v_{F}$ defined by Eq. (25). Certainly,
the centroid has also the nonzero effective mass (cf. Ref. photonPRA )

$M=\frac{\sqrt{\mathfrak{E}^{2}-v_{F}^{2}\langle
p_{y}^{2}\rangle}}{v_{F}^{2}}=\frac{\langle p_{x}^{2}\rangle^{1/2}}{v_{F}}.$
(26)

##  VI Multiwave states of graphene electrons in a static electric field

The following analysis shows that graphene electrons in a static electric
field can be in multiwave states. We can suppose that the static and
nonuniform electric field is negligible far from the field source and a
graphene electron is in a multiwave HG state in this area. Such a graphene
beam is accelerated or decelerated by the electric field. Far from the field
source, the potential energy vanishes. Importantly, the _total_ energy of each
partial beam remains unchanged. As a result, the beam remains coherent and
non-spreading. If the potential energy in any fixed point with the radius
vector $\displaystyle\bm{R}_{0}$ is equal to $\displaystyle
U_{0}=eA_{0}(\bm{R}_{0})$, it becomes equal to $\displaystyle
U=eA_{0}(\bm{R})=e[A_{0}(\bm{R}_{0})-\bm{E}(\bm{R}_{0})\cdot\bm{r}]$
($\displaystyle\bm{r}=\bm{R}-\bm{R}_{0}$) in a nearby point with the radius
vector $\displaystyle\bm{R}$. The simplest situation takes place when the
electric field is collinear to the beam direction ($\displaystyle y$ axis). In
this case, $\displaystyle p_{x}$ is not affected by the electric field. The
paraxial approximation results in

$\begin{array}[]{c}\epsilon=p=\sqrt{p_{x}^{2}+p_{y}^{2}}\approx
p_{y}+\frac{p_{x}^{2}}{2p}.\end{array}$ (27)

For graphene electrons, the Fermi velocity is equivalent to $\displaystyle c$
for photons and is here omitted. We suppose that $\displaystyle p_{y}>0$ and
$\displaystyle\bm{r}=(x,y)$. Since $\displaystyle p_{x}$ is small, the last
term in Eq. (13) can be neglected. As a result, the energy in positive-energy
states is given by

$\begin{array}[]{c}\mathfrak{E}=p+U_{0}-E_{y}y.\end{array}$ (28)

When $\displaystyle y=0$, $\displaystyle p=p_{0}=\mathfrak{E}-U_{0}$. In the
general case, $\displaystyle p=p_{0}+E_{y}y$. Therefore, Eq. (27) takes the
form

$\begin{array}[]{c}2(p_{0}+E_{y}y)^{2}\approx
2(p_{0}+E_{y}y)p_{y}+p_{x}^{2}.\end{array}$ (29)

Since $\displaystyle E_{y}=const$, the operator expressions $\displaystyle
yp_{y}$ and $\displaystyle p_{y}y$ [$\displaystyle p_{y}=-i\partial/(\partial
y)$] are equivalent for the next transformation. This equation can be
presented in the form

$\begin{array}[]{c}p_{x}^{2}+2(p_{0}+E_{y}y)p_{y}-4p_{0}E_{y}y-2E_{y}^{2}y^{2}\approx
2p_{0}^{2}.\end{array}$ (30)

When we denote $\displaystyle k_{0}=p_{0}/\hbar$ and omit
$\displaystyle\hbar$, the operator form of Eq. (30) reads

$\begin{array}[]{c}\left[-\frac{\partial^{2}}{\partial
x^{2}}-2i(k_{0}+E_{y}y)\frac{\partial}{\partial
y}-2(2k_{0}+E_{y}y)E_{y}y\right]\Phi_{FW}\\\ \approx
2k_{0}^{2}\Phi_{FW}.\end{array}$ (31)

The substitution $\displaystyle\Phi_{FW}=\exp{(ik_{0}y)}\Psi$ brings the
corresponding paraxial equation

$\begin{array}[]{c}\left[\frac{\partial^{2}}{\partial
x^{2}}+2i(k_{0}+E_{y}y)\frac{\partial}{\partial
y}+2(2k_{0}+E_{y}y)E_{y}y\right]\Psi=0.\end{array}$ (32)

A comparative analysis of Eqs. (22) and (32) shows that wave functions of
graphene electrons in the free (2+1)-space and a static electric field are
similar and have the form (23) in both cases. However, in the latter case the
parameters $\displaystyle w(y),\,R(y)$, and $\displaystyle\Phi_{G}(y)$ are not
defined by Eq. (23) and differently depend on $\displaystyle y$.

Thus, graphene electrons in a static electric field can be in multiwave states
which characterize non-spreading coherent beams. As a result, a static
electric field can accelerate and decelerate HG beams of graphene electrons.

##  VII Summary

The relativistic FW transformation has been used for an advanced description
of free and interacting planar graphene electrons. It has been proven that the
initial Dirac equation in (2+1)-space should be based on $\displaystyle
4\times 4$ Dirac matrices but not on $\displaystyle 2\times 2$ Pauli matrices.
New properties of graphene electrons in free (2+1)-space and external fields
have been found. The exact FW Hamiltonian of a graphene electron in uniform
and nonuniform magnetic field has been derived. It has been shown that spin
effects are rather important and the spin takes the values $\displaystyle\pm
1/2$. The exact energy spectrum agreeing with experimental data GrapheneMiller
; GraphenePCheng  and exact FW wave eigenfunctions have been obtained. These
eigenfunctions describe multiwave states in (2+1)-space. It has been proven
that HG beams exist even in the free space. Two extraordinary properties of
graphene electrons in the HG states are acquiring nonzero effective masses
dependent on a quantum number and moving with group velocities which are less
than the carrier velocity. In a static electric field, graphene electrons also
can be in multiwave HG states defining non-spreading coherent beams. This
means that HG beams of graphene electrons can be accelerated and decelerated.

## Acknowledgments

The author acknowledges the support by the Chinese Academy of Sciences
President’s International Fellowship Initiative (Grant No. 2019VMA0019) and
hospitality and support by the Institute of Modern Physics of the Chinese
Academy of Sciences.

## References

  * (1) A. K. Geim, K. S. Novoselov, The rise of graphene, Nature Materials 6, 183 (2007).
  * (2) A. J. Silenko, Foldy-Wouthuysen transformation for relativistic particles in external fields, J. Math. Phys. 44, 2952 (2003).
  * (3) A. J. Silenko, General method of the relativistic Foldy-Wouthuysen transformation and proof of validity of the Foldy-Wouthuysen Hamiltonian, Phys. Rev. A 91, 022103 (2015).
  * (4) A. J. Silenko, Energy expectation values of a particle in nonstationary fields, Phys. Rev. A 91, 012111 (2015).
  * (5) L. D. Landau, Diamagnetismus der Metalle, Z. Phys. 64, 629 (1930).
  * (6) L. D. Landau, E. M. Lifshitz, _Quantum Mechanics. Non-Relativistic Theory_ , 3rd ed. (Pergamon Press, Oxford, 1977), pp. 458-461.
  * (7) D. L. Miller, K. D. Kubista, G. M. Rutter, M. Ruan, W. A. de Heer, P. N. First, J. A. Stroscio, Observing the Quantization of Zero Mass Carriers Observing the Quantization of Zero Mass Carriers in Graphene, Science 324, 924 (2009).
  * (8) P. Cheng _et al._ , Landau Quantization of Topological Surface States in Bi2Se3, Phys. Rev. Lett. 105, 076801 (2010).
  * (9) A. K. Geim, Graphene: Status and Prospects, Science 324, 1530 (2009).
  * (10) L. L. Foldy, S. A. Wouthuysen, On the Dirac Theory of Spin 1/2 Particles and Its Non-Relativistic Limit, Phys. Rev. 78, 29 (1950).
  * (11) Liping Zou, Pengming Zhang, and A. J. Silenko, Position and spin in relativistic quantum mechanics, Phys. Rev. A 101, 032117 (2020).
  * (12) A. J. Silenko, Classical limit of relativistic quantum mechanical equations in the Foldy-Wouthuysen representation, Pis’ma Zh. Fiz. Elem. Chast. Atom. Yadra 10, 144 (2013) [Phys. Part. Nucl. Lett. 10, 91 (2013)].
  * (13) A. J. Silenko, Foldy-Wouthuysen transformation and semiclassical limit for relativistic particles in strong external fields, Phys. Rev. A 77, 012116 (2008).
  * (14) Y. Sucu, N. Ünal, Exact solution of Dirac equation in 2+1 dimensional gravity, Found. Phys. 48, 052503 (2007).
  * (15) C. Koke, C. Noh, and D. G. Angelakis, Dirac equation on a square waveguide lattice with site-dependent coupling strengths and the gravitational Aharonov-Bohm effect, Phys. Rev. D 102, 013514 (2020).
  * (16) R. Jackiw, Fractional charge and zero modes for planar systems in a magnetic field, Phys. Rev. D 29, 2375 (1984).
  * (17) Y. E. Lozovik, S. P. Merkulova and A. A. Sokolik, Collective electron phenomena in graphene, Phys.-Usp. 51, 727 (2008).
  * (18) V. B. Berestetskii, E. M. Lifshitz, and L. P. Pitayevskii, Quantum Electrodynamics, 2nd ed. (Pergamon, Oxford, 1982).
  * (19) P. A. Horváthy, L. Martina, P. Stichel, Comments on spin-orbit interaction of anyons, Mod. Phys. Lett. A 20, 1177 (2005).
  * (20) K. M. Case, Some Generalizations of the Foldy-Wouthuysen Transformation, Phys. Rev. 95, 1323 (1954).
  * (21) W. Tsai, Energy eigenvalues for charged particles in a homogeneous magnetic field – an application of the Foldy-Wouthuysen transformation, Phys. Rev. D 7, 1945 (1973).
  * (22) A. J. Silenko, Connection between wave functions in the Dirac and Foldy-Wouthuysen representations, Phys. Part. Nucl. Lett. 5, 501 (2008).
  * (23) A. J. Silenko, High precision description and new properties of a spin-1 particle in a magnetic field, Phys. Rev. D 89, 121701(R) (2014).
  * (24) Liping Zou, Pengming Zhang, and A. J. Silenko, Paraxial wave function and Gouy phase for a relativistic electron in a uniform magnetic field, J. Phys. G: Nucl. Part. Phys. 47, 055003 (2020).
  * (25) Liping Zou, Pengming Zhang, and A. J. Silenko, General quantum-mechanical solution for twisted electrons in a uniform magnetic field, Phys. Rev. A 103, L010201 (2021).
  * (26) F. D. M. Haldane, Model for a Quantum Hall Effect without Landau Levels: Condensed-Matter Realization of the “Parity Anomaly”, Phys. Rev. Lett. 61, 2015 (1988).
  * (27) H. Kogelnik and T. Li, Laser Beams and Resonators, Appl. Opt. 5, 1550 (1966).
  * (28) A. E. Siegman, _Lasers_ (University Science Books, Sausalito, 1986).
  * (29) F. Pampaloni, J. Enderlein, Gaussian, Hermite-Gaussian, and Laguerre-Gaussian beams: A primer, arXiv:physics/0410021 (2004).
  * (30) A. J. Silenko, Pengming Zhang, and Liping Zou, Relativistic quantum-mechanical description of twisted paraxial electron and photon beams, Phys. Rev. A 100, 030101(R) (2019).