# Modeling and Forecasting COVID-19 Cases using Latent Subpopulations

Roberto Vega
Department of Computing Science
University of Alberta
Alberta Machine Intelligence Institute
<EMAIL_ADDRESS>
&Zehra Shah
Department of Computing Science
University of Alberta
Alberta Machine Intelligence Institute
&Pouria Ramazi
Department of Mathematics & Statistics
Brock University
&Russell Greiner
Department of Computing Science
University of Alberta
Alberta Machine Intelligence Institute
<EMAIL_ADDRESS>

###### Abstract

Classical epidemiological models assume homogeneous populations. There have
been important extensions to model heterogeneous populations, when the
identity of the sub-populations is known, such as age group or geographical
location. Here, we propose two new methods to model the number of people
infected with COVID-19 over time, each as a linear combination of latent sub-
populations – i.e., when we do not know which person is in which sub-
population, and the only available observations are the aggregates across all
sub-populations. Method#1 is a dictionary-based approach, which begins with a
large number of pre-defined sub-population models (each with its own starting
time, shape, etc), then determines the (positive) weight of small (learned)
number of sub-populations. Method#2 is a mixture-of-$M$ fittable curves, where
$M$, the number of sub-populations to use, is given by the user. Both methods
are compatible with any parametric model; here we demonstrate their use with
first (a) Gaussian curves and then (b) SIR trajectories. We empirically show
the performance of the proposed methods, first in (i) modeling the observed
data and then in (ii) forecasting the number of infected people 1- to 4-weeks
in advance. Across 187 countries, we show that the dictionary approach had the
lowest mean absolute percentage error and also the lowest variance when
compared with classical SIR models and moreover, it was a strong baseline that
outperforms many of the models developed for COVID-19 forecasting.

##  1 Introduction

Classical epidemiological models, such as the Susceptible-Infected-Removed
model (SIR), have been extensively used to study how diseases spread. These
models work under the assumption of a constant, homogeneous population, and so
their performance degrades when analyzing heterogeneous populations – e.g.,
when the population consists of communities that have different behaviours, or
that are geographically separated (Calvetti et al., 2020). In particular, the
heterogeneity of disease transmission has a significant effect on
compartmental model dynamics, and ignoring this heterogeneity has has been
shown to over-estimate the height of the epidemic peak and the number of
infected individuals at the end of the epidemic (Dolbeault and Turinici, 2020,
2021). Similarly, Ellison (2020) describes how using homogeneous models to
describe an infectious disease in a heterogeneous population might lead to
wrong conclusions in terms of disease forecasting, herd immunity, and the
impact of government policies to control the spread of the disease. Hence it
is crucial to develop models that consider population heterogeneity. This
paper proposes two such methods that model the observed number of people
infected by COVID-19 as a linear combination of independent, unknown, sub-
populations. It also empirically demonstrates that these methods can
effectively model, and forecast, the number of infected people, across 187
countries.

We first propose a dictionary-based approach, that uses a predefined set of
plausible curves, that each describes the number of infected people over time;
see Figure 1. It then finds a linear combination of a small (learned) subset
of these curves that best explain the observed data. Each of the included
curves (those with a non-zero coefficient in this linear combination) can be
interpreted as an independent sub-population. We explored two types of prior
curves: either each is a (a) normalized Gaussian probability density
functions, or is an (b) SIR curve with its own learned parameters and initial
states.

Figure 1: The blue signal on the right is reconstructed as a linear
combination of the orange signals in the dictionary on the right.

The second method is mixture of (a fixed number, $M$, of) fittable curves,
where $M$, the number of sub-populations, is determined by the user. We tested
a mixture of $M=3$ Gaussian curves, and a mixture of $M=3$ SIR models. For the
mixture of SIR models, we modified the classical SIR model to allow for shifts
of the starting time-point of each SIR trajectory via a Kronecker delta
function. We evaluate these models based on both fit to the data, and also the
accuracy (mean absolute percentage error: MAPE), for predicting the number of
infected people 1– to 4-weeks in advance. Note this approach may help us to
capture the general trend of how the number of infected people change over
time, and so might be very useful when the objective is to describe the main
behaviour of the data.

To motivate our approach, consider the three simulated independent “SIR sub-
populations" shown by the dashed lines in Figure 2(a), whose orange line
represents the sub-population in a neighborhood A, while the other two sub-
populations are both located in a neighborhood B, but have different
behaviours – green being more socially active than red. But as the number of
infected people are reported only at the state level (not the neighborhood
level), we only observe the sum of the three sub-populations – the solid blue
line. This is the only information used by the algorithms modeling the number
of infected people – note they do not know have access to the number of sub-
populations in the data, nor the number of infections within each.

Figure 2: Simple example showing the rationale behind our proposed approach.
(a) Simulated data consisting of three sub-populations, and the observed data.
(b) Best fit of a classical (single) SIR model. (c) Curve reconstructed by
using Gaussian dictionaries; using (d) the coefficients associated with each
curve in the Gaussian dictionary. (e) Curve reconstructed by the mixture of
$M$ = 2 Gaussians. (f) Curve reconstructed by the mixture of $M$ = 5
Gaussians.

Figure 2(b) shows the curve corresponding to the best fit of a classical SIR
model. Note that it does a poor job describing the data, because the SIR model
(and most of its extensions) are unimodal, meaning they cannot accurately fit
multi-modal data.

Our proposed approach models the observed data as a combination of several
sub-populations. Since the number of sub-populations is unknown, the
dictionary approach uses as many sub-populations as needed to create a curve
that closely matches the observations (in terms of mean squared error). Figure
2(c) shows how a Gaussian dictionary appropriately models the data, albeit
using 114 Gaussian curves, each with the coefficient shown in Figure 2(d). The
mismatch between the number of real sub-populations and the number of sub-
populations estimated by the model arises in part because the observed data is
in fact not a mixture of Gaussians. In general, dictionary approaches can
accurately model the observed data, provided they have a sufficiently large
and diverse set of possible curves. However, this may require a large number
of curves, which can make the resulting combination difficult to interpret and
also can mean the forecasting version will overfit.

Figures 2(e) and (f) show the curves reconstructed by the mixture of $M=2$ and
$M=5$ fittable Gaussians. While the mixtures do capture the main trend of the
data, they clearly have a higher error than the dictionary approach. Of
course, as the number of curves, $M$, increases, the error decreases. Despite
showing a higher error, as this approach lets the user control the number of
curves to be used for modeling the data, it potentially allows an easier
analysis over the identity of the curves. For example, it would be easier to
interpret Figure 2(e) to mean that one curve belongs to a neighborhood A and
the other to neighborhood B, than trying to do a similar analysis for the 114
curves selected from the Gaussian dictionary (Figure 2(d) ).

###  1.1 Related work

Most of the efforts in modeling the evolution of epidemics with heterogeneous
populations fall under the umbrella term _meta-population models_ (Calvetti et
al., 2020; Tan et al., 2021), which divide the individuals into a specified
number of sub-populations, each representing either a spatial area (a country,
a city, a neighbourhood, etc.), or a known variable of interest (age, sex,
etc). They then use a compartmental model, such as an SIR or one of its
extensions, for each sub-population, and perhaps also include some additional
terms that describe the dynamics among the sub-populations (Calvetti et al.,
2020; Contreras et al., 2020; Gerasimov et al., 2021; Duan and Jin, 2022). The
sub-populations are usually allowed to have different interaction rates (Albi
et al., 2021). These models implicitly assume that the infections start at the
same time across the different sub-populations. The quality of the overall
model is tied to the quality of the selected sub-populations, and to the
accuracy in the modeling of the interaction among them (Kong et al., 2022).
Our work extends those results by allowing latent sub-populations, which can
start at different times and can be of any parametric family.

There are some studies that use latent variables to model the uncertainty
related to the reporting of newly infected people over time (Bartolucci and
Farcomeni, 2022; Shi et al., 2022). Here, the latent variables represent the
real number of new infections, while the reported quantities are the
observations associated with such latent variables, perhaps modeled using
Hidden Markov Models. Some studies additionally incorporate the possibility of
having observed sub-populations at the geographical level to account for
interactions between populations (Bartolucci and Farcomeni, 2022). Similarly,
Tan et al. (2021) described a state-space model where latent variables
represent the compartments of an SEIRD (Susceptible, Exposed, Infected,
Recovered, Deceased) model, while the observations are only the confirmed new
infections and cumulative number of deaths. They used a meta-population
framework that models the spread of the disease across separate geographical
sub-populations (e.g., two different cities). Each sub-population might have
different transition parameters, but their estimation requires the individual
observations for each sub-population – i.e., the daily number of infections in
city 1 and in city 2. Of course, this is a different type of latency, as this
is still assuming a single homogeneous population, but where the values for
some compartments are not observed. By contrast, we propose the use of latent
variables that represent the different sub-populations. Therefore, the
observed data is a single time-series that is the _sum_ of the time-series of
the individual sub-populations.

A different approach for modeling heterogeneous populations are the agent-
based models, which model a sample of the individuals in the population.
Nepomuceno et al. (2016) describe the general theoretical framework of an
agent-based modeling approach for infectious disease dynamics, where each
agent is an individual (potentially with different features based on their
age, location, activities, etc.), and the aggregation of individual behaviors
is an estimate of the population dynamics as a whole. Meta-population models
can be adjusted to approximate the dynamics of complex agent-based models
(Zachreson et al., 2022). However, the usability of agent-based models is
limited due to their computational complexity.

Very relevant to our approach is the work of Sarker et al. (2022), who also
consider latent sub-populations, where the learner can only observe a single
time-series that is the aggregate of the different sub-populations. They model
the number of infected people over time as a mixture of Gaussians, and propose
using a peak-finding heuristic to determine the number of sub-populations
(each modeled as a Gaussian) to use in their mixture. Our work extends theirs
as we additionally show how to model an epidemic using a mixture of SIRs, and
propose a general dictionary approach that allows us to model the number of
infected people as a combination of any parametric model.

##  2 Materials and Methods

We compare three different methods to model the dynamics of the number of
people infected with COVID-19: (1) The traditional SIR model with learned
parameters. (2) A dictionary approach that contains a set of predefined
‘candidate’ curves that can be combined to describe the COVID-19 dynamics
(where each curve describes a sub-population) and an optimization algorithm
determines which sub-populations to include (and with what weight). (3) A
mixture approach that combines a predetermined number of parameterized curves
(e.g., 2 or 3), which uses an optimization algorithm to determine the
parameters for each curve.

Importantly, methods 2 and 3 assume that the sub-populations are latent –
i.e., we only observe the total number of infected people in the entire
population, and not the number of infected people in each sub-population.
Also, the ‘identity’ of the sub-populations is a hidden variable – i.e., we do
not know which people are included in each sub-population. This allow us to
model sub-populations that are not even based on observed variables, such as
sex, age and geographical location (Kong et al., 2022). This means each sub-
population could perhaps represent people with some unobserved behaviour
(e.g., about vaccination, or masking) or some other set of characteristics.

###  2.1 Traditional SIR model

The SIR (Susceptible - Infected - Removed) model is a mathematical model
commonly used to describe the dynamics of epidemics (Kermack and McKendrick,
1927). It divides a population of $N$ individuals into 3 disjoint groups: (S)
Susceptible = the people who have not been infected, but can potentially get
the disease; (I) Infected = the people who are able to transmit the virus to
the susceptible people; (R) Removed = the people who recovered from the
infection, or who died. This model assumes that the population size is
constant ($N=S+I+R$) and is described by the set of differential equations:

$\begin{gathered}\frac{dS}{dt}\ =\
-\frac{\beta\,S(t)\,I(t)}{N},\qquad\frac{dI}{dt}\ =\
\frac{\beta\,S(t)\,I(t)}{N}\ -\ \gamma\,I(t),\qquad\frac{dR}{dt}\ =\
\gamma\,I(t)\end{gathered}$ (1)

where $S(t),I(t)$ and $R(t)$ represents the number of people in the
susceptible, infected and removed compartments, respectively, at time $t$.
$\beta$ represents the transmission rate (how quickly people move from
susceptible to infected), and $\gamma$ represents the removed rate (how
quickly infected people recover or die).

Following Vega et al. (2022), we approximate the differential equations of the
SIR model with their discrete counterpart:

$\begin{array}[]{llrr}S_{t}\ \ =&S_{t-1}&-\
\displaystyle\beta\frac{S_{t-1}\,I_{t-1}}{N}&\\\\[5.59721pt] I_{t}\ \
=&I_{t-1}&+\ \displaystyle\beta\frac{S_{t-1}\,I_{t-1}}{N}&-\
\gamma\,I_{t-1}\\\ R_{t}\ \ =&R_{t-1}&&+\ \gamma\,I_{t-1}\end{array}$ (2)

and estimate the optimal parameters $\beta^{*}$ and $\gamma^{*}$ by solving
the optimization problem of Equation 3:

$\begin{gathered}(\beta^{*},\
\gamma^{*})\quad=\quad\operatorname*{arg\,min}_{\beta,\gamma}\sum_{t=1}^{n}\left(I_{t}-\left(\beta\,\frac{S_{t-1}I_{t-1}}{N}\
-\ \gamma\,I_{t-1}\ +\ I_{t-1}\right)\right)^{2}\end{gathered}$ (3)

where $S_{t},I_{t},R_{t}$ are the number of individuals in the susceptible,
infected and removed compartments, respectively, at time $t$;
$S_{t-1},I_{t-1},R_{t-1}$ represent the number individuals in each compartment
at time $t-1$, $N\ =\ S_{t}+I_{t}+R_{t}$ is the time-invariant total number of
people, and $n$ represents the total number of points in the time-series used
for training. In our experiments every time-step spans a week.

Note that Equation 3 minimizes the error between the observations and
predictions over just the number of people infected, but not over the numbers
of susceptible or removed. This is because our purpose is to make predictions
over the infections, whose number is several orders of magnitude smaller than
the number of susceptibles.

###  2.2 Dictionary-based models

A dictionary, $D\in\mathbb{R}^{S\times T}$, is a matrix that contains $S$
time-series, each with $T$ time-points. We can then reconstruct an observed
time-series (e.g., the number of people infected of COVID-19) as a linear
combination of the elements in the dictionary:

$x=D^{T}\theta$ (4)

where $\theta\in\mathbb{R}^{S}$ contains the coefficient of each element in
the dictionary. Figure 1 depicts a simple example with a dictionary, $D$, that
contains 4 time-series (shown in orange, on the left). The vector of weights
$\theta=\begin{bmatrix}2&0&4&0\end{bmatrix}^{T}$ indicates the coefficients of
the linear combination of the signals in the dictionary used to reconstruct
the signal $x$ (shown in blue, on the right).

We can think of the entries in the dictionary as candidate trajectories
describing the number of people infected over time for a specific sub-
population. These entries are computed a priori and should include plausible
trajectories that we might expect to observe. There are different ways of
computing these curves. For example, we can create a set of Gaussian
probability density functions with different parameters for the mean,
variance. Alternatively, we can simulate a set of SIR models with different
parameters (initial number of people in each compartment at a specified time,
and parameters $\beta$ and $\gamma$), and then include in the dictionary the
curve corresponding to the number of infected people for that specific sub-
population, $I^{r}=[I^{r}_{0},\,I^{r}_{1},\,\dots,\,I^{r}_{T}]$.

We can then estimate the total number of infected people at each time
$I=[I_{0},\,I_{1},\,\dots,\,I_{k}]$ as a linear combination of the entries in
the dictionary $(I_{t}\ =\ \sum_{r=1}^{M}\theta_{r}\,I^{r}_{t})$, which we can
compare with the observed values. The vector of weights $\theta$ determines
which curves should be included in the linear combination (when
$\theta_{r}\neq 0)$, as well as their coefficients.

Given an observed time-series, $x$, and a dictionary with candidate
trajectories, $D$, we can estimate the parameter $\theta$ of Equation 4 by
solving the optimization problem:

$\begin{gathered}\theta^{*}\quad=\quad\operatorname*{arg\,min}_{\theta}||x\,-\,D^{T}\theta||^{2}_{2}\
\ +\ \ \lambda\,||\theta||^{2}_{2}\\\ \mbox{s.t.}\ \qquad\theta_{i}\geq
0\end{gathered}$ (5)

where $\lambda$ is a regularization parameter used to decrease the chance of
overfitting. Note that Equation 5 requires _non-negative_ coefficients for the
candidate curves in the dictionary, as the entries in the dictionaries
represent plausible trajectories over the number of infected people for
different sub-populations, which can only be added, but not subtracted – i.e.,
the number of people infected in a sub-population cannot be negative.

Equation 5 is a standard convex optimization problem with linear constraints.
The indexes of the non-zero coefficients of the parameter $\theta$ represent
the indexes of the sub-populations included in the model. Note that the number
of non-zero indexes is automatically determined by the algorithm, so there is
no need to manually identify the number of sub-populations to be included.

For our experiments, we compared the performance of two different
dictionaries: (1) $D_{G}$, whose entries are normalized probability density
functions (PDF) of a Gaussian. We normalize each curve by dividing the entire
PDF by its maximum. The dictionary of Gaussians, $D_{G}$, contains 390
trajectories, obtained by using all possible combinations of a set of means
$\mu\in\\{0,2,4,\dots,52\\}$ and a set of variances
$\sigma^{2}\in\\{1^{2},3^{2},5^{2},\dots,29^{2}\\}$. Intuitively, the mean
$\mu$ represents the week at which we expect the peak of the number of
infected people in a sub-population, while the standard deviation $\sigma$
represents the speed at which the number of infected people increases and
decreases. If required, we could make the Gaussian distributions non-symmetric
by controlling the skewness of the distribution (Azzalini and Capitanio,
1999).

(2) We also consider the $D_{SIR}$ dictionary, whose entries are each
generated by taking the number of people in the infected compartment of a SIR
model initialized with different parameters. For each, we need to specify the
parameters $\beta$ and $\gamma$ of the traditional SIR model, and also the
initial number of susceptible, infected and removed at the beginning of the
simulation. We also need to create time-shifts in the traditional SIR model –
_e.g. , _ to specify that the number of infected people will start increasing
after 10 weeks – i.e., become non-zero after week 10.

To do this, our model specifies that the initial number of people in the
infected and recovered compartments is zero, and then modifies the discrete-
time SIR model of Equation 2 as follows:

$\begin{array}[]{llrrr}S_{t}&=\ S_{t-1}&-\ \ S_{t-1}\,I_{t-1}\,\beta&-\ C\
\delta(t-k)\\\ I_{t}&=\ I_{t-1}&+\ \ S_{t-1}\,I_{t-1}\,\beta&+\ C\
\delta(t-k)&-\ \gamma\,I_{t-1}\\\ R_{t}&=\ R_{t-1}&&&+\ \
\gamma\,I_{t-1}\end{array}$ (6)

where $\delta(\cdot)$ is the Kronecker delta function, $k$ is the ‘start time’
for this sub-population, and $C$ is a parameter that specifies how many people
will be incorporated into this infected compartment at time $k$.

Using this model of Equation 6, we created a dictionary $D_{SIR}$ with 546
curves obtained by using all possible combinations of the parameters
$S_{0}\in\\{10^{4},10^{5},10^{6}\\},\ \beta\in\\{0.3,0.4,\dots 0.9\\},\
k\in\\{0,2,\dots,50\\}$, where we also set $\gamma=0.5,\ I_{0}=0,\ R_{0}=0,\
C=100$. Similarly to the dictionary with Gaussians, we normalized the curves
of the infected population by dividing over the maximum value of each curve,
so the peak of each curve is always one.

###  2.3 Fixed size mixture of fittable curves model

The dictionary-based models have the advantage of automatically finding the
number of sub-populations to include in the model, as well as the coefficient
associated with each sub-population; however, it might be useful to instead
just specify a priori the number, $M$, of sub-populations that we want to
include but here allow the system to fit the associated parameters.
Controlling the number of non-zero coefficients in Equation 5 requires a
0-norm regularization term, which is intractable for large models (e.g.,  when
there are a large number of curves in a dictionary).

Instead, we can use different parametric models to fit the observed data. For
example, given a sequence of T time-points, $x=[0,1,\dots,T]$:

$f(\,x,\,\theta,\,\mu,\,\sigma\,)\ \ =\ \ \theta_{0}\ +\
\sum_{m=1}^{M}\theta_{m}\
\exp\left(-\frac{(x-\mu_{m})^{2}}{2\sigma_{m}^{2}}\right)$ (7)

Then, given an observation of the number of infected people over time,
$I=[I_{1},I_{2},\dots,I_{T}]$, we find the parameters of the Gaussian mixture
that solve the optimization problem:

$\theta^{*},\mu^{*},\sigma^{*}\quad=\quad\operatorname*{arg\,min}_{\theta,\mu,\sigma}||y\
-\ f(x,\theta,\mu,\sigma)||_{2}^{2}$ (8)

Alternatively, we can model the number of infected people by taking the
Infected compartment of a mixture-of-$M$ fittable SIR models. Each SIR model
in the mixture follows the set of Equations 6. Therefore, the number of people
infected at time $t$ is given by:

$f(S_{t-1},I_{t-1},\beta,\gamma,C,k)\ =\
\sum_{m=1}^{M}S_{t-1}^{m}I_{t-1}^{m}\beta^{m}+C^{m}\delta(t-k^{m})-\gamma^{m}I_{t-1}^{m}+I_{t-1}^{m}$
(9)

There are two important things to note in Equation 9. First, each of the terms
$S_{t-1},I_{t-1},\beta,\gamma,C$ and $k$ is a vector of length $M$, whose
$i$-th entry represents the value of that parameter for the $i$-th sub-
population. Second, we do not observe the values of any of these vectors, but
instead only observe the aggregate quantity $I_{t}\ =\ \sum_{m}I_{t,m}$. That
is, we observe the total reported number of infected people, but we do not
know how many infected are in each sub-population.

Since Equation 9 is recursive, we can compute a time-series of length T if we
have the values $S_{t=0},\ I_{t=0},\ \beta,\ \gamma,\ C$ and $k$. Given a
sequence of T observations $I=[I_{1},I_{2},\dots,I_{T}]$ that represent the
number of reported infections over time, we assume that $I_{t=0}=[0]^{M}$ and
we can find the remaining parameters by solving the optimization problem:

$S_{0}^{*},\beta^{*},\gamma^{*},C^{*},k^{*}\ =\
\operatorname*{arg\,min}_{S_{0},\beta,\gamma,C,k}||y-\hat{y}(S_{0},\beta,\gamma,C,k,T)||_{2}^{2}$
(10)

where $\hat{y}(S_{0},\beta,\gamma,C,k)$ is a time-series of length T obtained
by recursively applying Equation 9 $T$ times.

Although the optimization problems of Equations 8 and 10 are not convex, we
can find a locally optimal solution using optimization techniques such as
generalized simulated annealing (Xiang and Gong, 2000). This optimization
technique requires specifying a tuple (lower boundj, upper boundj) that sets
the range of possible valid values for the $j^{th}$ variable to be optimized.
For optimizing Equation 8, we used $\mu\in(0,50),\sigma\in(1,6)$ and
$\theta\in(0,3\times 10^{5})$, and for optimizing Equation 10, we used
$S\in(0,10^{8})$, $\beta\in(0,1)$, $\gamma\in(0,1)$, $C\in(0,10^{3})$, and
$k\in(0,50)$. For both equations we set $M=3$ to illustrate the performance of
the algorithm. In a real application, the number of curves will change for
every geographical location; however, note that the number of parameters in
the search space increases exponentially with $M$. Therefore, we only
recommend this method for small values of $M$. For larger values, the
dictionary approach is more efficient.

Figure 3: Reconstruction of the number of weekly number of people infected in
Canada during 52 weeks, from 30 July 2020 to 29 July 2021. The blue line
represents the observed data, while the orange line represents the
reconstruction obtained by using a mixture of 2 Gaussians.

##  3 Results

We evaluated the performance of the methods described in Section 2, in two
different tasks: (T1) Accurate modelling of the number of infected people, and
(T2) Forecasting of the number of infected people 1 to 4 weeks in advance. For
T1 (modelling the pandemic) we use the same data for training and evaluation,
but T2 (Forecasting) requires a training set for learning the parameters, and
a disjoint set to assess the performance. Here, for each $t$, we predicted the
number at time $t+1$ based on the data observed from time $0,\dots,t$. For
both tasks, we evaluated the learned model using the mean absolute percentage
error (MAPE):

$\mbox{MAPE}\ \ =\ \ \frac{1}{n}\
\sum_{t=1}^{n}\left|\frac{A_{t}-F_{t}}{A_{t}}\right|$ (11)

where $n$ is the total number of time-points, $A_{t}$ is the actual value
observed at time $t$, and $F_{t}$ is the forecasted value at time $t$.

For both tasks, we used the publicly available COVID-19 Data Repository by the
Center for Systems Science and Engineering at Johns Hopkins University (Dong
et al., 2020), and followed the pre-processing pipeline suggested by Vega et
al. (2022). We restricted our analysis to the time period between 30 July 2020
and 29 July 2021 – selecting this time frame because it includes two or more
“COVID waves", and also because 187 countries were consistently reporting the
number of new people infected during that period.

###  3.1 Modeling task (T1)

Figure 4 shows the distribution of the MAPE over all 187 countries, for the
modeling task, for each of the models considered: the classical SIR model, the
dictionary models based on Gaussian distribution (resp., SIRs), as well as the
mixture models of Gaussian distributions (resp., SIRs). Table 1 provides the
corresponding numerical values.

Figure 4: Visualization of the distributions of the modeling task, in terms of mean absolute percentage error (MAPE), for the different modeling methods. Table 1: MAPE statistics of the different methods in the modeling task. | Mean |  Std | Min | 25% | 50% | 75% | Max
---|---|---|---|---|---|---|---
Gauss Dict | 3.5 |  4.7 | 0.1 | 1.1 | 1.9 | 4.1 | 38.5
SIR Dict | 7.5 |  13.3 | 0.1 | 1.8 | 3.8 | 7.7 | 108.9
SIR | 8.7 |  8.2 | 0.1 | 3.6 | 6.5 | 11.5 | 64.5
Mix Gauss | 11.3 |  13.9 | 0.2 | 3.3 | 7.0 | 12.7 | 76.0
Mix SIR | 13.7 |  16.3 | 0.3 | 4.6 | 10.5 | 16.5 | 121.2

To compute the MAPE, we use the entire time-series as an input for each of the
methods and learn the parameters that optimize Equations 3, 5, 8, or 10,
depending on the method selected. Then, we used the time-series
$x_{0},x_{1},\dots,x_{t}$, along with the learned parameters, to predict
$x_{t+1}$ for $t=1,\dots,52$. We applied this procedure individually to each
of the 187 countries in the dataset.

Figure 4 and Table 1 show that the dictionary based approaches not only have
the lowest median MAPE, but also the smallest variance. Both dictionary based
approaches show a better performance than the standard SIR model. Finally, the
mixtures of M fittable curves (for $M=3$) have a higher median error and
higher variance than traditional SIR models.

While both the dictionary-based methods and the mixture of $M$ fittable curves
compute mixtures of curves, they differ in that dictionary-based approaches
automatically decide on the number of curves to use, as well as their
respective coefficients. For example, Figure 5 shows the coefficients learned
by the algorithm based on Gaussian dictionaries; the non-zero coefficients,
which here is 32, are the number of selected sub-populations. By contrast, the
mixture of $M$ fittable curves approach requires the user to explicitly set
the number of mixtures to use.

Figure 5: Visualization of the weights learned by the Gaussian dictionary
approach in data from Canada.

Figure 6 depicts examples of curves generated, for modeling the number of
infected people over time, by each of the 5 methods, for 6 different
countries. We see that the Gaussian dictionary approach is able to accurately
model trajectories that are multimodal – i.e., it can properly model different
waves. It also can smooth trajectories that have high-frequency oscillations –
e.g., see Colombia. This smoothing effect is more visible for the dictionaries
based on SIR (see the behaviour of the first 15 weeks in Colombia, Egypt, and
South Korea).

While the traditional SIR model performs fairly well for uni-modal
trajectories – like Egypt, Finland, and Guatemala – it struggles with
multimodal time-series, such as Canada, Colombia and South Korea. When there
multiple waves, the model basically produces (an almost) delayed version of
the original data. This occurs because the SIR is unimodal by definition, so
most of the time the prediction for $x_{t+1}$ is the value very close to the
value observed at $x_{t}$.

Finally, the mixture of $M$ fittable curves have a higher error, on average,
than the other approaches. While these models do capture the general trend of
the data, as they have a limited number of curves (3 in this case), they are
unlikely to capture the minor variations of the original data. Also, note that
the mixture of SIR models always start at the minimum possible value, so it
fails to properly model cases where the minimum number of infected people is
not at the beginning of the training sequence (see the cases of Colombia and
Egypt).

Despite their inferior performance, the mixtures of M fittable curves might be
useful for understanding and explaining the dynamics of the number of infected
people over time. For example, in the cases of Canada and Colombia in Figure
6, it might be easier to explain the number of infected people in terms of 2
or 3 waves, respectively, than to give precise meaning to all the non-zero
coefficients depicted in Figure 5. Also, note that the “ground truth” is the
reported number of infections over time, which in general is only an
approximation of the actual number of infected people over time (which is in
fact unknown). Therefore, in some situations it might be desirable to model
the pandemic with a method that ignores the small changes in the time-series,
which might be just noise.

Figure 6: Examples of the resulting graphs produced by the different models in
the reconstruction task. The plots show the weekly number of infected people
in different countries.

Figure 7: Visualization of the distributions of the predictions 1-to-4 weeks
in advance, in terms of mean absolute percentage error (MAPE), for the
different predictions methods.

###  3.2 Forecasting task (T2)

Here, we compare 6 different approaches for forecasting the number of infected
people 1- to 4-weeks in advance: the classical SIR model, the dictionary
approaches based on Gaussian distributions and SIRs, the mixture of $M=3$
fittable SIR models and $M=3$ Gaussian models, and the “same as last observed
week” method (SLOW) method. Figure 6 plots the distribution of the MAPEs for
each approach for forecasting the number of infected people, 1- to 4-weeks in
advance; this is over all 187 countries. Their numerical values are shown in
Tables 2 to 5. As the values for the mixture of Gaussians are so much higher
than the other methods, we did not include them in Figure 7 for an easier
visualization.

Unlike the modeling task, which learns the model’s optimal parameters after
observing the _entire_ dataset, in the forecasting task, each approach
forecasts the number of infected people in time $t+1$, $x_{t+1}$ (and also
$x_{t+2},x_{t+3}$ and $x_{t+4})$ based on only the models learned from the
observations $x_{0},\dots,x_{t}$. We repeat this process iteratively from
$t=5$ to $t=48$.

Here, we included SLOW as a baseline method. While this method is very naive,
Vega et al. (2022) showed that this approach did quite well – indeed, better
than many other approaches to this task; see the results of the COVID-19
Forecasting Hub (Cramer et al., 2021), which is the official forecasting
project of the Center for Disease Control and Prevention (CDC) in the United
States.

Similar to the reconstruction task, the dictionary-based approaches outperform
the classical SIR model in terms of MAPE. Their variance is also the smallest
across the different methods used in this comparison and they are the only
approaches that perform better than the baseline model (SLOW). The mixture of
$M$ fittable curves, again, have the worst performance with the highest
variance.

As expected, the MAPE increases with longer forecasting horizons. While the
median MAPE of the Gaussian dictionary is 6% for 1-week forecasting, it
increases to 11.8%, 17.1%, and 21.7% when forecasting 2, 3, and 4 weeks ahead,
respectively. The other methods follow a similar behaviour. Figure 8 shows an
example of the predictions made by the different methods for the Canadian
dataset. Note how the dictionary based approaches tend to make more accurate
predictions than the SIR, SLOW, and mixture models.

Table 2: Prediction results 1 week in advance | Mean |  Std | Min | 25% | 50% | 75% | Max
---|---|---|---|---|---|---|---
Gauss Dict | 8.1 |  8.7 | 0.0 | 3.2 | 6.0 | 10.1 | 71.0
SIR Dict | 10.8 |  16.8 | 0.0 | 3.7 | 7.2 | 12.3 | 152.6
SIR | 15.9 |  16.2 | 0.0 | 6.3 | 12.0 | 21.2 | 128.5
Mix Gauss | 31.2 |  39.5 | 0.0 | 7.7 | 16.0 | 38.5 | 247.4
Mix SIR | 165.2 |  1350.7 | 0.2 | 6.1 | 11.7 | 23.8 | 18127.3
SLOW | 8.7 |  8.7 | 0.0 | 3.6 | 6.7 | 11.2 | 68.9
Table 3: Prediction results 2 weeks in advance | Mean |  Std | Min | 25% | 50% | 75% | Max
---|---|---|---|---|---|---|---
Gauss Dict | 14.8 |  14.5 | 0.0 | 6.5 | 11.8 | 19.1 | 128.5
SIR Dict | 16.5 |  21.7 | 0.0 | 6.6 | 11.5 | 20.3 | 174.9
SIR | 24.1 |  25.7 | 0.0 | 8.7 | 16.8 | 30.4 | 200.0
Mix Gauss | 228.1 |  492.2 | 0.1 | 22.2 | 54.9 | 173.8 | 3484.4
Mix SIR | 1E9 |  2E10 | 0.3 | 10.5 | 20.4 | 42.1 | 2E11
SLOW | 15.9 |  15.7 | 0.0 | 6.4 | 12.1 | 20.7 | 125.5
Table 4: Prediction results 3 weeks in advance | Mean |  Std | Min | 25% | 50% | 75% | Max
---|---|---|---|---|---|---|---
Gauss Dict | 20.8 |  20.4 | 0.0 | 9.3 | 17.1 | 26.9 | 182.8
SIR Dict | 22.3 |  27.2 | 0.0 | 9.1 | 16.5 | 26.8 | 229.1
SIR | 33.2 |  36.1 | 0.0 | 12.1 | 22.8 | 43.1 | 270.9
Mix Gauss | 1222.1 |  3381.4 | 0.8 | 38.6 | 111.9 | 519.4 | 23244.3
Mix SIR | 4E23 |  6E24 | 0.3 | 18.1 | 33.7 | 72.1 | 8E25
SLOW | 22.5 |  22.7 | 0.0 | 8.7 | 16.9 | 29.6 | 185.8
Table 5: Prediction results 4 weeks in advance | Mean |  Std | Min | 25% | 50% | 75% | Max
---|---|---|---|---|---|---|---
Gauss Dict | 25.4 |  23.2 | 0.0 | 11.8 | 21.7 | 33.5 | 220.5
SIR Dict | 27.3 |  29.7 | 0.0 | 11.9 | 21.5 | 33.6 | 223.0
SIR | 44.7 |  56.8 | 0.0 | 14.6 | 28.4 | 53.6 | 461.3
Mix Gauss | 3506.1 |  10909.3 | 2.1 | 42.1 | 153.1 | 1016.6 | 71112.5
Mix SIR | 4E52 |  5E53 | 0.4 | 27.4 | 50.7 | 119.7 | 7E54
SLOW | 28.3 |  27.7 | 0.0 | 11.0 | 21.1 | 37.3 | 234.1

Figure 8: Examples of the curves predicted by the different methods for Canada
1-to-4 weeks in advance. (GT = ground truth).

##  4 Discussion

Despite their popularity, SIR models and their extensions are usually
unimodal, which limits their usability in time-series with several infection
waves. These models also assume an homogeneous population. These assumptions
appear invalid in many real life scenarios. Here we presented a method to
model heterogeneous populations as a linear combination of several independent
sub-populations, where the details of these sub-populations (e.g., which
person is in which sub-population) is unknown. This approach naturally allows
us to model different waves of the pandemic. We offer two alternatives to set
the number of subpopulations to include in the model: to manually define the
number of sub-populations (mixture of $M$ fittable curves model), or to use an
optimization algorithm that automatically determines how many sub-populations
to use, among a set of candidate trajectories (Dictionary based approach).

Dictionary based approaches showed a better performance, in terms of median
and variance of the MAPE, on both the modeling and the forecasting task;
however their performance strongly depends on the curves used to create the
dictionary. Also, as shown in Figure 5, the optimization algorithm selects a
relatively high number of curves. Although we could interpret each of these
curves as a sub-population, it would be difficult to identify the exact nature
of them. Of course, we never claim that any of our approaches are learning
causality, so even when they accurately model the dynamics of the number of
infected people over time, there is no way of determining if in fact those are
the real number of sub-populations present in the data. Despite this
limitation in interpretability, these models proved to be more accurate than
the rest in the short-term forecasting task.

The mixture-of-$M$ curves, on the other hand, had the worst performance, being
even lower than traditional single SIR models. Nevertheless, the mixture-
of-$M$ curves might prove useful when the objective is to describe the general
trend of the number of infected people over time. Figure 6 shows that the
mixture models ignore the high-frequency oscillations in the data, identifying
only the “main waves” present in the observed time-series. Since users have
control over the number of sub-populations to model, it might be easier to
assign a meaning to the different sub-populations in the model.

There have been several earlier attempts to include sub-populations in the SIR
models; however, in most of these attempts the infections in the different
sub-populations “start” at the same time. Here, we propose the use of the
Kronecker delta in the SIR equations to allow for the modelling of a
combination of SIR’s that start at different points in time.

Importantly, the dictionaries and the mixture-of-$M$ fittable curves are
general methods for building mixtures of models. We presented examples based
on Gaussian distributions and SIR models, and showed their effectiveness over
187 countries; however, this is a general approach that is compatible with any
other parametric model, such as SEIR and other variants of the classical SIR
model. We expect that this approach can be used to model not only COVID-19,
but other infectious diseases as well. To facilitate the exploration of this
approach, the code for reproducing these results and for implementing the
different approaches here presented is publicly available at
https://github.com/rvegaml/SubpopulationModels

## References

  * Albi et al. [2021] G. Albi, L. Pareschi, and M. Zanella.  Control with uncertain data of socially structured compartmental epidemic models.  _Journal of Mathematical Biology_ , 82(7):1–41, 2021.
  * Azzalini and Capitanio [1999] A. Azzalini and A. Capitanio.  Statistical applications of the multivariate skew normal distribution.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 61(3):579–602, 1999.
  * Bartolucci and Farcomeni [2022] F. Bartolucci and A. Farcomeni.  A spatio-temporal model based on discrete latent variables for the analysis of covid-19 incidence.  _Spatial Statistics_ , 49:100504, 2022.
  * Calvetti et al. [2020] D. Calvetti, A. P. Hoover, J. Rose, and E. Somersalo.  Metapopulation network models for understanding, predicting, and managing the coronavirus disease covid-19.  _Frontiers in Physics_ , 8:261, 2020.
  * Contreras et al. [2020] S. Contreras, H. A. Villavicencio, D. Medina-Ortiz, J. P. Biron-Lattes, and Á. Olivera-Nappa.  A multi-group seira model for the spread of covid-19 among heterogeneous populations.  _Chaos, Solitons & Fractals_, 136:109925, 2020.
  * Cramer et al. [2021] E. Y. Cramer, Y. Huang, Y. Wang, E. L. Ray, M. Cornell, J. Bracher, A. Brennen, A. J. Castro Rivadeneira, A. Gerding, K. House, D. Jayawardena, A. H. Kanji, A. Khandelwal, K. Le, J. Niemi, A. Stark, A. Shah, N. Wattanachit, M. W. Zorn, N. G. Reich, and U. C.-. F. H. Consortium.  The united states covid-19 forecast hub dataset.  _medRxiv_ , 2021.  doi: 10.1101/2021.11.04.21265886.
  * Dolbeault and Turinici [2020] J. Dolbeault and G. Turinici.  Heterogeneous social interactions and the covid-19 lockdown outcome in a multi-group seir model.  _Mathematical Modelling of Natural Phenomena_ , 15:36, 2020\.
  * Dolbeault and Turinici [2021] J. Dolbeault and G. Turinici.  Social heterogeneity and the covid-19 lockdown in a multi-group seir model.  _Computational and Mathematical Biophysics_ , 9(1):14–21, 2021.
  * Dong et al. [2020] E. Dong, H. Du, and L. Gardner.  An interactive web-based dashboard to track covid-19 in real time.  _The Lancet infectious diseases_ , 20(5):533–534, 2020.
  * Duan and Jin [2022] M. Duan and Z. Jin.  The heterogeneous mixing model of covid-19 with interventions.  _Journal of theoretical biology_ , page 111258, 2022.
  * Ellison [2020] G. Ellison.  Implications of heterogeneous sir models for analyses of covid-19.  Technical report, National Bureau of Economic Research, 2020.
  * Gerasimov et al. [2021] A. Gerasimov, G. Lebedev, M. Lebedev, and I. Semenycheva.  Covid-19 dynamics: a heterogeneous model.  _Frontiers in Public Health_ , 8:558368, 2021.
  * Kermack and McKendrick [1927] W. O. Kermack and A. G. McKendrick.  A contribution to the mathematical theory of epidemics.  _Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character_ , 115(772):700–721, 1927.
  * Kong et al. [2022] L. Kong, M. Duan, J. Shi, J. Hong, Z. Chang, and Z. Zhang.  Compartmental structures used in modeling covid-19: a scoping review.  _Infectious diseases of poverty_ , 11(1):1–9, 2022.
  * Nepomuceno et al. [2016] E. G. Nepomuceno, R. H. C. TAKAHASHI, and L. A. AGUIRRE.  Individual-based model (ibm): an alternative framework for epidemiological compartment models.  _Brazilian Journal of Biometrics_ , 34(1):133–162, 2016.
  * Sarker et al. [2022] A. Sarker, A. Jadbabaie, and D. Shah.  Unifying epidemic models with mixtures.  _arXiv preprint arXiv:2201.04960_ , 2022.
  * Shi et al. [2022] R. Shi, H. Zhang, and Q. Morris.  Pan-code: Covid-19 forecasting using conditional latent odes.  _Journal of the American Medical Informatics Association_ , 2022.
  * Tan et al. [2021] Y. Tan, D. Iii, M. Ndeffo-Mbah, and U. Braga-Neto.  A stochastic metapopulation state-space approach to modeling and estimating covid-19 spread.  _Mathematical Biosciences and Engineering: MBE_ , 18(6):7685–7710, 2021.
  * Vega et al. [2022] R. Vega, L. Flores, and R. Greiner.  Simlr: Machine learning inside the sir model for covid-19 forecasting.  _Forecasting_ , 4(1):72–94, 2022.
  * Xiang and Gong [2000] Y. Xiang and X. Gong.  Efficiency of generalized simulated annealing.  _Physical Review E_ , 62(3):4473, 2000.
  * Zachreson et al. [2022] C. Zachreson, S. Chang, N. Harding, and M. Prokopenko.  The effects of local homogeneity assumptions in metapopulation models of infectious disease.  _Royal Society Open Science_ , 9(7):211919, 2022\.


# Enhancement in Breaking of Time-reversal Invariance in the Quantum Kicked
Rotor

Ramgopal Agrawal<EMAIL_ADDRESS>Akhilesh Pandey
<EMAIL_ADDRESS>Sanjay Puri<EMAIL_ADDRESS>School of Physical
Sciences, Jawaharlal Nehru University, New Delhi – 110067, India.

###### Abstract

We study the breaking of time-reversal invariance (TRI) by the application of
a magnetic field in the quantum kicked rotor (QKR), using Izrailev’s finite-
dimensional model. There is a continuous crossover from TRI to time-reversal
non-invariance (TRNI) in the spectral and eigenvector fluctuations of the QKR.
We show that the properties of this TRI $\rightarrow$ TRNI transition depend
on $\alpha^{2}/N$, where $\alpha$ is the chaos parameter of the QKR and $N$ is
the dimensionality of the evolution operator matrix. For $\alpha^{2}/N\gtrsim
N$, the transition coincides with that in random matrix theory. For
$\alpha^{2}/N<N$, the transition shows a marked deviation from random matrix
theory. Further, the speed of this transition as a function of the magnetic
field is significantly enhanced as $\alpha^{2}/N$ decreases.
(To appear in Physical Review E.)

##  I Introduction

The quantum kicked rotor (QKR) has been a prototype model for quantum chaos.
It was first introduced by Casati et al. Casati  _et al._ (1979) as a quantum
analog of the classical kicked rotor (CKR). An important paradigm of quantum
chaos states that the spectral properties of quantum chaotic systems should
follow the predictions of random matrix theory (RMT) Bohigas  _et al._ (1984);
Brody  _et al._ (1981); Mehta (2004); Haake (2010); Beenakker (1997); Guhr
_et al._ (1998). For time-independent systems, Gaussian ensembles of random
matrices are appropriate. On the other hand, for time-periodic systems,
circular random matrix ensembles are more useful. Izrailev Izrailev (1986,
1990) introduced a model for the evolution operator of the QKR in
$N$-dimensional Hilbert space, and found that the quasi-energy fluctuations
agree with RMT. The discrete symmetries present in the QKR classify the
spectral statistics according to the universality classes of RMT. In case the
QKR maintains time-reversal invariance (TRI), the statistics is described by
circular orthogonal ensembles (COE). The TRI is broken by applying a strong
magnetic field, and the spectral statistics then becomes consistent with that
of circular unitary ensembles (CUE). However, at moderate perturbation
strengths, it is found that the spectral fluctuations exhibit intermediate
statistics, and the TRI is said to be partly broken Pandey  _et al._ (1993);
Lenz and Haake (1991).

Over the decades, breaking of TRI in quantum chaotic systems has received
special attention in the literature Beenakker (1997); Dupuis and Montambaux
(1991); Bohigas  _et al._ (1995); Pluhar  _et al._ (1995). According to the
von Neumann-Wigner theorem, TRI-breaking is accompanied by increased level
repulsion in the spectrum Mehta (2004); Haake (2010). Berry and Robnik Berry
and Robnik (1986) studied chaotic billiards with a single line of magnetic
flux as a perturbation, and confirmed this level repulsion. The eigenvector
properties also change during the breaking of TRI Rosenzweig  _et al._ (1968);
French  _et al._ (1985); Fal’ko and Efetov (1994); Chung  _et al._ (2000).
When TRI is preserved, the absolute-squared eigenvector component has a
$\chi^{2}_{1}$-distribution. This becomes a $\chi^{2}_{2}$-distribution when
TRI is completely broken Brody  _et al._ (1981); Rosenzweig  _et al._ (1968).
(We will provide the functional forms of the $\chi^{2}_{1}$\- and
$\chi^{2}_{2}$-distributions later.) In the crossover regime, the distribution
is intermediate to $\chi^{2}_{1}$ and $\chi^{2}_{2}$. This TRI to time-
reversal non-invariance (TRNI) transition is driven by a control parameter
(say $\lambda$), which measures the strength of the applied field. The
influence of $\lambda$ depends on the dimensionality $N$. In the semi-
classical limit ($N\rightarrow\infty$), the transition is discontinuous in
$\lambda$. One identifies a local transition parameter, which smoothens the
transition for systems with large $N$. This parameter is further useful in
comparing the transition in different systems. It depends on the near-diagonal
matrix elements of the TRI-breaking perturbation and the mean level spacing of
the unperturbed spectrum French  _et al._ (1985, 1988a). In the context of
classical RMT, Pandey and Mehta Mehta and Pandey (1983); Pandey and Mehta
(1983) have investigated the spectral transition from Gaussian orthogonal
ensembles (GOE) to Gaussian unitary ensembles (GUE). Using Dyson’s Brownian
motion model Mehta (2004); Dyson (1962), they studied the continuous crossover
from GOE $\rightarrow$ GUE. Later, Pandey and Shukla Pandey and Shukla (1991);
Pandey (2004) extended the idea of Brownian motion to unitary ensembles and
studied the COE $\rightarrow$ CUE transition.

In the present paper, we investigate the effect of the TRI $\rightarrow$ TRNI
transition on the spectral and eigenvector fluctuations of the QKR. We report
results from a detailed numerical study of this transition. The motivation for
this work comes from our recent proposal of finite-range Coulomb gas (FRCG)
models for random matrix ensembles, and their applicability to a finite-
dimensional matrix model of the QKR Pandey  _et al._ (2017, 2020); Kumar  _et
al._ (2020). The CKR is characterized by a parameter $\alpha$, with the system
being highly chaotic for $\alpha\gg 1$. The relevant parameter in the QKR is
$\alpha^{2}/N$, where $N$ is the size of the evolution operator matrix Casati
_et al._ (1990); Pandey  _et al._ (2017). In our recent work Pandey  _et al._
(2017, 2020); Kumar  _et al._ (2020), we have proposed that the spectral
statistics of the QKR with chaos parameter $\alpha$ is described by an FRCG
model of range $d=\alpha^{2}/N$. These FRCG models are generalizations of the
Dyson’s Coulomb gas model to the case where eigenvalues have a restricted
interaction range.

Let us next consider the TRI $\rightarrow$ TRNI transition in the QKR. For
$\alpha^{2}/N\gtrsim N$, this is described well by the COE $\rightarrow$ CUE
transition of RMT Pandey  _et al._ (1993); Shukla and Pandey (1997); Shukla
(1996). This is because the TRI-breaking perturbation, when represented in the
diagonal basis of the TRI-preserving part, belongs to the Gaussian RMT
ensembles. However, when $\alpha^{2}/N<N$, the perturbation is no longer
completely random. There has been almost no investigation of this crossover
when $\alpha^{2}/N<N$. This is the open problem which we address in this
paper. As we will show later, the perturbation matrix takes the form of a
power-law random band matrix (PRBM) Mirlin  _et al._ (1996); Mirlin (2000),
whose bandwidth depends on the parameter $\alpha^{2}/N$. The significant
perturbation matrix elements lie inside the band, and determine the nature of
the transition.

The main observations of our study are as follows: (a) For small values of
$\alpha^{2}/N$, the TRI $\rightarrow$ TRNI transition in the spectral
fluctuations of QKR does not obey the usual COE $\rightarrow$ CUE transition.
Further, the speed of the transition is enhanced with the decrease in
$\alpha^{2}/N$. (b) The transition in one-point eigenvector statistics also
does not obey the COE $\rightarrow$ CUE transition for any fixed value of
$\alpha^{2}/N<N$. As in the case of spectral fluctuations, the transition is
enhanced with decrease in $\alpha^{2}/N$.

This paper is organized as follows. In Sec. II, we discuss TRI-breaking
transitions in RMT using a Brownian motion model and introduce the local
transition parameter. In Sec. III, we briefly review the QKR and discuss its
spectral properties. In Sec. IV, we study the matrix elements of the TRI-
breaking perturbation in the QKR. The transition in spectral fluctuations of
the QKR is studied in Sec. V. In Sec. VI, the one-point eigenvector
fluctuations are analyzed. Finally, in Sec. VII, we summarize this paper and
discuss our results.

##  II  Breaking of time-reversal invariance: Local Transition Parameter

The TRI-breaking transition in time-periodic quantum systems is driven by an
external magnetic field of strength $\lambda$. The time-evolution operator of
the system can be expressed as

$\mathcal{U}(\lambda)=\mathcal{U}_{0}~{}{\rm e}^{i\lambda V},$ (1)

where $\mathcal{U}_{0}$ is the unperturbed part and $V$ is the TRI-breaking
perturbation.

It is known that quantum systems possessing TRI and/or exact symmetries are
often modeled by the stationary ensembles of classical RMT. In such systems,
the partial violation of TRI and/or an exact symmetry results in a crossover
between local fluctuation properties. The system, during this crossover,
follows intermediate Brownian motion ensembles (BME) Dyson (1962); Pandey and
Shukla (1991). In the context of quantum maps, which we study in this paper,
the relevant ensembles are the unitary BME. These are parameter-dependent
ensembles of $N$-dimensional unitary matrices $U(\tau)$, where $\tau$ plays
the role of fictitious time. The evolution of $U(\tau)$ under Brownian motion
is described as follows:

$U(\tau+\delta\tau)=U(\tau)~{}{\rm e}^{i\sqrt{\delta\tau}~{}V(\tau)}.$ (2)

Here, $\delta\tau$ is the infinitesimal time increment. At each $\tau$, a
random perturbation is introduced via an $N$-dimensional Hermitian random
matrix $V(\tau)$. In Dyson’s model for circular ensembles, the perturbation
$V(\tau)$, when represented in $U(\tau)$-diagonal basis, belongs to Gaussian
ensembles of random matrices. On the other hand, in FRCG models, $V(\tau)$
takes the form of band matrices Pandey  _et al._ (2017, 2020); Kumar  _et al._
(2020). The Brownian motion of matrix elements in Eq. (2) implies that the
quasi-energies of the parameter-dependent ensembles execute Brownian motion
under two-body Coulomb gas repulsion.

Under a nonlinear rescaling between the fictitious time $\tau$ and the TRI-
breaking parameter $\lambda$ with $\tau=\lambda^{2}$ and
$\delta\tau=(\delta\lambda)^{2}$, Eq. (2) is appropriate to model TRI-breaking
in the quantum system. As stated earlier, the transition in terms of $\lambda$
is dependent on the dimension $N$ and is discontinuous for
$N\rightarrow\infty$ Berry and Robnik (1986); Mehta and Pandey (1983). To
obtain a smooth transition, one has to parameterize the transition in terms of
a local transition parameter. This parameter is obtained from perturbation
theory French  _et al._ (1985, 1988a) as follows:

$\Lambda=\frac{\lambda^{2}~{}v^{2}}{D^{2}}.$ (3)

Here, $D$ is the mean level spacing of the quasi-energy spectrum in the
unperturbed case. Further, $v^{2}$ denotes the variance of near-diagonal
matrix elements of the TRI-breaking perturbation $V(\lambda)$ in the diagonal
basis of the TRI-part. The transition in nearest-neighbor fluctuations is
complete when $\lambda v$ is nearly equal to the mean level spacing $D$, i.e.,
$\Lambda\sim\mathcal{O}(1)$. Eq. (3) demonstrates the strong dependence of the
transition on the perturbation matrix elements.

##  III  Quantum Kicked Rotor

In the presence of a magnetic field of strength $\lambda$, the Hamiltonian for
the CKR is given by

$H=\frac{(p-\lambda)^{2}}{2}+\alpha\cos(\theta+\theta_{0})\sum_{n=-\infty}^{\infty}\delta(t-nT),$
(4)

where $\theta$ and $p$ are position and momentum variables. The parameter
$\alpha$ is the kicking strength, which introduces non-integrability in the
system. Therefore, it is known as the chaos parameter. The period between two
consecutive kicks is denoted as $T$. The time-reversal
($p\rightarrow-p,t\rightarrow-t$) and parity
($\theta\rightarrow-\theta,p\rightarrow-p$) invariance are broken by
parameters $\lambda$ and $\theta_{0}$, respectively.

The corresponding quantum dynamics can be understood in terms of a quantum
map, known as the QKR-map, which describes the evolution of a state
$|\psi_{t}\rangle$ as follows:

$|\psi_{t+T}\rangle=\mathcal{U}|\psi_{t}\rangle.$ (5)

Here, $\mathcal{U}$ is the evolution operator defined over one period of time.
It can be written as $\mathcal{U}=BG$, where

$B=\exp\left[-\frac{i\alpha\cos(\theta+\theta_{0})}{\hbar}\right]$ (6)

describes instantaneous kicking. Further,

$G=\exp\left[-\frac{iT(p-\lambda)^{2}}{2\hbar}\right]$ (7)

corresponds to the free motion of the rotor between two consecutive kicks. The
position and momentum variables $\theta$ and $p$ are now operators.

When dealing with TRI-breaking in quantum maps, the symmetric form of
$\mathcal{U}$ is more convenient to work with. This is related to other
representations by a unitary transformation, i.e.,

$\mathcal{U}=B^{1/2}GB^{1/2}.$ (8)

The chosen form of $\mathcal{U}$ is manifestly TRI at $\lambda=0$. For any
non-zero $\lambda$, $\mathcal{U}$ is TRNI. Thus, the momentum operator $p$ is
the TRI-breaking perturbation in the system. If we write $\mathcal{U}$ in a
basis in which $B$ is diagonal, i.e., position basis (see Appendix A for
details), we obtain the matrix elements:

$\begin{split}\mathcal{U}_{qq^{\prime}}=&\exp\left[-i\frac{\alpha}{2\hbar}\left\\{\cos\left(\theta_{q}+\theta_{0}\right)+\cos\left(\theta_{q}^{\prime}+\theta_{0}\right)\right\\}\right]\\\
&\times\sum_{l=-\infty}^{\infty}\exp\left[-i\frac{T}{2\hbar}(l\hbar-\lambda)^{2}\right]\exp\Big{[}i~{}l(\theta_{q}-\theta_{q}^{\prime})\Big{]}.\end{split}$
(9)

Since the momentum basis for the QKR is unbounded, the eigenvalues of the
infinite matrix $\mathcal{U}$ depend on the continuous parameters $\theta_{q}$
and $\theta_{q}^{\prime}$. For a discrete spectrum Manos and Robnik (2013);
Izrailev (1990), one is often interested in the level fluctuations. However,
there are infinite quasi-energies in the interval $(0,2\pi)$. Therefore, the
mean level spacing is zero, and the mean level density is infinite, i.e., the
fluctuations are Poisson. To observe level repulsion analogous to that in
bounded time-independent quantum chaotic systems Bohigas  _et al._ (1984);
Haake (2010), we need a finite-dimensional model of the QKR. Such a model was
first formulated by Izrailev Izrailev (1986). (We will discuss later that,
even in a finite-dimensional model, the spectral fluctuations can be Poisson.)
In the Izrailev model, an additional periodic boundary condition is imposed on
the momentum basis, i.e., the corresponding classical phase space becomes
bounded in both position and momentum
*[RegardlessofwhetherthephasespaceoftheCKRiscylindrical(boundedinonevariable)ortoroidal(boundedinbothvariables);
itsstatisticalpropertiesatlargekickingstrengthsareunchanged.See~][]ott1993chaos.
(There also exist other possible ways to approximate Eq. (9), e.g., truncation
of the infinite matrix $\mathcal{U}$ Izrailev (1990).) In Izrailev’s model,
the representation of $\mathcal{U}$ is chosen in a generic basis with a finite
number of levels (say, $N$). In the limit $N\rightarrow\infty$, this model
approaches the original infinite-dimensional QKR in Eq. (9). (See also the
recent work of Manos-Robnik Manos and Robnik (2013, 2015).)

To employ periodic boundary conditions in the momentum basis, we consider the
momentum eigenstates to obey the following condition: $|l\rangle=|l+N\rangle$,
where the label $l$ is an integer. This gives discrete values for the
eigenvalues of $\theta$:

$\theta_{m}=\frac{2\pi m}{N},~{}~{}~{}~{}~{}m=0,\pm 1,\pm 2,\ldots.$ (10)

Note that $\theta_{m}$ can take $N$ discrete values in the period
$\theta\rightarrow\theta+2\pi$. Now the position basis, as well as the
momentum basis, are periodic and both have a finite number of basis vectors
$N$. In this case, the evolution operator $\mathcal{U}$ in position basis is
the $N$-dimensional matrix:

$\begin{split}\mathcal{U}_{mn}=&\frac{1}{N}\exp\left[-i\frac{\alpha}{2}\left\\{\cos\left(\frac{2\pi
m}{N}+\theta_{0}\right)+\cos\left(\frac{2\pi
n}{N}+\theta_{0}\right)\right\\}\right]\\\
&\times\sum_{l=-N_{1}}^{N_{1}}\exp\left[-i\left(\frac{l^{2}}{2}-\lambda l-2\pi
l~{}\frac{(m-n)}{N}\right)\right].\end{split}$ (11)

Here, $N_{1}=(N-1)/2$, and the indices
$m,n=-N_{1},-N_{1}+1,\ldots,N_{1}-1,N_{1}$. Without loss of generality, we
have set $T=\hbar=1$. For $N\rightarrow\infty$, the model in Eq. (11)
converges to the infinite-dimensional model in Eq. (9) (see Appendix A).

For the CKR, we know that the dynamics depends entirely on $\alpha$. If
$\alpha\gg 1$, the stroboscopic map of the CKR is fully chaotic. The spectral
properties of the corresponding quantum map are expected to follow RMT
predictions Bohigas  _et al._ (1984), but this statement is only partly true.
As discussed above, the original QKR (which has an infinite-dimensional matrix
for the time-evolution operator $\mathcal{U}$) has Poisson quasi-energy
statistics for all $\alpha$. This occurs due to the localization of
eigenstates of the evolution operator in the momentum basis. This phenomenon
is known as dynamical localization, and is analogous to Anderson localization
in the 1-$d$ infinite lattice with disorder Fishman  _et al._ (1982); Anderson
(1958). The traces of dynamical localization also appear in the
$N$-dimensional model of the QKR with finite but large $N$. This is analogous
to Anderson localization in the 1-$d$ finite lattice Casati  _et al._ (1990).
Therefore, the quantum dynamics actually depends on two parameters: $\alpha$
and $N$. For applicability of RMT, apart from the requirement of chaotic
classical dynamics, the eigenstates of the evolution operator should also be
fully extended in the momentum basis. This second requirement depends on the
parameter $\alpha^{2}/N$ Pandey  _et al._ (2017); Casati  _et al._ (1990).
When $\alpha^{2}/N\gtrsim N$, all the eigenstates are extended. However, when
$\alpha^{2}/N$ decreases below $N$, the eigenstates start to become localized.
This gives rise to level crossing in the spectrum, and the spectral statistics
is then intermediate to Poisson and RMT.

We emphasize that the dynamics of the CKR is invariant under change in
parameters $\lambda$ and $\theta_{0}$ Pandey  _et al._ (1993); Shukla and
Pandey (1997). However, the corresponding quantum dynamics is not invariant.
The quantum dynamics has TRI for $\lambda=0$, but this is destroyed for
$\lambda\neq 0$. In the latter case, if $\theta_{0}=0$, TRI remains broken but
an additional TP-invariance holds (TP is the product of time-reversal and
parity operators). This TP-invariance of the QKR is broken by setting
$\theta_{0}\neq 0$. Further, when both $\lambda=0$ and $\theta_{0}=0$, the
parity operator commutes with the Hamiltonian and both can be diagonalized in
a common basis. On the other hand, for $\theta_{0}\neq 0$, the parity operator
no longer commutes with the Hamiltonian.

The spectral properties of the evolution operator can be analyzed in terms of
the number variance statistics $\Sigma^{2}(r)$, which is defined as the
variance of the number of quasi-energies $n$ in an interval of size $rD$,
where $D$ is the mean level spacing ($=2\pi/N$):

$\Sigma^{2}(r)=\overline{n^{2}}-\overline{n}^{2}.$ (12)

Here, the bar denotes the spectral or ensemble average. To numerically
calculate $\Sigma^{2}(r)$ from Eq. (12), we first find the quasi-energy
spectrum by diagonalizing $\mathcal{U}$ in Eq. (11). The quasi-energies are
then unfolded using the mean level density $N/2\pi$. We consider a range of
values for $\alpha^{2}/N$, with $\alpha\gg 1$ and a fixed $N=2001$.
Statistical averages are performed on an ensemble of $50$ independent spectra,
obtained by locally varying $\alpha$ in the range $[\alpha-5,\alpha+5]$. This
slight variation of $\alpha$ is adequate to ensure the statistical
independence of the spectra. The central value of $\alpha$ will be used to
label the corresponding spectra.

We first discuss the spectral properties of QKR when TRI is preserved
($\lambda=0$), and TRI is fully broken ($\lambda\simeq 1$). In both cases,
$\theta_{0}$ is fixed to a nonzero value: $\theta_{0}=\pi/(2N)$. In Fig. 1, we
plot $\Sigma^{2}(r)$ vs. $r$ for various values of $\alpha^{2}/N$. In Fig.
1(a), we show the case with $\lambda=0$. Since TRI is maintained, one expects
the number variance to follow COE predictions for
$\alpha^{2}/N\sim\mathcal{O}(N)$ Izrailev (1986), which is seen in the figure.
For small values of $\alpha^{2}/N<N$, one sees a clear deviation from COE in
the direction of the Poisson result. As $\alpha^{2}/N$ is reduced, the
deviation from COE is stronger at larger values of $r$. Note that, for any
value of $\alpha^{2}/N<N$, $\Sigma^{2}(r)$ attains the RMT-limit (COE, in this
case) in the interval size $r\ll\alpha^{2}/N$ and becomes Poisson-like
($\Sigma^{2}\sim r$) for $r\gg\alpha^{2}/N$ Pandey  _et al._ (2017); Kumar
_et al._ (2020).

In Fig. 1(b), we show analogous results for $\lambda=0.9$. In this case, TRI
is completely broken due to a strong magnetic field, and $\Sigma^{2}(r)$ vs.
$r$ for $\alpha^{2}/N\sim\mathcal{O}(N)$ coincides with the CUE result. For
$\alpha^{2}/N<N$, $\Sigma^{2}(r)$ is higher than the CUE value and approaches
the Poisson result, similar to Fig. 1(a). Comparing the data sets at a fixed
value of $\alpha^{2}/N$ in Figs. 1(a)-(b), we conclude that the spectrum is
more rigid in the TRNI case.

##  IV  Perturbation matrix elements

In the previous section, we discussed the spectral properties of the QKR for
zero and strong field strengths. However, when TRI is broken by the
application of moderate field strengths ($\lambda\ll 1$), the local
fluctuations exhibit intermediate statistics. We need to calculate the local
transition parameter $\Lambda$ (introduced in Sec. II) to characterize the
crossover obtained by varying field strength. It depends on the near-diagonal
matrix elements of the TRI-breaking perturbation, which is the momentum
operator $p$. For $\alpha^{2}/N\gtrsim N$, the off-diagonal elements of $p$ in
$\mathcal{U}(\lambda=0)$-diagonal basis are Gaussian random variables with
zero mean and uniform variance, i.e., $ip$ is a typical member of GOE Pandey
_et al._ (1993); Shukla and Pandey (1997). The variance of these matrix
elements depends linearly on $N$, $v^{2}\simeq\hbar^{2}N/12$. Then, the
transition parameter $\Lambda$ using Eq. (3) is given by

$\Lambda=\frac{\lambda^{2}N^{3}}{48\pi^{2}},$ (13)

where the mean level spacing of the quasi-energies of $\mathcal{U}(\lambda=0)$
is $D=2\pi/N$, and $\hbar$ is set to unity.

The above form of $\Lambda$ is valid for $\alpha^{2}/N\gtrsim N$. When
$\alpha^{2}/N$ decreases below $N$, the operator $p$ in the diagonal basis of
$\mathcal{U}(\lambda=0)$ is not fully random. In fact, it has a banded
structure with significant matrix elements lying inside a band around the
diagonal. The matrix elements are Gaussian distributed with mean zero, but the
variance depends on the distance from the diagonal.

In this section, we numerically investigate the perturbation matrix elements.
For a value of $\alpha^{2}/N$ (with fixed $N$ = 2001), we first write the
operator $p$ in $\mathcal{U}(\lambda=0)$-diagonal basis. The obtained matrix
of $p$ is imaginary and anti-symmetric, as the diagonal basis of $\mathcal{U}$
in Eq. (8) is time-reversal invariant at $\lambda=0$. To obtain good accuracy
in the statistics, we construct an ensemble of 50 independent matrices of $p$
as described in Sec. III. The matrix elements are denoted as $p_{ij}$. We have
$\overline{\overline{p_{ij}}}=0$, where the double bar denotes an ensemble
average and an average over elements at the same distance $L=|i-j|$ from the
diagonal. In Fig. 2(a), we plot the variance of off-diagonal matrix elements
$Var(L)$ vs. $L$ for different values of $\alpha^{2}/N$. For $\alpha^{2}/N\sim
1$, $Var(L)$ begins to fall as $L$ increases from $1$. For higher values of
$\alpha^{2}/N$, $Var(L)$ is nearly constant inside a band around the diagonal
with $|i-j|<b$, and decays for $|i-j|>b$. We define $b$ as the value of $L$
where $Var(L)$ decays to half of $Var(1)$. The bandwidth $b$ increases
linearly with $\alpha^{2}/N$, and $Var(L=1)\sim N/\alpha^{2}$ [see Figs. 2(b)
and 2(c)]. The fall of $Var(L)$ outside the band for a fixed value of
$\alpha^{2}/N<N$ is linear on the log-log scale of the figure, showing that it
has a power-law decay with $L$. (This power-law behavior will be further
analyzed shortly.) Now, looking at data for $\alpha^{2}/N\sim\mathcal{O}(N)$
in the same figure, we see that the band has covered the entire matrix,
$b\simeq N/2$. Therefore, the matrix $p$ becomes fully random with $Var(L)$
fluctuating around the theoretical value ($\simeq N/12$) for all $L$.

In Fig. 2(d), we plot the scaled variables $y=Var(L)/Var(1)$ against
$x=(L-1)/b$ on a log-log scale. This yields an excellent data collapse. The
solid line denotes the best-fit to the collapsed data, and has the form

$y=\frac{1}{1+x^{m}},$ (14)

with exponent $m\simeq 1.35$. At large $x$, Eq. (14) has a power-law decay:
$y\sim x^{-m}$. This peculiar behavior of $Var(L)$ shows that the TRI-breaking
perturbation, when represented in $\mathcal{U}(\lambda=0)$-diagonal basis,
belongs to a more general class of random matrices, known as power-law random
band matrices (PRBM) Mirlin  _et al._ (1996); Mirlin (2000). In PRBM, the
variance of off-diagonal elements decays as a power-law with distance from the
diagonal.

Table 1: The variance $v^{2}$ of significant matrix elements around the diagonal for different values of $\alpha^{2}/N$ with $N$ = 2001. $\alpha^{2}/N$ |  $v^{2}$
---|---
5 | 8227.92
10 | 3822.59
25 | 1515.06
50 | 778.40
100 | 428.63

We next discuss the main theme of this paper. Fig. 2(b) shows that $Var(L)$
near the diagonal (at $L=1$) is higher for smaller values of $\alpha^{2}/N$.
The reason behind this increase is the enhancement of overlaps between
neighboring eigenstates of $\mathcal{U}(\lambda=0)$ with a decrease in
$\alpha^{2}/N$. The matrix elements $p_{ij}$ fluctuate strongly when $p$ is
represented in such highly correlated neighboring eigenstates. Therefore, the
variance of $p_{ij}$ has higher values near the diagonal. One can see from Eq.
(3) that, even for a small perturbation strength $\lambda$, the transition
parameter $\Lambda$ becomes large due to an increase in the variance.
Therefore, the TRI $\rightarrow$ TRNI transition in terms of $\Lambda$ is
enhanced for small values of $\alpha^{2}/N$. (We remind the reader that the
TRI $\rightarrow$ TRNI transition only coincides with the COE $\rightarrow$
CUE transition in the limit $\alpha^{2}/N\sim N$.)

In Table 1, we present values of variance $v^{2}$ for different
$\alpha^{2}/N$. We compute $v^{2}$ as the variance of all elements of matrix
$p$ up to the bandwidth $b$. These values will be used to obtain the
transition parameter $\Lambda$ in the following discussion.

##  V  Transition in Spectral Fluctuations

When TRI-breaking takes place in the QKR, the resistance of levels to
crossings in the spectrum becomes greater, i.e., the rigidity of the spectrum
is increased. In the previous section, we have predicted an enhancement in the
TRI $\rightarrow$ TRNI transition for small $\alpha^{2}/N$, which can be
observed via a rapid increase in rigidity of the underlying spectrum. For
this, we investigate the transition curves of the number variance
$\Sigma^{2}(r)$, defined in Eq. (12), as a function of the transition
parameter $\Lambda$ for different values of $\alpha^{2}/N$.

The TRI $\rightarrow$ TRNI transition curves for two-point fluctuations in
Gaussian ensembles were obtained in Mehta and Pandey (1983); Pandey and Mehta
(1983). These are also valid for the COE $\rightarrow$ CUE transition in the
large-$N$ limit Pandey and Shukla (1991):

$\Sigma^{2}(r,\Lambda)=\Sigma^{2}(r,\infty)-2\int_{0}^{r}\mathrm{d}s~{}(r-s)\left[~{}Y_{2}(s,\Lambda)-Y_{2}(s,\infty)~{}\right].$
(15)

Here, $Y_{2}(s,\Lambda)$ is the two-level cluster function for the transition
ensemble:

$Y_{2}(s,\Lambda)=Y_{2}(s,\infty)-\left[\int_{0}^{1}\mathrm{d}x~{}\mathrm{e}^{2\pi^{2}\Lambda
x^{2}}x~{}\sin(\pi
xs)\int_{1}^{\infty}\mathrm{d}y~{}\mathrm{e}^{-2\pi^{2}\Lambda
y^{2}}~{}\frac{\sin(\pi ys)}{y}\right].$ (16)

In Fig. 3, we present numerical results for the transition curves of
$\Sigma^{2}(r,\Lambda)$ vs. $\Lambda$ for various values of $\alpha^{2}/N$ in
interval sizes $r=1$ and $2$. The transition parameter $\Lambda$ for different
values of $\lambda$ is calculated from Eq. (3). The analytical result for the
COE $\rightarrow$ CUE transition is obtained by numerically integrating Eqs.
(15)-(16).

Let us first consider Fig. 3(a) for $r=1$. When $\alpha^{2}/N$ is small
($=5$), one can see a smooth transition in terms of $\Lambda$, but lying above
the COE $\rightarrow$ CUE transition curve. The limit $\Lambda=0$ corresponds
to the TRI-preserved case, and $\Lambda=1$ to the completion of the TRI
$\rightarrow$ TRNI transition. For large $\alpha^{2}/N\gtrsim 25$, the
transition curves of $\Sigma^{2}(1,\Lambda)$ follow the COE $\rightarrow$ CUE
transition. In Fig. 3(b), we examine the transition curves for
$\Sigma^{2}(2,\Lambda)$. We observe again that the transition in terms of
$\Lambda$ is smooth. However, due to the larger interval size, the curves for
small $\alpha^{2}/N$ are now more separated from the COE $\rightarrow$ CUE
curve. For large $\alpha^{2}/N$, the transition curves again coincide with the
COE $\rightarrow$ CUE curve. We have examined the transition curves of
$\Sigma^{2}(r,\Lambda)$ in even larger intervals of size $r=5,10$ (not shown
here). For a fixed $r$, the TRI $\rightarrow$ TRNI transition reduces to the
COE $\rightarrow$ CUE transition in the limit $\alpha^{2}/N\gg r$. However, at
a fixed value of $\alpha^{2}/N<N$, the transition starts deviating
significantly from the COE $\rightarrow$ CUE transition for
$r\sim\alpha^{2}/N$.

The contribution from the variance of significant perturbation matrix elements
in $\Lambda$ is dependent on the parameter $\alpha^{2}/N$ and is large for
small $\alpha^{2}/N$ (see Table 1). Therefore, the required perturbation
strength $\lambda$ for a fixed value of $\Lambda$ is smaller for smaller
$\alpha^{2}/N$ — see Eq. (3).

Table 2: The perturbation strength $\lambda$ where $\Sigma^{2}(r,0)$ decays to $[\Sigma^{2}(r,0)+\Sigma^{2}(r,1)]/2$ for different values of $\alpha^{2}/N$ with $N$ = 2001. $\alpha^{2}/N$ | $\lambda\times 10^{5}$
---|---
$r$ = 1 |  $r$ = 2
5 | 0.46 | 0.69
10 | 0.83 | 1.04
25 | 1.32 | 1.65
50 | 1.85 | 2.32
$\mathcal{O}(N)$ | 3.84 | 5.16

Let us obtain a more quantitative understanding with the present scenario. We
obtain the values of $\lambda$ where the number variance has decayed from
$\Sigma^{2}(r,0)$ to $[\Sigma^{2}(r,0)+\Sigma^{2}(r,1)]/2$ for the data sets
in Figs. 3(a)-(b). In Table 2, we list these values of $\lambda$ for different
$\alpha^{2}/N$. For both $r=1$ and $2$, the value of $\lambda$ needed to drive
transitions from $\Lambda=0$ to $\Lambda=1$ becomes less with decrease in
$\alpha^{2}/N$. Therefore, we conclude that the TRI $\rightarrow$ TRNI
transition in spectral fluctuations is enhanced when $\alpha^{2}/N$ is
decreased below $N$.

##  VI  Eigenvector Fluctuations

We next study the eigenvectors of the evolution operator. Since the transition
in eigenvector fluctuations is also governed by $\Lambda$ French  _et al._
(1985, 1988a), an enhancement in the TRI $\rightarrow$ TRNI transition should
also show up in eigenvector fluctuations. We will work in a basis where the
eigenvectors of $\mathcal{U}(\lambda=0)$ are real. Then, the speed of the TRI
$\rightarrow$ TRNI transition describes the rate at which an eigenvector
becomes complex with increase in $\lambda$ from $\lambda=0$.

Let us consider the one-point eigenvector fluctuation measure, i.e., the
probability distribution of the absolute square of an eigenvector component
$y=|c_{n}|^{2}$, with $c_{n}$ denoting the $n\textsuperscript{th}$ component
of an eigenvector. In RMT, the probability distribution function (normalized
to unit mean) of $y$ has the form Brody  _et al._ (1981); Rosenzweig  _et al._
(1968)

$P(y)=\frac{\left(\nu/2\right)^{1/2}}{\Gamma\left(\nu/2\right)}~{}y^{\nu/2-1}\rm{e}^{-\nu
y/2},$ (17)

where $\nu=1$ for COE and 2 for CUE, respectively. The function $P(y)$ changes
from a $\chi^{2}_{1}$-distribution with large variance ($\sigma^{2}=2$) to a
$\chi^{2}_{2}$-distribution with smaller variance ($\sigma^{2}=1$) during the
COE $\rightarrow$ CUE transition. The variance $\sigma^{2}$ in terms of
$\Lambda$ continuously changes from $2$ to $1$ in this transition French  _et
al._ (1988a). We stress that RMT results for eigenvector statistics in the
crossover regime are still not fully available. There are a few preliminary
works, but these do not precisely capture the intermediate regime (see Sommers
and Iida (1994); Chung  _et al._ (2000) and references therein).

In this section, we numerically investigate the eigenvector fluctuations for
different values of $\alpha^{2}/N$, where $N=2001$. Note that the eigenvector
statistics is basis-dependent. Therefore, we chose a generic basis, i.e., the
position basis. We first calculate the eigenvectors of the evolution operator
matrix $\mathcal{U}$ in Eq. (11) using the standard LAPACK routines Anderson
_et al._ (1999). Then, we find the absolute squares $|c_{n}|^{2}$ from all the
eigenvectors of a single matrix and normalize them to mean unity. For
eigenvector statistics, a single realization of the $\mathcal{U}$-matrix
(yielding $2001\times 2001$ eigenvector components) is sufficient to achieve
good accuracy in data. To compare the probability distributions $P(y)$ for
different $\alpha^{2}/N$, we plot $P(\log_{10}(y))$ vs. $\log_{10}(y)$. The
RMT results for the COE $\rightarrow$ CUE transition are obtained numerically
by using an interpolating ensemble of large Gaussian random matrices.

Table 3: (i) Values of $\Lambda$ at $\lambda=1.719\times 10^{-5}$ for different values of $\alpha^{2}/N$ with $N=2001$. (ii) Corresponding values of $\sigma^{2}(\lambda=0)-\sigma^{2}(\lambda=1.719\times 10^{-5})$. $\alpha^{2}/N$ |  $\Lambda$ |  $\sigma^{2}(\lambda=0)-\sigma^{2}(\lambda=1.719\times 10^{-5})$
---|---|---
5 | 0.25 | 0.92
10 | 0.11 | 0.80
25 | 0.05 | 0.52
50 | 0.02 | 0.35
$\mathcal{O}(N)$ | 0.005 | 0.12

In Fig. 4(a), we plot $P(\log_{10}(y))$ against $\log_{10}(y)$ for different
$\alpha^{2}/N$ at $\lambda=0$. The data sets for small $\alpha^{2}/N$
($=5,10$) show a small departure from the $\chi^{2}_{1}$-distribution around
the peak – this region is amplified in the inset. On the other hand, for
$\alpha^{2}/N\geq 25$, the data shows an excellent agreement with the
$\chi^{2}_{1}$-distribution (COE). In Fig. 4(b), we show data for different
$\alpha^{2}/N$ at $\lambda=0.9$. In this case, all data sets are in good
agreement with the $\chi^{2}_{2}$-distribution (CUE). This shows that the
eigenvectors become complex random states (irrespective of $\alpha^{2}/N$)
when TRI is broken by a strong magnetic field. We stress that the one-point
eigenvector fluctuations capture only the symmetries present in the evolution
operator.

In Fig. 4(c), we show the corresponding data for a small perturbation strength
$\lambda=1.719\times 10^{-5}$. In this case, we expect an intermediate
behavior between $\chi^{2}_{1}$ and $\chi^{2}_{2}$, depending on the value of
$\Lambda$. Fig. 4(c) shows that the data sets for smaller values of
$\alpha^{2}/N$ approach closer to the limiting $\chi^{2}_{2}$-distribution.
Therefore, the eigenvectors at large $\alpha^{2}/N$ require larger
perturbations to reach the $\chi^{2}_{2}$ or CUE-limit. For the data in Fig.
4(c), we present in Table 3 the following: (i) the values of the transition
parameter $\Lambda$ from Eq. (3), and (ii) the difference between the
variances of the distributions at $\lambda=0$ and $\lambda=1.719\times
10^{-5}$. For $\alpha^{2}/N=5$, $\Lambda\simeq 0.25$ with a large decrease in
$\sigma^{2}\simeq 0.92$. For $\alpha^{2}/N\sim\mathcal{O}(N)$, $\Lambda\simeq
0.005$ and the change in $\sigma^{2}$ is relatively minor $\simeq 0.12$. This
confirms that the eigenvector fluctuations for small $\alpha^{2}/N$ are more
sensitive to $\lambda$. The reason behind this behavior is again the increased
variance of perturbation matrix elements, which causes the strong admixing of
imaginary amplitude in unperturbed (real) eigenvectors.

In Fig. 5, we examine the transition in terms of $\Lambda$. We plot
$\sigma^{2}$ vs. $\Lambda$ for different values of $\alpha^{2}/N$. The solid
curve denotes the COE $\rightarrow$ CUE transition. The transition curves for
$\alpha^{2}/N<N$ show a departure from the COE $\rightarrow$ CUE curve. (For
$\alpha^{2}/N\sim\mathcal{O}(N)$, the transition curve is in good agreement
with the COE $\rightarrow$ CUE curve.) Further, due to the increment in
$v^{2}$, the speed of the transition is enhanced for small $\alpha^{2}/N$. In
the inset of the same figure, we plot the distribution $P(\log_{10}(y))$ vs.
$\log_{10}(y)$ for $\Lambda=0.5$. For $\alpha^{2}/N\sim\mathcal{O}(N)$, the
distribution is in very good agreement with the COE $\rightarrow$ CUE
transition. As $\alpha^{2}/N$ is decreased below $N$, the distribution shows a
deviation from the COE $\rightarrow$ CUE curve and becomes broader around the
mean. This fact can be appreciated in the main frame, where
$\sigma^{2}(\Lambda)$ is seen to be larger for smaller values of
$\alpha^{2}/N$.

##  VII  Summary and Discussion

Let us conclude this paper with a summary and discussion of our results. We
have investigated the transition from time-reversal invariance (TRI) to time-
reversal non-invariance (TRNI) in the spectral and eigenvector fluctuations of
the quantum kicked rotor (QKR) using Izrailev’s finite-dimensional model. The
local fluctuation properties of the QKR depend on the parameter
$\alpha^{2}/N$, where $\alpha$ is the chaos parameter and $N$ is the
dimensionality of the evolution operator matrix. The TRI in the QKR is broken
immediately by the application of a small magnetic field. However, there is a
continuous crossover from TRI $\rightarrow$ TRNI in the spectral and
eigenvector fluctuations. This is governed by a local parameter $\Lambda$,
which depends on the magnetic field and the near-diagonal matrix elements of
the TRI-breaking perturbation. For $\alpha^{2}/N\gtrsim N$, the TRI
$\rightarrow$ TRNI transition in the QKR is described well by the COE
$\rightarrow$ CUE transition of random matrix theory (RMT). This scenario
applies if the TRI-breaking perturbation, when represented in the diagonal
basis of the TRI-preserving part, becomes a member of Gaussian RMT ensembles.

However, for $\alpha^{2}/N<N$, the TRI-breaking perturbation no longer
satisfies this. Rather, it has the form of a power-law random band matrix
(PRBM), whose bandwidth grows linearly with $\alpha^{2}/N$. The power-law fall
of the variance of perturbation matrix elements outside the band has a decay
exponent $m\simeq 1.35$, which is the same for all values of $\alpha^{2}/N<N$.
(This value of $m$ is reminiscent of so-called sub-critical PRBM models
($1<m<2$) Mirlin  _et al._ (1996); Méndez-Bermúdez  _et al._ (2010).)

In addition to the PRBM structure, $Var(1)$ increases linearly with
$N/\alpha^{2}$. The increased variance for small $\alpha^{2}/N$ enhances the
sensitivity of $\Lambda$ to the perturbation strength, resulting in a faster
TRI $\rightarrow$ TRNI transition as a function of the magnetic field.

For spectral fluctuations, we investigated the transition curves of rigidity
measure [i.e., number variance $\Sigma^{2}(r,\Lambda)$, where $r$ is the size
of the spectral interval] by varying $\Lambda$ at fixed $r$. For small
$\alpha^{2}/N$, the transition curve deviates from the COE $\rightarrow$ CUE
transition. Further, the speed of this transition is enhanced for smaller
values of $\alpha^{2}/N$. As we increase $\alpha^{2}/N$, the COE $\rightarrow$
CUE transition is approached. The eigenvector fluctuations are studied via the
probability distributions $P(y)$ of the absolute square of the eigenvector
component $y$. For a small increment in $\lambda$, the variance $\sigma^{2}$
of $P(y)$ changes substantially for small $\alpha^{2}/N$. We also analyzed the
transition of variance $\sigma^{2}$ with respect to the parameter $\Lambda$.
For $\alpha^{2}/N<N$, the transition curve of $\sigma^{2}$ deviates from the
COE $\rightarrow$ CUE curve. For $\alpha^{2}/N\sim\mathcal{O}(N)$, the TRI
$\rightarrow$ TRNI transition in eigenvectors is consistent with the COE
$\rightarrow$ CUE transition. Again, the speed of the transition in the
eigenvector fluctuations is higher for smaller values of $\alpha^{2}/N$ below
$N$.

A similar enhancement in the TRI $\rightarrow$ TRNI transition was predicted
by French et al. French  _et al._ (1987, 1988b) in the context of a nuclear
Hamiltonian. These authors modeled the TRI-preserving and TRI-breaking parts
of the Hamiltonian as $m$-particle embedded-GOE ensembles with $k$-body
interactions. In the limit $m\gg k$, the near-diagonal variance of the
perturbation matrix was seen to increase. This resulted in an enhancement of
the transition. However, a semi-classical treatment of this enhanced
transition remains an open problem.

###### Acknowledgements.

We would like to thank Avanish Kumar for discussions.

##  Appendix A Time-evolution Operator $\mathcal{U}$ in Infinite-dimensional
Position Basis

Consider the form of the operator $\mathcal{U}$ in Eq. (8). The matrix
elements of $\mathcal{U}$ in position basis
$\\{|q\rangle:\theta|q\rangle=\theta_{q}|q\rangle\\}$ are given by

$\mathcal{U}_{qq^{\prime}}=\langle q|B^{1/2}GB^{1/2}|q^{\prime}\rangle.$ (18)

The position basis for the QKR is periodic $(|q+2\pi\rangle=|q\rangle)$. Then,
the orthonormality and completeness relations in the position basis
$\\{|q\rangle\\}$ can be defined as Albert  _et al._ (2017); Raynal  _et al._
(2010)

$\langle
q|q^{\prime}\rangle=2\pi~{}\delta^{(2\pi)}(q-q^{\prime}),\quad\mbox{and}\quad\int_{(2\pi)}\frac{dq}{2\pi}~{}|q\rangle\langle
q|=\mathbb{I},$ (19)

where the superscript $(2\pi)$ on the $\delta$-function indicates
$2\pi$-periodicity. The integration for variable $q$ is over a $2\pi$
interval.

Due to the periodicity of the position basis, the eigenvalues of the momentum
operator are discrete, i.e., $p|l\rangle=l\hbar|l\rangle$; $l=0,\pm 1,\pm
2,\ldots$. The orthonormality and completeness relations in the momentum basis
are expressed as

$\langle
l|l^{\prime}\rangle=\delta_{l,l^{\prime}},\quad\mbox{and}\quad\sum_{l=-\infty}^{l=\infty}|l\rangle\langle
l|=\mathbb{I}.$ (20)

Further, the inner products between the basis vectors of the discrete momentum
basis and the continuous position basis are defined as

$\langle q|l\rangle={\rm e}^{il\theta_{q}},\quad\langle l|q\rangle={\rm
e}^{-il\theta_{q}}.$ (21)

By using the completeness relations from Eqs. (19)–(20), the matrix element
$\mathcal{U}_{qq^{\prime}}$ in Eq. (18) can be written as

$\mathcal{U}_{qq^{\prime}}=\frac{1}{(2\pi)^{2}}\sum_{l,l^{\prime}=-\infty}^{\infty}\int_{(2\pi)}\int_{(2\pi)}dr~{}dr^{\prime}~{}\langle
q|B^{1/2}|r\rangle~{}\langle r|l\rangle~{}\langle
l|G|l^{\prime}\rangle~{}\langle l^{\prime}|r^{\prime}\rangle~{}\langle
r^{\prime}|B^{1/2}|q^{\prime}\rangle.$ (22)

Here, the continuous variables $q$, $q^{\prime}$, $r$ and $r^{\prime}$ label
the position basis, and the integers $l$ and $l^{\prime}$ label the momentum
basis.

By using Eq. (7) for $G$ in Eq. (22), and simplifying the expression using
Eqs. (19)–(21), we obtain

$\mathcal{U}_{qq^{\prime}}=B^{1/2}(\theta_{q})B^{1/2}(\theta_{q}^{\prime})\sum_{l=-\infty}^{l=\infty}{\rm
e}^{-iT(l\hbar-\lambda)^{2}/(2\hbar)}~{}{\rm
e}^{il(\theta_{q}-\theta_{q}^{\prime})}.$ (23)

Further, recalling the definition of $B$ from Eq. (6), the matrix elements
$\mathcal{U}_{qq^{\prime}}$ take the form given in Eq. (9) of the main text.

## References

  * Casati  _et al._ (1979) G. Casati, B. V. Chirikov, F. M. Izrailev, and J. Ford, in  _Stochastic Behavior in Classical and Quantum Hamiltonian Systems_ , edited by G. Casati and J. Ford (Springer Berlin Heidelberg, Berlin, Heidelberg, 1979) pp. 334–352.
  * Bohigas  _et al._ (1984) O. Bohigas, M. J. Giannoni, and C. Schmit, Phys. Rev. Lett. 52, 1 (1984).
  * Brody  _et al._ (1981) T. A. Brody, J. Flores, J. B. French, P. A. Mello, A. Pandey, and S. S. M. Wong, Rev. Mod. Phys. 53, 385 (1981).
  * Mehta (2004) M. L. Mehta,  _Random Matrices_, Pure and Applied Mathematics, Vol. 142 (Elsevier, 2004).
  * Haake (2010) F. Haake,  _Quantum Signatures of Chaos_ , Springer Series in Synergetics, Vol. 54 (Springer-Verlag Berlin Heidelberg, 2010).
  * Beenakker (1997) C. W. J. Beenakker, Rev. Mod. Phys. 69, 731 (1997).
  * Guhr  _et al._ (1998) T. Guhr, A. Muller–Groeling, and H. A. Weidenmuller, Physics Reports 299, 189 (1998).
  * Izrailev (1986) F. M. Izrailev, Phys. Rev. Lett. 56, 541 (1986).
  * Izrailev (1990) F. M. Izrailev, Physics Reports 196, 299 (1990).
  * Pandey  _et al._ (1993) A. Pandey, R. Ramaswamy, and P. Shukla, Pramana 41, L75–L81 (1993).
  * Lenz and Haake (1991) G. Lenz and F. Haake, Phys. Rev. Lett. 67, 1 (1991).
  * Dupuis and Montambaux (1991) N. Dupuis and G. Montambaux, Phys. Rev. B 43, 14390 (1991).
  * Bohigas  _et al._ (1995) O. Bohigas, M. J. Giannoni, A. M. O. de Almeida, and C. Schmit, Nonlinearity 8, 203 (1995).
  * Pluhar  _et al._ (1995) Z. Pluhar, H. Weidenmuller, J. Zuk, C. Lewenkopf, and F. Wegner, Annals of Physics 243, 1 (1995).
  * Berry and Robnik (1986) M. V. Berry and M. Robnik, Journal of Physics A: Mathematical and General 19, 649 (1986).
  * Rosenzweig  _et al._ (1968) N. Rosenzweig, J. Monahan, and M. Mehta, Nuclear Physics A 109, 437 (1968).
  * French  _et al._ (1985) J. B. French, V. K. B. Kota, A. Pandey, and S. Tomsovic, Phys. Rev. Lett. 54, 2313 (1985).
  * Fal’ko and Efetov (1994) V. I. Fal’ko and K. Efetov, Phys. Rev. B 50, 11267 (1994).
  * Chung  _et al._ (2000) S.-H. Chung, A. Gokirmak, D.-H. Wu, J. S. A. Bridgewater, E. Ott, T. M. Antonsen, and S. M. Anlage, Phys. Rev. Lett. 85, 2482 (2000).
  * French  _et al._ (1988a) J. French, V. Kota, A. Pandey, and S. Tomsovic, Annals of Physics 181, 198 (1988a).
  * Mehta and Pandey (1983) M. L. Mehta and A. Pandey, Journal of Physics A: Mathematical and General 16, 2655 (1983).
  * Pandey and Mehta (1983) A. Pandey and M. L. Mehta, Communications in Mathematical Physics 87, 449–468 (1983).
  * Dyson (1962) F. J. Dyson, Journal of Mathematical Physics 3, 1191 (1962).
  * Pandey and Shukla (1991) A. Pandey and P. Shukla, Journal of Physics A: Mathematical and General 24, 3907 (1991).
  * Pandey (2004) A. Pandey, Phase Transitions 77, 835 (2004).
  * Pandey  _et al._ (2017) A. Pandey, A. Kumar, and S. Puri, Phys. Rev. E 96, 052211 (2017).
  * Pandey  _et al._ (2020) A. Pandey, A. Kumar, and S. Puri, Phys. Rev. E 101, 022217 (2020).
  * Kumar  _et al._ (2020) A. Kumar, A. Pandey, and S. Puri, Phys. Rev. E 101, 022218 (2020).
  * Casati  _et al._ (1990) G. Casati, I. Guarneri, F. Izrailev, and R. Scharf, Phys. Rev. Lett. 64, 5 (1990).
  * Shukla and Pandey (1997) P. Shukla and A. Pandey, Nonlinearity 10, 979 (1997).
  * Shukla (1996) P. Shukla, Phys. Rev. E 53, 1362 (1996).
  * Mirlin  _et al._ (1996) A. D. Mirlin, Y. V. Fyodorov, F.-M. Dittes, J. Quezada, and T. H. Seligman, Phys. Rev. E 54, 3221 (1996).
  * Mirlin (2000) A. D. Mirlin, Physics Reports 326, 259 (2000).
  * Manos and Robnik (2013) T. Manos and M. Robnik, Phys. Rev. E 87, 062905 (2013).
  * Ott (1993) E. Ott,  _Chaos in Dynamical Systems_ (Cambridge University Press, 1993).
  * Manos and Robnik (2015) T. Manos and M. Robnik, Phys. Rev. E 91, 042904 (2015).
  * Fishman  _et al._ (1982) S. Fishman, D. R. Grempel, and R. E. Prange, Phys. Rev. Lett. 49, 509 (1982).
  * Anderson (1958) P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Sommers and Iida (1994) H.-J. Sommers and S. Iida, Phys. Rev. E 49, R2513 (1994).
  * Anderson  _et al._ (1999) E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_ , 3rd ed. (Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999).
  * Méndez-Bermúdez  _et al._ (2010) J. A. Méndez-Bermúdez, V. A. Gopar, and I. Varga, Phys. Rev. B 82, 125106 (2010).
  * French  _et al._ (1987) J. B. French, V. K. B. Kota, A. Pandey, and S. Tomsovic, Phys. Rev. Lett. 58, 2400 (1987).
  * French  _et al._ (1988b) J. French, V. Kota, A. Pandey, and S. Tomsovic, Annals of Physics 181, 235 (1988b).
  * Albert  _et al._ (2017) V. V. Albert, S. Pascazio, and M. H. Devoret, J. Phys. A: Math. Theor. 50, 504002 (2017).
  * Raynal  _et al._ (2010) P. Raynal, A. Kalev, J. Suzuki, and B.-G. Englert, Phys. Rev. A 81, 052327 (2010).

Figure 1: Number variance $\Sigma^{2}(r)$ vs. $r$ for various values of
$\alpha^{2}/N$ (with $N=2001$) for (a) $\lambda=0$, and (b) $\lambda=0.9$. The
dashed and solid curves denote the Poisson and RMT results, respectively.
Figure 2: (a) Plot of $Var(L)$ vs. $L$ on a log-log scale for different values
of $\alpha^{2}/N$. The solid line denotes the semi-classical result:
$Var(L)\simeq N/12$ with $N=2001$. (b) Plot of $Var(1)$ vs. $\alpha^{2}/N$ on
a log-log scale. The solid line denotes the slope $-1$. (c) Plot of $b$ vs.
$\alpha^{2}/N$ on a log-log scale. The solid line denotes the slope $+1$. (d)
Plot of $Var(L)/Var(1)$ vs. $(L-1)/b$ on a log-log scale for the data sets in
(a). The dashed curve is the best-fit to Eq. (14) with $m\simeq 1.35$. Figure
3: Transition curves of number variance, $\Sigma^{2}(r,\Lambda)$ vs.
$\Lambda$, for different values of $\alpha^{2}/N$ when (a) $r=1$, and (b)
$r=2$. The parameter $\Lambda$ is estimated using Eq. (3). The solid curves
are theoretical results for the COE $\rightarrow$ CUE transition. The dashed
lines indicate the COE and CUE limits. Figure 4: Probability distribution
$P(\log_{10}(y))$ vs. $\log_{10}(y)$ for various values of $\alpha^{2}/N$. The
values of $\lambda$ are as follows: (a) $\lambda=0$. The solid curve denotes
the $\chi^{2}_{1}$-distribution. The inset magnifies the peak region. (b)
$\lambda=0.9$. The solid curve denotes the $\chi^{2}_{2}$-distribution. (c)
$\lambda=1.719\times 10^{-5}$. The solid and dot-dashed curves denote the
$\chi^{2}_{1}$ and $\chi^{2}_{2}$-distributions, respectively. Figure 5:
Transition curves of the variance, $\sigma^{2}(\Lambda)$ vs. $\Lambda$, for
various values of $\alpha^{2}/N$. The solid curve denotes the COE
$\rightarrow$ CUE transition. The dashed lines indicate the COE and CUE
limits. The inset plots $P(\log_{10}(y))$ vs. $\log_{10}(y)$ for the same
values of $\alpha^{2}/N$, and a fixed value of $\Lambda=0.5$. The solid curve
in the inset is the distribution at $\Lambda=0.5$ for the COE $\rightarrow$
CUE transition.