# Dynamical clustering of U.S. states reveals four distinct infection patterns
that predict SARS-CoV-2 pandemic behavior

Joseph L. Natale  Halıcıoğlu Data Science Institute, University of California,
San Diego, La Jolla, CA 90293, USA  Varun Viswanath  Department of Electrical
and Computer Engineering, University of California, San Diego, La Jolla, CA
90293, USA  Oscar Trujillo Acevedo  División de Ingeniería Bioquímica,
Instituto Tecnológico Superior of Atlixco, C. Heliotropo No.1201, Unidad 8
Norte Nueva Xalpatlaco, Vista Hermosa, 74210 Atlixco, Pue.  Sophia Pérez
Giottonini  Colegio de Bachilleres del Estado de Sonora, plantel Villa de
Seris, Hermosillo, SON, México.  Sandy Ihuiyan Romero Hernández  Instituto
Tecnológico Superior de Atlixco, C. Heliotropo No.1201, Unidad 8 Norte Nueva
Xalpatlaco, Vista Hermosa, 74210 Atlixco, Pue., México  Diana G. Cruz Millán
Department of Health Sciences, Universidad de Sonora, campus Cajeme, Cd.
Obregón, Son., México  A. Montserrat Palacios-Puga  Department of Biochemical
and Chemical Engineering, Instituto Tecnológico de Tijuana, Tecnológico
Nacional de México, Tijuana, B.C, México  Ammar Mandvi  Department of
Emergency Medicine, Department of Family Medicine, University of California,
San Diego, CA, USA  Brian M. Khan  Department of Emergency Medicine,
Department of Internal Medicine, University of California, San Diego, CA, USA
Martin Lilik  Department of Software Engineering, Park Media, Semarang,
Indonesia  Jay Park  Department of Emergency Medicine, University of
California, San Diego, La Jolla, CA 92037  Benjamin L. Smarr  Department of
Bioengineering and Halıcıoğlu Data Science Institute, University of
California, San Diego, La Jolla, CA 90293, USA

(November 19, 2021)

###### Abstract

The SARS-CoV-2 pandemic has so far unfolded diversely across the fifty United
States of America, reflected both in different time progressions of infection
“waves” and in magnitudes of local infection rates. Despite a marked diversity
of presentations, most U.S. states experienced their single greatest surge in
daily new cases during the transition from Fall 2020 to Winter 2021. Popular
media also cite additional similarities between states – often despite
disparities in governmental policies, reported mask-wearing compliance rates,
and vaccination percentages. Here, we identify a set of robust, low-
dimensional clusters that 1) summarize the timings and relative heights of
four historical COVID-19 “wave opportunities” accessible to all 50 U.S.
states, 2) correlate with geographical and intervention patterns associated
with those groups of states they encompass, and 3) predict aspects of the
“fifth wave” of new infections in the late Summer of 2021. In particular, we
argue that clustering elucidates a negative relationship between vaccination
rates and subsequent case-load variabilities within state groups. We advance
the hypothesis that vaccination acts as a “seat belt,” in effect constraining
the likely _range_ of new-case upticks, even in the context of the Summer
2021, variant-driven surge.

##  I Introduction

The SARS-CoV-2 pandemic has claimed millions of lives globally – and directly
affected hundreds of millions more – doing so in a highly unequal manner Van
Dorn et al. (2020); Chaudhry et al. (2020). Early foci for COVID-19, such as
northern Italy, coincided with populations with high proportions of elderly,
susceptible people Ciminelli and Garcia-Mandicó (2020); the availability of
hospital beds, treatments, and clinical personnel are again to be a major
problem in many geographical areas Weissman et al. (2020); Sen-Crowe et al.
(2021); and new variants are emerging at different rates in different locales
Alfaro et al. (2021); Chadha et al. (2021).

Internal to the United States of America (at the time of writing, that nation
with the highest recorded numbers of both “total” COVID-19 cases, and
COVID-19-related deaths) disparities also permeate local infection rates
Loomba et al. (2021); Hamman (2021). This has prompted a number of explicit
comparisons, at the state level Xian et al. (2020); Solís et al. (2021). In
principle, such collations should be of great interest to public health
decision-makers: to the extent that different states can be expected to behave
in categorically different ways, information that supports early planning for
interventional responses can be gleaned by mapping the states’ diverse
progressions over time.

Publicly available data render unambiguous that all fifty U.S. states saw
their greatest, sustained incidence of new cases during the country’s last
winter season (2020 - 2021) for Disease Control et al. (2021). Furthermore,
media outlets have directed attention to other, ostensible similarities: in
the case of Florida and California, news coverage emphasized an attainment of
“similar COVID-19 case rates” Curley (2021), or even allegedly “identical
outcomes” News (2021), despite the drastically different measures taken by the
independently-operating, state-level governments at the time of reporting.
Still, in complex systems it is always possible that multiple parts or
subsystems exhibit temporarily superficial similarities without sharing their
underlying dynamics Strogatz (2018). Whether or not the U.S. state infection
histories fall into universal similarity classes with relevance to pandemic
planning – and, indeed, questions of how reliably government policy, and even
vaccination rates, determine their ultimate case statistics – remains
unanswered by these prior analyses.

In this manuscript we consider one specific approach to demonstrate that the
historical “infection trajectories” for all fifty U.S. states can be grouped
according to some reasonable measure of similarity, such that resulting groups
form a low-dimensional set (i.e., one much smaller than the number of
individual states), and also prove useful for predicting the present, or even
future, of this pandemic – not merely its past. As such groupings are
established, any relevant information supporting early planning around
responses to the states’ diverse presentations can be contextualized and made
available to public health decision-makers.

##  II Methods

###  II.1 Data Acquisition and Preparation

####  II.1.1 Infection Data

The data referenced in the subsequent paragraphs were aggregated from sources
developed by Johns Hopkins and Oxford Universities. The infection data,
reported at daily resolution, were obtained from Johns Hopkins Repository of
COVID-19 data underlying the analysis and visualizations created by the Johns
Hopkins Centers for Civic Impact for the Coronavirus Resource Center (CRC).
These included reports of _total_ (confirmed and probable) _cases_ , expressed
as cumulative counts for each U.S. state. In order to transform these
cumulative records into “daily new” infection counts for each state, we
subtracted each day’s _total cases_ record from the previous day’s value, and
then replaced their raw counts with the corresponding moving (14-rolling)
averages attained by each day considered. In so doing, we therefore neglected
the effects of deaths and recoveries on the “daily new” case numbers; we
justified this approximation by noting that a) the latter were often as small
as a full order of magnitude less than the total cases on a given day, and b)
isolated “negative” days (those for which deaths and recoveries presumably
exceeded the new, emerging case count) would be absorbed, without a strong
effect, by the aforementioned temporal averaging. Training data were amassed
between April 2020 to July 2021.

Since the different U.S. states encompass populations with markedly different
sizes and densities, but our temporal picture seeks to distill all the states’
histories into a form for which their synchrony and structure can be compared
on equal footing, we normalized the aforementioned, raw estimates of the
states’ daily new infection counts by their respective maximum values (all
observed to occur in the same time frame, Winter 2020 - 2021). This yielded a
new set of case records for every state, with possible values between 0 and 1.
Treating each as a time series, and applying 14-day moving (rolling) averages
between April 13, 2020 and July 2, 2021 – a total of 446 calendar days –
defined our U.S. “daily new” _infection trajectories_. An example trajectory,
for New York state, is shown in Fig. 1, in blue.

Daily new infection records for dates prior to April 13, 2020 are depicted for
temporal context in Section III, but were not used for developing or training
any models: in these early stages of the pandemic, we could not be confident
that the reporting standards across states were of uniform or sufficiently
high quality Wu et al. (2020). For instance, a lack of PCR testing supplies
might have led to a greater proportion of “probable,” vs. “confirmed,” cases
so early on. Some numerical artifacts resulting from performing 14-day
averages on data that begins midway through the initial surge of cases in
Spring 2020 were removed from the figures in Section III, by adjusting the
plotting limits across the respective x-axes. During the course of this work,
data were re-scraped to extend between March 6, 2020 and August 3, 2021.

####  II.1.2 Definition of five “wave opportunities”

The example infection trajectory for New York (in Fig. 1) includes periods of
“low” relative case loads, punctuated by several dramatic – and, persistent –
elevations. A central question in the present work is whether the temporal
patterns associated with these periodic surges in the daily new cases fall
into stereotypes across states. As a starting point for comparing the time
progressions realized by different states, we modify a previously established
convention, dividing the pandemic’s history between the Spring of 2020 and
Spring 2021 into consecutive “waves” Hallas et al. (2020). Here, we
reconceptualize these divisions in terms of annual, calendar seasons:

  1. 1.

First wave: an initial case accumulation, roughly between March and May
(“Spring”) of 2020

  2. 2.

Second wave: a subsequent resurgence, roughly between July and September
(“Summer”) of 2020

  3. 3.

Third wave: the aforementioned “Winter peak,” lasting from around November
2020 through March 2021, but reaching its apex between December 2020 - January
2021

  4. 4.

Fourth wave: recalcitrant surge during the Spring of 2021 that was noted in,
but had not necessarily resolved by the time of publication of, Ref. Hallas et
al. (2020)

In addition, we consider a “fifth” wave, encompassing more recent case load
rises beginning in late July 2021 (Fig. 1). This upswing was reflected in the
daily new case records for all fifty U.S. states, but the same was not
universally true of the previous four waves. To emphasize explicitly that not
all states participated equally in each of those four, we adopt the
nomenclature _wave opportunity_ to refer to the potentially-eventful spans of
time with persistent elevations. As one example, New York expressed strongly
the “first” wave opportunity, but the “second wave” seen in other parts of the
country, as defined above, was so small relative to the former that New York
effectively progressed without participation in it. To describe the climbing
infection rates in the Winter months, we refer here not to “New York’s second
wave,” but to “New York’s (degree of) participation in the third U.S. wave
opportunity” (see again Fig. 1).

Since we are excluding data before April 13, 2020, and our rolling averages
are fully-defined only after fourteen points have been observed (i.e., from
April 26, 2020 on), the first wave opportunity is not considered in our
analyses.

####  II.1.3 Ancillary Data Sources

In addition to infection data, we aggregated additional data on vaccine
administration rates and government policy “stringencies” during the pandemic.
We obtained vaccination data from repositories published by the Johns Hopkins
University Coronavirus Resource Center Repository of COVID-19 data underlying
the analysis and visualizations created by the Johns Hopkins Centers for Civic
Impact for the Coronavirus Resource Center (CRC), and curated by Our World in
Data Repository of COVID-19 data curated by Our World In Data
(https://ourworldindata.org/). We incorporated government policy data compiled
by the Blavatnik School of Government, University of Oxford Hallas et al.
(2020).

The percentage of “fully-vaccinated” individuals residing in each state – our
chosen proxy for progress toward a local “herd immunity” status – was
estimated by, first, extracting the number of people recorded as having
received _two doses of either two-series_ (Pfizer-BioNTech or Moderna), or
_one dose of the single-shot_ (J&J/Jannsen), vaccine against COVID-19
asvailable in the U.S., _at each calendar day_ between February 1, 2021 and
July 22, 2021, inclusive. This time range, distinct from (but still
overlapping with) that of our infection data, consisted of 172 days.

We aggregated tabulated counts by state, and then divided by the number of
residents in each state. Although U.S. state population levels do change over
time, we made the tacit approximation of slow change over the course of the
pandemic: dividing by the “2021” population counts given in Ref. Review (2021)
allowed us to arrive at a cumulative percentage of fully-vaccinated residents
in each state as a time series – a _vaccination-rate_ analog of our infection
trajectories. The vaccination-rate trajectory for New York is presented as a
red curve, normalized between 0 and 1, in Fig. 1.

A “Stringency Index” score for comparing government policies across states was
developed by previous authors, to summarize the “strictness of ‘lockdown
style’ closure and containment policies that primarily restrict people’s
behavior” Hallas et al. (2020). Stringency scores could range, in principle,
from a minimum value of 0 to a maximum of 100 for a given state; rescaling the
stringency time series values provided a third, government-_policy trajectory_
, normalized to the same vertical scale as the infection and vaccine
trajectories (between 0 and 1). Policy trajectories began much earlier than
daily-infection and vaccination-rate trajectories, on January 1, 2020 (before
the “start” of the pandemic in the U.S.) and ended on July 2, 2021 – a total
of 549 days. A gray curve denotes the policy trajectory in Fig. 1.

###  II.2 Correlation-Based Clustering of Infection Data

Given the above considerations, we hypothesized that i) if placed on some even
footing, the $n=50$ U.S. “ _daily new infection_ ” _trajectories_ should
collapse into a set $\mathcal{C}$ of archetypal forms with cardinality
$|\mathcal{C}|\ll n$. That is, we set out to verify whether the records for
all 50 states could be summarized by – and compressed into – a handful of
discrete _clusters_ , each associated with a particular timing and structure
(relative peak heights) for all those distinct “wave opportunities”
encompassed by the available data, between April 2020 and July 2021. We then
investigated ii) whether one such partitioning might prove useful, not only
for summarizing the states’ dynamics in the _past_ , but also for predicting
aspects of the _subsequent_ , ongoing surge, attributed largely to the
B.1.617.2 (“Delta”) variant Allen et al. (2021).

“Daily new infection” trajectory values spanning the time range between April
12th, 2020 and July 2nd, 2021 were clustered in Python 3, using the _SciPy_
Clustering package. Specifically, we performed hierarchical / agglomerative
clustering via the named “ _linkage_ ” function in the _hierarchy_ module of
_SciPy_ (v. 1.4.1). A total of 446 data points per state, each representing a
trajectory value for 1 calendar date, served as primary input to this
function. We used a “complete” (i.e., Farthest Point, or Voor Hees, Algorithm)
linkage method to compute cluster distances. In what follows, this distance
was measured by the (Pearson) _correlation_ (see community (2021)) between
respective pairs of time series.

Using the named “ _fcluster_ ” function of the same _SciPy_ module, we
obtained the set of “flat” clusters that would result from “slicing” the
discovered hierarchy at a specific _threshold_ value of our correlation-based
distance metric. This same process could be repeated to cluster other types of
data – e.g., percent-fully-vaccinated “trajectories”.

###  II.3 Principal Components Analysis (PCA)

####  II.3.1 Discovering “important” date ranges that set states apart

Principal Components Analysis was applied directly to the infection
trajectories, treating them as 50 observations of one dynamical system, and
the 446 daily new case values as their basis “features” (those to be linearly
combined). This was accomplished using the _decomposition_ module of the
_scikit-learn_ package (v. 0.22.2.post1) in Python 3. Using the named function
“ _pca_ ,” we extracted the first 50 principal components. Since these
principal components are essentially linear combinations of the input
features, each represents a weighted sum of contributions from the rolling,
“daily new infections” on different _dates_ ; therefore, we used their feature
importance rankings to reveal the subsets of _dates_ during the pandemic
associated with the greatest degrees of variation across all 50 U.S. states.

Feature importance, in the context of one given PCA-derived eigenvector (or,
“principal component”), was scored according to absolute values: the vector
components with highest (relative) absolute values were ranked as the most
important. In order to arrive at a more all-encompassing set of dates that
took into account the importances from the first _several_ principal
components $k\in\\{1,\ldots k_{0}\\}$, we also took a step back from the usual
abstraction of principal components as (ideally, orthogonal) eigenvectors, and
treated each instead as a set of general ‘importance-amplitude” vectors,
having components for each of the 446 dates. Summing these importance-
amplitude vectors, component-wise, resulted in a new curve that we used to
represent “overall,” daily importances.

In order to select an appropriate value for $k_{0}$, we examined the
cumulative percentage of “explained variance” associated with the first
several principal components $k\in\\{1,2,\dots 10\\}$; these values are
reported in Table I, along with the “errors” corresponding to reconstructions
of the 50 original observations – our trajectories – that were built using all
principal components from 1 to $k_{0}$, inclusive. The “error” at each time
point was measured via the mean absolute deviation between the original and
reconstructed trajectory value, averaged over states.

The cutoff for choosing $k_{0}$ in terms of the cumulatively explained
variance was set at a value of $90\%$. In the context of the “overall” scheme
above, dates were considered “important” if the absolute value of a given
cumulative importance curve exceeded the mean value of that curve. An upper
limit of $k=50$ (the minimum between the number of observations, 50, and the
number of features, 466) principal components was permitted by the “ _pca_ ”
function.

The results for the most important date ranges inferred via PCA are depicted
in Section IVc, contextualized by their superposition over all $n=50$ original
state trajectories, with a common axis. There is no intrinsic guarantee that
PCA should discover date ranges that are contiguous, for a generic time series
decomposed in this manner.

####  II.3.2 Developing “minimal models” for U.S. state trajectories

In the course of performing trajectory reconstructions, we verified that those
reconstructions built using just the first, $k=1$ principal component
collapsed to a single curve. In other words, PCA lead to a description of the
system that was identical and equivalent for all $n$ states, at this
“coarsest” level of resolution. Isolating just the most prominent peaks within
this curve (“find_peaks” function, signal processing toolbox, _SciPy_ v.1.4.1,
with the “prominence” parameter set to the minimum-observed value of that
single, $k=1$ reconstruction curve, 0.05) allowed us to _compress_ our “$k=1$
reconstruction” into a handful of values representing its most pronounced
excursions above its minimum. The dates associated with these maxima were then
cross-referenced with those of the greatest overall “importance,” according to
the first $k_{0}$ components cumulatively (see Sec. III.3.3).

Since the $k=1$, reconstructed trajectory represented a single, unifying
summary of all $n$ original trajectories, we referred to this curve as the
“master model.” The master model served as our _minimal_ model of the
archetypal daily new infection trajectory, for an arbitrary U.S. state.

We quantified the errors associated with assuming that all states followed a
given master-model behavior exactly. These included _i)_ the extent to which
the full, $n=50$ master model failed to capture nuances in the $n$ individual
trajectories from which it was inferred, and _ii)_ how well the aforedescribed
compression, consisting exclusively of the master model _maxima_ , captured
information conveyed by the master itself. To do so, we computed the mean
absolute deviations between _i)_ the 446-day, $n=50$ master and all the
individual state trajectories used in its creation, and _ii)_ the identified
maxima in the master and the corresponding trajectory heights at the
associated dates for the individual states. For both cases, we report the
averages of these deviations, across all included states, in Section III.3.3.

To characterize improvements afforded by invoking our dynamical clustering, we
repeated this process – creation and compression of a master model, and an
evaluation of the errors accrued when naïvely assuming this archetypal
behavior for all the states it summarizes – separately for each cluster, for
both the 446-day and maxima-only cases.

###  II.4 Significance of stringency and vaccination trajectory differences
among clusters

_Averages_ of government policy stringency trajectories and vaccination-rate
trajectories (introduced in Section II.1) were created for each cluster having
at least two members by taking the means, across all cluster members, over
time. That is, we separated states according to a learned cluster partitioning
$\mathcal{C}$ and computed the mean, standard deviation and standard error (of
the mean) over the constituent states for each cluster separately, by day.

To develop a framework for discerning whether the $\sim|\mathcal{C}|$
resulting policy (or vaccination-rate) “average trajectories” exhibited
statistically significant differences from one another – as well as to reduce
overfitting and more reasonably approximate sample independence – we chose not
to compare from all daily stringency values, but rather made the approximation
that the trajectory values from the _first day of each month_ could themselves
serve as independent observations on a given cluster-system. We then performed
a Krusal-Wallis test, with the number of samples equal to the number
$\sim|\mathcal{C}|$ of clusters with at least two state members, drawing
observations from the first day of each month between May 2020 and May 2021
(inclusive, so as to associate a total of exactly 13 observations with each
group).

As detailed below, vaccination trajectories can be shown to become (visually)
distinguishable only after a certain calendar date in the history of the
pandemic. We first performed a dynamical clustering above for all vaccination-
rate trajectories in the manner described in Sec. II.2 for infection
trajectories; we then created “average vaccination trajectories” for all the
(previously identified) infection-trajectory clusters with at least two states
as members, and identified where the vaccination-trajectory slopes were
changing by studying their first and second derivatives (built-in “ _diff_ ”
function in Python 3’s _NumPy_ and _Pandas_). Finally, we performed a Kruskal-
Wallis test to determine whether the vaccination-rate values on July 22, 2021
– a date by which the rates for all states were as well-separated as they
would be by the end of our data – differed significantly, between the
aforementioned infection-derived clusters.

We report both test statistics $H$ and $p$-values for $5\%$ significance in
Section III.3.4, as well as the results of a post-hoc test (Python
implementation of Conover’s test as part of the scikit-posthocs package, using
“holm” step-down method with Sidak adjustments to $p$-values) to infer which
of the distinct pairs of groups are statistically distinguishable.

###  II.5 Relationship between vaccination rates and emerging, new-case rates

All the methods discussed to this point have been aimed at establishing that
all $n=50$ U.S. daily new infection trajectories fall into a low-dimensional
set, $\mathcal{C}$; in Section III, we infer an appropriate value for
$|\mathcal{C}|$ (i.e., a partitioning that sufficiently distinguishes coarse
differences in “wave” histories). We also wish to assess the _usefulness_ of a
cluster set $\mathcal{C}$ for efficiently summarizing U.S. history – and its
potential grounding in reality, beyond its practical use for data compression
– according to how it recapitulates patterns in geography, government
stringency, and vaccination rates.

In Section III.3.4, we see the vaccination-rate trajectories – percentages of
fully-vaccinated individuals represented in each U.S. state’s population –
diverge from one another markedly (and then remain divergent) only following
the late Spring of 2021. Furthermore, those vaccines against COVID-19
available in the period described by our study – April 2020 to July 2021 –
require several weeks’ time to accomplish even partial seroconversion Lombardi
et al. (2021); Sette and Crotty (2021). We therefore studied the degree to
which the average vaccination rates within each cluster, _approximately two
weeks before the end_ of our data, predicted daily-new-case statistics in the
corresponding clusters on _the last full day_ of data.

We performed tabulations of both the July 22, 2021 vaccination rates and the
August 3, 2021 “daily new infection” records, by cluster, according to the set
$\mathcal{C}$ learned and fixed in Section III.2.1. We then investigated the
relationship between within-cluster vaccination rates and _i)_ the within-
cluster _mean_ of states’ infection trajectory values and _ii)_ the within-
cluster _standard error of the mean_ for these (relative) new-case loads, at
an approximate two-week lag, by means of both Spearman rank-order correlations
and two linear regressions. For both types of analyses, we considered only
states that belonged to clusters encompassing at least two states, omitting
“singleton” clusters that consisted of only a solitary state (see Section
III.2.1). We read off the output “ _r_ ”of the “ _linregress_ ” function in
_Scipy_ ’s Statistical functions module as the $r^{2}$ value, or linear
(Pearson) correlation coefficient, to further quantify linear trends.

##  III Results

###  III.1 Example data and relevant time frames

The state of New York exhibited a strong participation in the first, third,
and fourth wave opportunities (Fig. 1), with essentially no reprieve between
the latter two; New York was also among those states exhibiting a relative
quiescence during the second wave opportunity. The number of new infections
during the fourth wave opportunity (Spring 2021) were comparable to those
observed exactly one year prior, during the Spring of 2020. Only the initial
weeks of the fifth wave opportunity are depicted, since the infection data
source terminates on August 3, 2021.

Figure 1: Infection, vaccination, and policy trajectories for the state of New
York.

The vaccination trajectory for New York begins rising after the new-infection
peak in the Winter of 2020 - 2021, coincident with the onset of vaccine
administration to the public. It changes slope appreciably only during the
daily new infections peak that was realized in the context of the fourth wave
opportunity (April 2021). Government policy stringency information exists
prior to the first wave opportunity, but ends prior to the fifth, in our data
sets.

###  III.2 Initial clustering and ocular analysis

####  III.2.1 State’s unique trajectories fall into 4 major clusters

As the cluster threshold is varied from its minimum (0) to its maximum
possible value (1), the initial hierarchical clustering of our infection
trajectories converges rapidly to resolve just a handful of large clusters. At
threshold 0.2, we resolve only 8 multi-state clusters (Fig. 2a); approaching
threshold 0.43 and beyond, clusters tend to grow in size, while
$|\mathcal{C}|$ decreases concomitantly. Once the threshold values reach
$0.6$, representing a temporal cross-correlation of at least
$1-0.6\rightarrow\mbox{}40\%$ between any pair of constituent trajectories
internal to a cluster, there remain only three big clusters and two outlying
states – California and Hawaii both form their own, standalone, singleton
“clusters.”

The 0.43 threshold represented the highest threshold for which non-singleton
clusters were still of comparable size. In particular, whereas the three large
clusters at the 0.6 threshold are composed of 12, 12, and 24 states, those at
threshold 0.43 have sizes 12, 12, 8, and 16. Although it is possible for the
distinct temporal patterns characteristic of the 50 states to fall into
differently-sized clusters, we expect that power can be maximized by comparing
the statistics of uniformly-sized clusters (the size distribution associated
with highest entropy Cover (1999)); by this heuristic, we chose to work
predominantly with the clusters resolved at this intermediate threshold value.
To visualize salient disparities in the states’ temporal dynamics that
separate the corresponding, $|\mathcal{C}|=4$ major clusters at this “locked”
threshold, we plot mean infection trajectories for each cluster, color-colored
by proximity in correlation-space (Fig. 2b); these color-coded data, taken
from the training set that ran from April 2020 to July 2021, were superimposed
upon gray curves representing the same mean trajectories, but updated Ref.
Repository of COVID-19 data underlying the analysis and visualizations created
by the Johns Hopkins Centers for Civic Impact for the Coronavirus Resource
Center (CRC) to run between March 2020 and August 2021.

####  III.2.2 Clusters differ in their wave timings and relative heights

Commonalities among the mean, cluster-wise trajectories include the
aforementioned Winter 2020 - 2021 maxima, as well as the presence of smaller
peaks in the Spring and Summer seasons. Although our correlation-based
clustering was intended primarily to address the _timing_ , or
synchronicities, of realized wave opportunities, it also highlights some
salient differences in the magnitudes, or _relative heights_ , of those
opportunities. The distinguishing properties follow.

Cluster 1 alone saw a major, initial surge of new cases in Spring of 2020 –
the first wave opportunity – while the states in Clusters 2 and 3 saw their
first collective spike only weeks later, during the second wave opportunity in
the Summer of 2020. Cluster 4 is distinguished from Clusters 1-3 by its
effective lack of either a Spring or a Summer 2020 surge, with the third,
Winter wave opportunity playing host to its first major elevation. At a finer
resolution, this principal, Winter maximum is observed to concentrate early
(vicinity of December) for Cluster 4, with only a modest resurgence of new
infections shortly thereafter; other clusters exhibited a different, more
multiply-peaked structure, as well as comparatively _later_ peaks, in that
season. The closest pair in correlation-space – Clusters 2 and 3 – share most
of their coarse features, with aspects of this multiply-peaked Winter setting
the latter apart from the former. All U.S. states were observed to have
partial resurgences in Spring 2021, and every single state experienced an
upswing in Summer 2021 (not explicitly depicted), with pronounced variability
even _within_ clusters by August (Fig. 2b).

Figure 2: Clustering results for the U.S. state infection trajectories. a)
Dendrogram summarizing hierarchical clustering. b) Trajectory mean (standard
error) within each cluster; any discrepancies between training data source
(color) and testing data (grayscale) are due to the use of independent
sources. Outliers (California and Hawaii) are not depicted. c) Geographical
distribution of cluster memberships at the different thresholds indicated by
gray shading in top panel.

Collectively, these patterns in the realized heights for each wave opportunity
– that is, the (relatively) medium-sized first wave, essentially absent second
wave, and very large third wave expressed by Cluster 1; then, the absent first
wave, medium-sized second wave, and very large third wave for Clusters 2 and
3; and then, the absent first and second waves, but large third wave for
Cluster 4 – provide some degree of intuitive explanation for why these
particular clusters were resolved. We quantify this intuition by focusing on
specific, “important dates,” in Section III.3.3.

###  III.3 Correlative and Predictive Cluster Validations

####  III.3.1 Dynamical clusters turn out to be geographically related

Even given the robust set of qualitative dynamical patterns teased out by
means of ocular analysis in the previous section, the learned partitioning
$\mathcal{C}$ still eludes practical explanation: are there other, non-
dynamical patterns associated with the states that were grouped together that
recapitulate the learned clustering? In order to visualize the spatial
arrangement of those states whose trajectories clustered together at various
thresholds, as one such pattern, we colored the states on a U.S. geopolitical
map, according to cophenetic similarity (correlation-space distance) as in
Fig. 2b.

Adjacent, or neighboring, U.S. states were much more likely to be clustered
together (Fig. 2c). With a few exceptions (including southern states, at the
0.43 threshold, that are separated by land but still connected by water), the
majority of the clusters at all three thresholds unite states that lie along
relatively uninterrupted, existing geographical paths. As those clusters
resolved at small thresholds merge to form larger clusters at larger
thresholds, the latter begin to reveal clear correspondences to established
U.S. geopolitical regions; these are visible at our “locked,” 0.43 threshold.

In particular, our Cluster 1 predominantly describes states that are
traditionally subsumed under “the Northeast”. Cluster 2 describes “the South,”
encompassing both the coastal southeastern border and several states along the
southern rim of the country. Cluster 4 represents “the Midwest,” only
exchanging the states of that region’s traditional, southernmost border for
the three on its western border. Notably, the two outliers – California and
Hawaii – are also _geographical_ extremes. The remaining Southern, Western,
and Central form one, comparably-sized cluster that (with the exception of
Colorado) sweeps a clean arc from roughly the latitudinal center of the
country to Washington State; this “Rest” of the states, comprising Cluster 3,
is geospatially contiguous (and, at higher thresholds, merges) with its most
proximal neighbor in correlation-space, Cluster 2.

It might be reemphasized that all clustering operations were done
_dynamically_ – that is, according to correlative similarities between the
respective time courses of “daily new case” loads associated with each U.S.
state. No explicit spatial or geographical clustering was considered; rather,
any geographical patterns emerged as _consequences_ of the dynamical
clustering described in Section IIIb.

####  III.3.2 Date ranges most explanatory for telling states apart coincide
with peak times for whole-U.S. infection “waves”

Analyzing the (cumulative) percentages of variability that can be “explained”
by incorporating exclusively the first $k$ principal components in
reconstructions of the daily new infection trajectories (Table I), it is
apparent that $k_{0}=6$ principal components suffice to capture over $90\%$ of
the total variance considered by PCA. Reconstructions using the first $k=10$
components produce trajectories almost indistinguishable from the originals
(not depicted).

| Principal Component
---
Number
$k$
| Cumulative
---
% of Variance
Explained
| Reconstruction
---
error for first
$k$ components
|
1 | .437 | .096 |  |
2 | .646 | .070 |  |
3 | .746 | .055 |  |
4 | .831 | .048 |  |
5 | .871 | .040 |  |
6 | .908 | .035 |  |
7 | .934 | .029 |  |
8 | .949 | .025 |  |
9 | .962 | .022 |  |
10 | .970 | .019 |  |
Table 1: Key results of PCA.

Feature importance analysis suggests that the key dates for explaining
trajectory variance encompassed long periods from November to December 2020,
and January to mid-February 2021, with weaker contributions in the Spring of
2021 and even late Summer 2020 – that is, all the “important dates” comprised
continuous _ranges_ – for the $k=1$ component (Fig. 3a). The same is true for
the successive components $k\in\\{2,3,4,5\\}$. For instance, the high-
prominence ranges for the $k=2$ component coincide with the time-continuous
Summer 2020 and Spring 2021 waves. Two lower-prominence peaks in the $k=2$
importance plot describe a more subtle “fine-tuning” to temporally briefer,
Winter fluctuations.

The overall importances imputed for each date, based on the first $k_{0}=6$
components, also manifest as continuous stretches in time. Above-average
values spanned only three, broad date ranges (Fig. 3b) across the second,
third, and fourth wave opportunities. These three ranges all corresponded
closely to established date ranges for whole-U.S. peaks in daily new infection
records (see discussion of the master curve in Sec. III.3.3 below, or Ref.
Hallas et al. (2020)), implying that the time periods surrounding _global_ ,
national peaks, were also those of the greatest _local_ , state-level
variation.

Figure 3: Results of applying PCA to temporal features (daily new infection
values) of the $n=50$ U.S. state trajectories. a) Feature importances for
individual PCA components. b) Cumulative ($k=1,\dots,6$), “overall” feature
importances. c) Original $n=50$ U.S. daily new infection trajectories, with
master model superimposed in black and “important” dates for explaining
variation as red vertical lines. d) Mean policy and vaccination trajectories,
by cluster, for context. Gray boxes indicate regions where daily values of the
curve reconstructed from the first $k_{0}$ eigenvectors exceeded their mean
across all encompassed days.

####  III.3.3 Clustering improves PCA-derived, minimal models of archetypal
U.S. infection trajectories

Projecting the $n$ individual state trajectories onto just the first, $k=1$,
component resulted in the same, archetypal curve for all states. This curve
was identified with the _mean trajectory_ over the $n=50$ time series,
demonstrating that PCA had ultimately learned to distinguish states by their
_deviations_ from the collective, average national behavior.

Excluding the Spring of 2020 for the reasons discussed in Sections II.1.1 \-
II.1.2, all the “wave opportunities” are represented by maxima in the master
curve (red vertical lines in Fig. 3c). At the prominence level used to detect
peaks here (Section II.3.2), the multiply-peaked, third wave opportunity
splits into two separate maxima for consideration.

Separating the $n$ states by cluster and inferring a master model for each
cluster separately recovers the “average” trajectories for each cluster (see
Fig. 2b). Reconstructing individual state trajectories using their respective
master models can reduce reconstruction error, from at least two perspectives:
1) the compressed, maxima-only representation and the full, 446-day
trajectories should be reconstructed more accurately when a master-model
construction is done _separately_ for each cluster, as opposed to all $n$
states _together_ (mean absolute deviations between the states’ original and
reconstructed trajectories should tend to decrease); 2) the error between an
original trajectory and the “tailored” master made for its containing cluster
should diminish, even if no explicit reconstruction is done (the mean absolute
deviations between any cluster-averaged trajectory and the original trajectory
for given state in the corresponding cluster should be smaller than that
between the original trajectory and the full, “$n=50$” master model).

Measurements for the latter (Table II), based on either the full, 446-day
curves and the 4 key, “important” dates discerned above confirm that errors
are smaller within clusters, but also suggest that encoding the 446-day
trajectories as a mere 4-vector of amplitudes at the important dates
accomplishes a fair _compression_ of states’ infection histories: the latter
are only slightly higher. Since all error values are of order $\sim 0.1$, and
cluster memberships are the only inputs needed to a create master model
(beyond the data themselves), specifying the cluster membership – or,
equivalently, its geographical region, according to Fig. 2c – amounts to a
prediction of that state’s normalized trajectory heights at each wave
opportunity with collective, net error of $10\%$ of the height of the
corresponding Winter peak, on average.

| Relevant
---
master model
for evaluation
| Error (MAD),
---
based on full
master model
| Error (MAD),
---
based on peaks
of master model
| Difference in
---
(MAD) error,
full $\rightarrow$ peaks
|
All states ($n=50$) | .0962 | .1514 | .0552 |  |
Cluster 1 ($n=12$) | .0624 | .0933 | .0309 |  |
Cluster 2  ($n=8$) | .0514 | .0694 | .0180 |  |
Cluster 3 ($n=16$) | .0611 | .0882 | .0271 |  |
Cluster 4 ($n=12$) | .0562 | .0906 | .0344 |  |
Table 2: Errors associated with average, state-wise deviation from master
models

####  III.3.4 Government stringencies & vaccination rates for different
clusters have different average behaviors over time

The geographical considerations above (Fig. 2c) support the notion that there
might be “real-world” correspondences to the four, abstracted clusters
identified in Sec. III.2.1. A real-world relevance for these clusters is
further supported by their ability to find discriminatory patterns in two
other metrics, neither of which was included in their definition: government
policy stringency scores and vaccination rates (shown for temporal context in
Fig. 3d; see Section II.1.3).

For all Clusters 1 - 4, the Spring of 2020 begins with a rapid increase in
stringency that reaches an all-time high around April, despite the fact that
not all states realized an equally strong “first wave” opportunity (see Fig.
2b). It is clear that _i)_ average stringency scores are highest for the
Northeast, and lowest within the Midwest, through time; _ii)_ the remaining
clusters are less distinct from each other than they are from Cluster 1;
_iii)_ Clusters 2 and 3, the most overlapping, sustained stringencies
_intermediate_ to those of the Northeast and Midwest for the majority of time.

Repeating this analysis for the vaccination-rate trajectories again highlights
the Northeastern states as the highest-scored cluster, and reemphasizes that
high degree of similarity between the Southeastern rim and the “Rest” of the
Southern, Central and Northwestern regions that we had initially established
via infection-trajectory clustering. In contrast to the overall pattern
observed for stringencies, these latter groups are the lowest-ranked for
vaccination.

Statistical hypothesis testing (four-sample Kruskal-Wallis, with 13
observations in each group) suggests that the monthly values for the mean
stringency scores ($H=17.6,~{}p=5\cdot 10^{-4}$) and vaccination rates
($H=31.9,~{}p=5\cdot 10^{-7}$) are _not_ identical across the four clusters.
Adjusted $p$-values for the various, possible pairings suggest that Cluster 1
was distinct from both Clusters 2 and 3 ($p=2\cdot 10^{-2}$) as well as from
Cluster 4 ($p=5\cdot 10^{-5}$), but that Clusters 2 and 3 were statistically
indistinguishable from each other ($p=9\cdot 10^{-1}$), and possibly even from
Cluster 4 ($p=1\cdot 10^{-1}$).

The vaccination trajectories can be observed to diverge appreciably only
around the start of April, and even then a hierarchical, correlation-based
clustering (Section II.2) unifies all $n=50$ trajectories at a “distance” of
$\sim 0.025$, a full order of magnitude smaller than that of the most tightly-
correlated merge (Cluster 2 in Fig. 2b) among the infection trajectories.
Also, numerical differentiation (not depicted) reveals that the only major,
dramatic slope changes occur around April 1 and May 1, 2021, so that later
divergences are largely explained via higher vaccine “adoption rates,” on
average, in the early stages of availability among states belonging to certain
clusters. All this evidence for a “low distinguishability” through the Spring
of 2021 and an apparent separability only later, in Summer 2021, was
consistent with the Kruskal-Wallis test for cluster-wise differences among
July 22, 2021 vaccination rates: $H=20.6,~{}p=1\cdot 10^{-4}$.

Among the four clusters containing at least two states each, Conover’s post-
hoc tests suggests that Cluster 1 alone was unambiguously distinguishable from
all the others (Cluster 2, $p=2\cdot 10^{-4}$; Cluster 3, $p=1\cdot 10^{-5}$;
Cluster 4, $p=9\cdot 10^{-3}$), and that Clusters 2 and 3 were completely
indistinguishable from one another ($p=0.968$). Post-hoc testing did not
indicate strong degrees of separation between Cluster 4 and Cluster 2
($p=0.234$) or Cluster 3 ($0.170$).

The fact that vaccination-rate trajectories do not differ significantly until
the concluding region of our time series indicates that the only vaccination
values potentially useful in discriminative _predictions_ are those from
Summer 2021.

###  III.4 Seatbelt Hypothesis: Within-cluster vaccination rate predicts
uncertainty in future new-case upticks

Correlation and regression analyses suggest monotonic relationships between
our cluster mean vaccination rates and several aspects of the “average
infection trajectories” associated with the corresponding clusters, at a two-
week lag. First, the mean values of the daily new infections in any given
cluster at the beginning of August, relative to the corresponding Winter 2021
peak values, are inversely related to the mean vaccination rate in that
cluster (Spearman rank correlation coefficient $-1$; Fig. 4). Emphatically,
this relationship is not strictly linear. For instance, several Midwestern
states – the cluster with the _second-lowest_ mean in daily new cases in early
August – recorded fewer cases than all the Northeastern states, the cluster
with the _lowest_ mean in daily new cases in August. Still, the Pearson
correlation coefficient associated with the linear-order, pairwise correlation
between mean vaccination rates and mean case loads at a two-week lag is
measured (via regression or direct computation) as $-0.8$, suggesting a
directional trend.

Beyond the means, the average, within-cluster July 22 vaccination rate also
closely tracks another quantity – namely, the _standard error_ of the mean in
the daily new infection trajectory values at the beginning of August.
Specifically, the fraction of fully-vaccinated individuals is proportional to
the “spread” of – or, “uncertainty” in – trajectory values realized across the
clusters’ constituent states. Concretely, for the most _highly-vaccinated_
group of states (Cluster 1, the Northeast) the range of observed infection-
trajectory upticks is the _smallest_ , and this range tends to increase with
decreasing vaccination coverage (Spearman rank correlation coefficient of
$-1$). Fitting a linear model to measure the Pearson correlation coefficient
(value $-0.74$) reveals a strongly negative trend: vaccination rates recorded
in late July are _positively_ correlated with a measure of case-load
predictability in early August (not depicted in Fig. 4).

Figure 4: “Seat belt” hypothesis: organizing the clustered regions in order of
decreasing vaccination rate reveals negative relationships with the mean and
standard error of new, emerging cases rates, at a two-week, lag, during Summer
2021.

Consider, as an analogy, the adoption of seat belts as a measure of protection
against severe injury. Seat belts do not insure against the acquisition of all
wounds, but are designed to be good at decreasing the probability of
fatalities. One cannot guarantee survival in the worst of accidents by using a
seat belt, nor is one ever bound to acquire severe injuries when failing to
activate the device – yet the _range_ of the appreciably likely injuries that
one might receive, in the context of common accident configurations, can be
constrained via wearing a seat belt. Similar is the pattern suggested by Fig.
4: among those states belonging to clusters for which vaccination was more
widespread in late July, the _range_ of relative trajectory heights two weeks
later, at the beginning of August, is observed to be smaller. In that short
time period during which vaccination trajectory values differed significantly
among state clusters, these values were also predictive of yet unseen,
“future” infection rates.

##  IV Discussion

Here, we demonstrate that dynamical clustering allows separation of all 50
states into a small number of clusters (we settled on four, with two
outliers). Meanwhile, Principal Components Analysis corroborated the existence
of several distinct “waves” of COVID-19 infections nationally. These
encompassed four key dates, which served in to build accurate compressions of
cluster “average trajectories” with little other information, showing that a
low-dimensional model can capture much of the variation in relative COVID-19
case loads over time.

The fact that our set of clusters $|\mathcal{C}|$ happened to fall into
established geographical regions, and showed correlations to state
stringencies and vaccination rates further supports the potential for the
clusters we discovered to be accurately reflecting reality. Models trained
using states’ data through June 2021 were useful in predicting standard
deviations of cluster case loads in the subsequent, “fifth wave” from July
through September 2021. One parsimonious explanation is that the vaccination
rates do correlate with case loads, but that greater usefulness lies in their
(positive) correlation with the _predictability_ of case loads in our period
of interest. Based on all these results, we developed the “seat belt”
hypothesis, founded on the comparison that, while seat belts don’t prevent
accidents, they do restrict the likely range of adverse outcomes. In this
vein, high vaccination rates may not prevent COVID-19 from affecting a state,
but states with lower vaccination rates might expect a higher potential
_ceiling_ for their daily new case loads.

Our results are consistent with previous analyses that identified four,
historical waves Hallas et al. (2020), and a role for geography in identifying
states with analogous pasts Times (2020), although our findings recapitulated
both without explicit geographic input to our models. Nevertheless, there are
several important considerations for interpreting our own results in
comparison with other analyses. First, we have normalized all case rates to
one peak at the end of 2020, so that our results relate to _relative change_
and not _absolute case loads_. That is, we clustered states not based on how
many cases they have recorded, but the synchronous changes-of-shape in their
infection patterns over time. Additionally, we developed models only of
reported _cases_ , and different results might have arisen had we focused on
other metrics, as in deaths or hospitalizations. Similarly, many more data
sources could conceivably be added to this model and might change the specific
predictions. We chose not to pursue a more exhaustive data-scraping approach,
instead choosing to focus on the capacity to make useful predictions with a
simple model without either wide-spread exploration or explicit notions of
causality, both of which could soak up arbitrary amounts of future time and
energy.

As yet more new data becomes available, it will be possible to “test” the
hypothesis that higher vaccination rates will translate to a restricted
“likely range” of new case loads generally. Here, too, future researchers must
exercise caution. For example, whereas our vaccination-rate trajectories
remain similar over long stretches of time, a more exhaustive analysis would
have to rule out that infection dynamics do not bifurcate into qualitatively
different behaviors if local vaccination rates reach a certain (even sub-herd
immunity) value. On this note, we emphasize that our state- and state cluster-
level results might be superseded by different insights at, for instance, the
level of U.S. counties. In addition to these considerations, as well as the
overwhelming likelihood that vaccination prevalence alone does not control
new-infection dynamics, there remains the possibility that other complex
influences – such as geographic spread from the COVID-19 pandemic” (9/22/20)
or climate patterns Bukhari et al. (2020) could conceivably translate to
different clusters “taking the lead” in emerging, new-case rates at different
times. We stress that our “seat belt” measurements were evaluated while U.S.
states were experiencing a clear upswing, ostensibly due to the spread of a
novel SARS-CoV-2 variant.

We emphasize again that, while our models do seem to provide useful
predictions and to correspond to other real-world features, they are not
causal models. The fact that state’ dynamics can be identified and tied into
clusters with few inputs is efficient, but does not resolve the issue of _why_
compact patterns exist. Explorations along these lines do remain important,
and we hope our model might be useful enabling others to add more variables in
search of causality. Effectively infinite other data sources could be explored
in this endeavor, but given the behavior of the clustering we observed, some
possibilities suggest themselves, including mask-use policies and climate.
While we do not assess political data University (2021) explicitly in our
model, the regional clustering and relative differences in stringency suggest
that they maybe worth exploring in these contexts as well.

While our pipeline is not the only existing approach to compressing infection-
rate information in such a way that is useful for prediction, it offers
distinct advantages. Several analyses elaborated previously, as in Ref Hallas
et al. (2020), focused predominantly upon tracking and categorizing variations
in state-level policies, and provided the basic terminology for describing the
three completed “waves” of COVID-19 cases in the United States (as well as
“concerns of a fourth wave” at the time of its publication in Spring 2021)
Hallas et al. (2020). This four-wave structure was discerned from observations
on the overall patterns of COVID-19 cases across the U.S (the sum of
contributions from all 50 states); instead, we _impute_ a data-driven number
of waves, and find that the encompassed dates agree well with those of Ref.
Hallas et al. (2020).

This latter sentiment is also echoed within two earlier studies that had
explored multiple clustering methods James and Menzies (2020a) and attempted
to classify U.S. states as being within their first or second “surge” James
and Menzies (2020b). Unlike our work here, the “optimality” of a given
partitioning $\mathcal{C}$ in Ref. James and Menzies (2020a) is not evaluated
on whether clusters predict external, ancillary variables, but by analytic
means. The analyses in Ref. James and Menzies (2020b) also focused on “daily
new cases,” which we chose over other metrics, such as cumulative cases or
daily death rates, as a snapshot of the evolving transmission situation within
the U.S. The authors of the present work were unaware of both works while
conducting the research described here, but our contribution can nonetheless
be considered an extension or generalization of the former: aside from the
benefit of a somewhat more retrospective look at the progression of the
pandemic, our emphasis on external validation provides for more direct
application of clustering to future decisions and research.

In conclusion, our model supports the idea that gross fluctuations the U.S.
states’ COVID-19 case rates are more regionally coordinated than is revealed
by a focus on raw case numbers alone. As a result, while others have reported
occasions on which California and Florida may have seen similar statistics
Curley (2021); Bump (2021), our findings suggest that a neighbor like Georgia
would be more likely to provide actionable information for what is likely to
happen to Florida (especially in the context of a new wave) than would
California or New York. It remains to be shown, in future work, whether these
same states would unite and cluster together on other grounds (e.g.,
demographic, political, or climatological), and if the specific groupings
discovered here serve useful in future infectious disease containment or
mitigation strategies.

We have demonstrated that the fifty U.S. “wave” histories fall into just four
distinct infection patterns, and that each pattern can be associated with a
specific geographical region. Grouping states by infection pattern revealed
that within-group vaccination rates predict key new-case statistics during the
early stages of the Summer 2021 surge.

###### Acknowledgements.

The authors would like to extend a special thanks to the CaliBaja Center for
Resilient Materials and Systems, and to those at the San Diego Supercomputer
Center, for hosting the ENLACE 2021 program.

## References

  * Van Dorn et al. [2020] Aaron Van Dorn, Rebecca E Cooney, and Miriam L Sabin.  COVID-19 exacerbating inequalities in the US.  _Lancet (London, England)_ , 395(10232):1243, 2020\.
  * Chaudhry et al. [2020] Rabail Chaudhry, George Dranitsaris, Talha Mubashir, Justyna Bartoszko, and Sheila Riazi.  A country level analysis measuring the impact of government actions, country preparedness and socioeconomic factors on OVID-19 mortality and related health outcomes.  _EClinicalMedicine_ , 25:100464, 2020.
  * Ciminelli and Garcia-Mandicó [2020] Gabriele Ciminelli and Sílvia Garcia-Mandicó.  COVID-19 in Italy: An Analysis of Death Registry Data.  _Journal of Public Health_ , 42(4):723–730, 2020\.
  * Weissman et al. [2020] Gary E Weissman, Andrew Crane-Droesch, Corey Chivers, ThaiBinh Luong, Asaf Hanish, Michael Z Levy, Jason Lubken, Michael Becker, Michael E Draugelis, George L Anesi, et al.  Locally Informed Simulation to Predict Hospital Capacity Needs During the COVID-19 Pandemic.  _Annals of internal medicine_ , 173(1):21–28, 2020.
  * Sen-Crowe et al. [2021] Brendon Sen-Crowe, Mason Sutherland, Mark McKenney, and Adel Elkbuli.  A Closer Look Into Global Hospital Beds Capacity and Resource Shortages During the COVID-19 Pandemic.  _journal of surgical research_ , 260:56–63, 2021.
  * Alfaro et al. [2021] Sophie Alfaro, Brendon Sen-Crowe, Mark McKenny, and Adel Elkbuli.  A closer look at U.S. COVID-19 vaccination rates and the emergence of new SARS-CoV-2 variants: It’s never late to do the right thing.  _Annals of Medicine and Surgery_ , 69:102709, 2021.
  * Chadha et al. [2021] Jatin Chadha, Lavanya Khullar, and Nidhi Mittal.  Facing the wrath of enigmatic mutations: a review on the emergence of severe acute respiratory syndrome coronavirus 2 variants amid coronavirus disease-19 pandemic.  _Environmental microbiology_ , 2021.
  * Loomba et al. [2021] Rohit S Loomba, Gaurav Aggarwal, Saurabh Aggarwal, Saul Flores, Enrique G Villarreal, Juan S Farias, and Carl J Lavie.  Disparities in case frequency and mortality of coronavirus disease 2019 (COVID-19) among various states in the United States.  _Annals of medicine_ , 53(1):151–159, 2021.
  * Hamman [2021] Mary K Hamman.  Disparities in COVID-19 mortality by county racial composition and the role of spring social distancing measures.  _Economics & Human Biology_, 41:100953, 2021.
  * Xian et al. [2020] Zhaoying Xian, Anshul Saxena, Zulqarnain Javed, John E Jordan, Safa Alkarawi, Safi U Khan, Karan Shah, Farhaan Vahidy, Khurram Nasir, and Prachi Dubey.  Racial and Ethnic Disparities in COVID-19 Infection and Mortality in the United States: A state-wise update.  _medRxiv_ , 2020.
  * Solís et al. [2021] Patricia Solís, Gautam Dasarathy, Pavan Turaga, Alexandria Drake, Kevin Jatin Vora, Akarshan Sajja, Ankith Raaman, Sarbeswar Praharaj, and Robert Lattus.  Understanding the spatial patchwork of predictive modeling of first wave pandemic decisions by US governors.  _Geographical Review_ , 111(4):592–615, 2021\.
  * for Disease Control et al. [2021] Centers for Disease Control, Prevention, et al.  Covid data tracker weekly review.  _Centers for Disease Control and Prevention_ , 2021.
  * Curley [2021] Christopher Curley.  Why California, Florida Have Similar Number of COVID-19 Cases, 2021.  URL https://www.healthline.com/health-news/why-do-california-and-florida-have-similar-covid-19-case-rates-the-answer-is-complicated.
  * News [2021] CBN News.  Freedom vs. Lockdowns: Data Shows FL, CA With Nearly Identical Outcomes in COVID Cases Despite Opposite Approach, 2021\.  URL https://cmsedit.cbn.com/cbnnews/us/2021/march/freedom-vs-lockdowns-data-shows-fl-ca-with-nearly-identical-outcomes-in-covid-cases.
  * Strogatz [2018] Steven H Strogatz.  _Nonlinear dynamics and chaos with student solutions manual: With applications to physics, biology, chemistry, and engineering_.  CRC press, 2018.
  * Repository of COVID-19 data underlying the analysis and visualizations created by the Johns Hopkins Centers for Civic Impact for the Coronavirus Resource Center [CRC] Repository of COVID-19 data underlying the analysis and visualizations created by the Johns Hopkins Centers for Civic Impact for the Coronavirus Resource Center (CRC).  John Hopkins University Coronavirus Resource Center.  https://github.com/govex/COVID-19, 2021.  [Online; accessed 15-July-2021].
  * Wu et al. [2020] Sean L Wu, Andrew N Mertens, Yoshika S Crider, Anna Nguyen, Nolan N Pokpongkiat, Stephanie Djajadi, Anmol Seth, Michelle S Hsiang, John M Colford, Art Reingold, et al.  Substantial Underestimation of SARS-CoV-2 Infection in the United States.  _Nature communications_ , 11(1):1–10, 2020.
  * Hallas et al. [2020] Laura Hallas, Ariq Hatibie, Saptarshi Majumdar, Monika Pyarali, Rachelle Koch, Andrew Wood, and Thomas Hale.  Variation in US states’ Responses to COVID-19.  _Blatavnik School Working Paper Series_ , 2020.  URL https://www.bsg.ox.ac.uk/research/publications/variation-us-states-responses-covid-19.
  * Repository of COVID-19 data curated by Our World In Data [https://ourworldindata.org/] Repository of COVID-19 data curated by Our World In Data (https://ourworldindata.org/).  COVID-19 Dataset by Our World in Data.  https://github.com/owid/covid-19-data, 2021.  [Online; accessed 15-July-2021].
  * Review [2021] World Population Review.  Us States - Ranked by population 2021, 2021.  URL https://worldpopulationreview.com/states.
  * Allen et al. [2021] Hester Allen, Amoolya Vusirikala, Joe Flannagan, Katherine A Twohig, Asad Zaidi, Dimple Chudasama, Theresa Lamagni, Natalie Groves, Charlie Turner, Christopher Rawlinson, et al.  Household Transmission of COVID-19 Cases Associated with SARS-CoV-2 Delta Variant (b. 1.617. 2): National case-control study.  _The Lancet Regional Health-Europe_ , page 100252, 2021.
  * community [2021] The SciPy community.  scipy.spatial.distance.pdist, 2021.  URL https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html#scipy.spatial.distance.pdist.
  * Lombardi et al. [2021] Andrea Lombardi, Giorgio Bozzi, Riccardo Ungaro, Simone Villa, Valeria Castelli, Davide Mangioni, Antonio Muscatello, Andrea Gori, and Alessandra Bandera.  Mini Review Immunological Consequences of Immunization With COVID-19 mRNA Vaccines: Preliminary Results.  _Frontiers in Immunology_ , 12:677, 2021.
  * Sette and Crotty [2021] Alessandro Sette and Shane Crotty.  Adaptive Immunity to SARS-CoV-2 and COVID-19.  _Cell_ , 2021.
  * Cover [1999] Thomas M Cover.  _Elements of information theory_.  John Wiley & Sons, 1999.
  * Times [2020] The New York Times.  Coronavirus in the U.S.: Latest Map and Case Count, 2020.  URL https://www.nytimes.com/interactive/2021/us/covid-cases.html.
  * from the COVID-19 pandemic” [9/22/20] Anthony Fauci: “Insights from the COVID-19 pandemic” (9/22/20).  Youtube, 2020.  URL https://www.youtube.com/watch?v=J38QN1vnSO0.
  * Bukhari et al. [2020] Qasim Bukhari, Joseph M Massaro, Ralph B D’agostino, and Sheraz Khan.  Effects of Weather on Coronavirus Pandemic.  _International journal of environmental research and public health_ , 17(15):5399, 2020.
  * University [2021] Oxford University.  Oxford Covid-19 Government Response Tracker (OxCGRT), 2021\.  URL https://github.com/OxCGRT/covid-policy-tracker.
  * James and Menzies [2020a] Nick James and Max Menzies.  Cluster-based Dual Evolution for Multivariate Time Series: Analyzing COVID-19.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 30(6):061108, 2020a.
  * James and Menzies [2020b] Nick James and Max Menzies.  COVID-19 in the United States: Trajectories and Second Surge Behavior.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 30(9):091102, 2020b.
  * Bump [2021] Phillip Bump.  One year in, how four states fared in combating the coronavirus, 2021\.  URL https://www.washingtonpost.com/politics/2021/03/09/one-year-how-four-states-fared-combatting-coronavirus/.


# Energy transport in harmonically driven segmented Frenkel-Kontorova lattices

M. Romero-Bastida  SEPI ESIME-Culhuacán, Instituto Politécnico Nacional, Av.
Santa Ana No. 1000, Col. San Francisco Culhuacán, Delegación Coyoacan,
Distrito Federal 04430, Mexico<EMAIL_ADDRESS>
###### Abstract

In this work we study the energy transport in a one-dimensional system
composed of two dissimilar Frenkel-Kontorova lattices connected by a time-
modulated coupling and in contact with two heat reservoirs operating at
different temperature by means of molecular dynamics simulations. There is a
value of the driving frequency at which the heat flux takes its maximum value,
a phenomenon termed as thermal resonance. Structural modifications in the
lattice strongly alter the way in which the external driving interacts with
the phonon bands. The overlap of the latter in the harmonic regime of the
model determine the frequency range wherein resonance emerges. Parameter
dependencies by which the incoming heat flux can be directed to either of the
heat reservoirs are examined as well. Our results may be conductive to further
developments in designing thermal devices.

###### pacs:

44.10.+i; 05.60.-k; 05.45.-a; 05.10.Gg

##  I Introduction

Structural stability in mesoscale and nanoscale devices is closely related to
heat generation within. Experiments measuring the heat generated in electric
current-carrying metal-molecule junctions determined that the generated heat
can be substantial Tsutsui08 ; Huang07 ; Tsutsui19  and can therefore threaten
the device’s integrity. Efficiently dissipating heat in such devices is thus
important and a problem that must be considered, especially with regard to
technological applications. An example of devices wherein the above issues are
relevant is that of microelectromechanical and nanoelectromechanical systems.
These are being developed for a host of nanotechnological applications, such
as highly sensitive mass Yang06 , spin Rugar04 , and charge detectors Andrew98
, as well as for basic research in the mesoscopic physics of phonons Volz16 .
In fact, an improved understanding of the manipulation and control of phonons
—that manifest themselves as heat at the nanoscale level— is necessary for
further progress in addressing the above mentioned problem Li12 .

Several models and mechanisms have been proposed to control or manipulate the
heat flux at the nanoscale. For a static thermal bias the most explored
control mechanism so far has been by means of tuning the structural asymmetry
and the degree of anharmonicity in tailored one-dimensional lattice
structures. The ensuing dependence on temperature of the power spectra of
dissimilar segments results in the phenomenon of thermal rectification, i.e.
asymmetrical heat flow, with substantial progress being achieved in the last
two decades Roberts11 . This phenomenon offers improved thermal management at
the nanoscale and its success can be gauged by the fact that, only shortly
after the first theoretical models of the heat rectifier Terraneo02  and
thermal memory Wang08  were developed, successful experimental realizations
were reported Chang06a ; Xie11 .

In order to obtain an even more flexible control of heat energy, comparable
with the richness available for electronics, one may utilize temporal
modulation that directs heat from one part of the device to another or to an
external reservoir by means of an applied external work. Models on the
mechanism of such a nanoscale heat pump have been proposed in systems, mostly
coupled anharmonic lattices, where there is no net thermal bias between the
two reservoirs Zhan09 ; Nianbei09 ; Ren10  and where the pump works against
the imposed static thermal gradient Broek06 ; Bao-quan10 ; Zhang11 .
Furthermore, models based on pumping phenomena have been proposed as moving
barriers in a cavity to pump phonons from a cold reservoir to a hotter one, or
driven two-level systems or molecular junctions in asymmetric contact with
phononic baths characterized by different spectral properties Arrachea2012 ;
Segal06 . In addition, other models employing quantum particle pumps that
differentiate and filter hot and cold particles have been proposed Rey07 .
From the above examples it can be inferred that among the necessary
prerequisites to run such heat machinery are nonlinearity, thermal noise,
unbiased nonequilibrium driving, and a symmetry-breaking mechanism. However,
other than the aforementioned information, externally driven energy transport
remains poorly understood. Furthermore, since its study is far from trivial,
contradictory results have been reported. For example, in Ref. Bao-quan10  it
was claimed that heat pumping appeared in the Frenkel-Kontorova (FK) lattice
under the influence of a periodic driving force, but later it was demonstrated
that such an effect is indeed absent Zhang11 .

In this work we reconsider the one-dimensional model consisting of two
dissimilar FK lattices connected together by a time-modulated harmonic
coupling under the influence of a static thermal bias previously studied in
Ref. Beraha15 . Our results seem to indicate that the interpretation of some
of the results in the aforementioned work may be questionable. More precisely,
the obtained resonant heat transport regime, i.e. maximization of heat flux
for a specific value of the external driving frequency, was explained by a
shift towards lower frequencies of the phonon bands in each segment, induced
by a transition from anharmonic to harmonic behavior as the resonant frequency
is approached. Nevertheless, we have observed that such a phonon band shift
can only be attained by the imposed structural asymmetries of the model. We
have introduced them in a systematic way by means of a scaling of the involved
parameters, first proposed in Ref. Li04a  to study the thermal rectification
properties of this model and extensively employed afterwards Hu06 ; Li12 .
Once the asymmetry is fixed, no shift in the phonon bands is observed
altogether; these are only weakly modified by the external drive. Furthermore,
we have determined that, for the temperature bias employed in Ref. Beraha15
and in the present work, the system is in its harmonic regime. The overlap of
the ensuing phonon bands determines the frequency range wherein thermal
resonance can manifest itself, in close analogy to the way such overlap is at
the origin of thermal rectification in this and similar systems. Thus we
unravel the underlying physical mechanism for such a resonance phenomenon.

This paper is organized as follows: in Sec. II the model system and
methodology are presented. Our results on the dependence of the thermal
resonance on the structural parameters of the model are reported in Sec. III.
The discussion of the results, as well as our conclusions, are presented in
Sec. IV.

##  II The Model

Our system consists of two segments ($L,R$) of nonlinear 1D lattices coupled
together by a harmonic spring with a time-modulated strength $k_{{}_{C}}(t)$,
as shown in Fig. 1. The equations of motion (EOM) for a given oscillator
within each segment can be written, in term of dimensionless variables, as
$\dot{q}_{i}=p_{i}/m_{i}$ and

$\displaystyle\dot{p}_{i}$ $\displaystyle=$ $\displaystyle
k_{{}_{L,R}}(q_{i+1}+q_{i-1}-2q_{i})-{V_{{}_{L,R}}\over 2\pi}\sin(2\pi q_{i})$
(1) $\displaystyle+$
$\displaystyle(\xi_{{}_{1}}-\gamma_{{}_{1}}p_{i})\,\delta_{i1}+(\xi_{{}_{N}}-\gamma_{{}_{N}}p_{i})\,\delta_{i{N}},$
(2)

where $N$ is the system size. $k_{{}_{L,R}}$ and $V_{{}_{L,R}}$ are the
harmonic spring constant and the amplitude of the FK on-site potential in each
segment, respectively; more precisely, in the above equations we employ
$k_{{}_{L}},V_{{}_{L}}$ if $i\in[1,n_{{}_{c}}-1]$ and $k_{{}_{R}},V_{{}_{R}}$
if $i\in[n_{{}_{c}}+2,N]$. In order to reduce the number of adjustable
parameters, we set $V_{{}_{R}}=\lambda V_{{}_{L}}$ and $k_{{}_{R}}=\lambda
k_{{}_{L}}$. Here we consider only the commensurate case where the on-site
potential assumes the same spatial periodicity as the lattice constant.
$\\{m_{i},q_{i},p_{i}\\}_{i=1}^{N}$ are the dimensionless mass, displacement,
and momentum of the $i$th oscillator; see the Appendix of Ref. Li12  for a
detailed procedure on how to construct such dimensionless variables. Fixed
boundary conditions are assumed ($q_{{}_{0}}=q_{{}_{{N}+1}}=0$). Henceforth we
will consider a homogeneous system, i.e., $m_{i}=1\,\,\forall\,i$. The
Gaussian white noise $\xi_{{}_{1,N}}$ has zero mean and correlation
$\langle\xi_{{}_{1,N}}(t)\xi_{{}_{1,N}}(t^{\prime})\rangle=2\gamma_{{}_{1,N}}k_{{}_{B}}T_{{}_{1,N}}m_{i}(\delta_{1i}+\delta_{{N}i})\delta(t-t^{\prime})$,
with $\gamma_{{}_{1,N}}$ (taken as $0.5$ in all computations hereafter
reported) being the coupling strength between the system and the left and
right thermal reservoirs operating at temperatures $T_{{}_{L}}=0.15$ and
$T_{{}_{R}}=0.05$, respectively; the system thus operates at a constant
average temperature value of $T_{{}_{0}}\equiv(T_{{}_{L}}+T_{{}_{R}})/2=0.1$.
The EOM for the last oscillator in the first segment ($L$) and the first one
in the second ($R$) are given by

$\displaystyle\dot{p}_{n_{{}_{c}}}$ $\displaystyle=$ $\displaystyle
k_{{}_{L}}(q_{{n_{{}_{c}}}-1}-q_{n_{{}_{c}}})+k_{{}_{C}}(t)(q_{{n_{{}_{c}}}+1}-q_{n_{{}_{c}}})$
(3) $\displaystyle-$ $\displaystyle{V_{{}_{L}}\over 2\pi}\sin(2\pi
q_{n_{{}_{c}}}),$ (4) $\displaystyle\dot{p}_{{n_{{}_{c}}}+1}$ $\displaystyle=$
$\displaystyle
k_{{}_{R}}(q_{{n_{{}_{c}}}+2}-q_{{n_{{}_{c}}}+1})+k_{{}_{C}}(t)(q_{n_{{}_{c}}}-q_{{n_{{}_{c}}}+1})$
(5) $\displaystyle-$ $\displaystyle{V_{{}_{R}}\over 2\pi}\sin(2\pi
q_{{n_{{}_{c}}}+1}),$ (6)

with $k_{{}_{C}}(t)=k_{{}_{0}}(1+\sin\omega t)$ being the time-modulated
amplitude of the harmonic coupling, which is an external driving with
frequency $\omega$. The aforementioned EOM were integrated with a stochastic
velocity-Verlet integrator with a time step of $\Delta t=0.005$ for a
production time interval of $2\times 10^{7}$ time units after a transient time
of $10^{8}$ time units.

Figure 1: (Color online) Sketch of our model system composed of two dissimilar
FK lattices connected by a time-modulated harmonic interaction and in contact
with two thermal reservoirs.

Once the nonequilibrium stationary state is attained, the local heat flux is
computed as ${J}_{i}=k_{{}_{L,R}}\langle\dot{q}_{i}(q_{i}-q_{i-1})\rangle$,
with $k_{{}_{L}}$ if $i\in[2,n_{{}_{c}}-1]$ and $k_{{}_{R}}$ if
$i\in[n_{{}_{c}}+1,N]$, and the local temperature as $T_{i}=\langle
p_{i}^{2}/m_{i}\rangle$; in both instances $\langle\cdots\rangle$ indicates
time average. In the stationary state the heat flux in each segment becomes
independent of the site, and, in order to improve the statistical precision of
our results, the mean heat flux $J_{{}_{L,R}}$ on each side of the lattice is
calculated as the algebraic average of ${J}_{i}$ over the number of
unthermostatted oscillators in each segment. Now, the rate of work $\dot{W}$
done by the external driving in the contact at $n_{{}_{c}}$ is dissipated into
the reservoirs, implying that

$\displaystyle\dot{W}=J_{{}_{L}}+J_{{}_{R}},$ (7)

where $J_{{}_{L,R}}$ are defined as positive when the heat flows into the
reservoirs.

##  III Results

In Fig. 2 we present the results of the dependence of heat fluxes
$J_{{}_{L}}$, $J_{{}_{R}}$, and the average $J$ as a function of the driving
frequency $\omega$ with $V_{{}_{L}}=5$, $V_{{}_{R}}=1$ $k_{{}_{L}}=1$,
$k_{{}_{R}}=0.2$, $k_{{}_{0}}=0.05$, and $n_{{}_{c}}=N/2$ for a lattice with
$N=32$ oscillators. Our results in Fig. 2(a) are qualitatively similar to
those reported in Ref. Beraha15 . The observed differences are clearly a
result of the way in which we chose to implement the structural asymmetry into
the system, which is that of Ref. Li04a  as previously mentioned. In the
adiabatic driving limit $\omega\rightarrow 0$ the heat flows from oscillator
$i=1$ to $i=N$, with $J_{{}_{R}}=-J_{{}_{L}}>0$, and thus the averaged net
power released to the system is zero. In the opposite limit
$\omega\rightarrow\infty$ the coupling oscillates very fast and converges to a
time average constant value $k_{{}_{0}}$ as if there is no driving. Thus the
relevant phenomenology occurs at intermediate $\omega$ values. First, for
$0.1<\omega\lesssim 0.26$ there is a small net power contribution released to
the system, and now $J_{{}_{R}}-|J_{{}_{L}}|=P>0$; thus, although heat still
flows from the hot to the cold reservoir, $J_{{}_{L}}$ and $J_{{}_{R}}$ have
different magnitudes. Then, within $0.26\lesssim\omega\lesssim 2$, the power
released in the contact region is dissipated into the reservoirs, since
$J_{{}_{L}}>0$. This phenomenon, due to the resonant interaction of the
external drive with the system’s intrinsic frequencies, is known as thermal
resonance Bao-quan10 ; Zhang11 ; Beraha15 . Its main feature is the
maximization of the heat flux $J_{\mathrm{m}}$ at a specific frequency,
$\omega_{\mathrm{m}}=0.6$ in this particular instance.

Figure 2: (Color online) Heat flux vs driving frequency $\omega$. The energy
currents through the left and right segments are $J_{{}_{L}}$ and
$J_{{}_{R}}$, respectively; $J\equiv(J_{{}_{L}}+J_{{}_{R}})/2$ is the average
current with $V_{{}_{L}}=5$, $k_{{}_{L}}=1$, $\lambda=0.2$, $k_{{}_{0}}=0.05$,
$N=32$, and $n_{{}_{c}}=N/2$. Maximum heat flux is at
$\omega_{\mathrm{m}}=0.6$. Continuous lines are a guide to the eye.

Next we will explore the relative contribution of low- and high-frequency
phonons to the heat released into the system through the contact region. In
Fig. 3 we plot the phonon spectra
$|\tau^{-1}\\!\\!\int_{{}_{0}}^{\tau}\\!\\!dt\dot{q}_{i}(t)\exp(-\mathrm{i}\Omega
t)|^{2}$ of the interface oscillators at the left and right sides of the
contact for both the adiabatic driving limit, $\omega=0.001$, and in the
thermal resonance regime, $\omega_{\mathrm{m}}=0.6$, which corresponds to the
maximum heat flux $J_{\mathrm{m}}$ of Fig. 2. It is evident that the spectra
in the latter regime have almost twice the spectral power as those in the
former, which results in a larger overlap in the low-frequency range of the
thermal resonance regime depicted in Fig. 3(b). The energy transport into the
reservoirs goes through the phonon channels determined by the imposed thermal
bias since the structure, unlike the magnitude, of the spectra remains largely
unchanged. The results in our Fig. 3 can be understood if we recall that there
is a critical temperature $T_{{}_{\mathrm{cr}}}\approx V/(2\pi)^{2}$ above
which the kinetic energy is large enough to overcome the on-site potential
barrier, hence the contribution of the on-site potential can be neglected.
Li04a . In our particular case we have
$T_{{}_{\mathrm{cr}}}^{(L)}=0.13<T_{{}_{L}}$ for $V_{{}_{L}}=5$ and
$T_{{}_{\mathrm{cr}}}^{(R)}=0.025<T_{{}_{R}}$ for $V_{{}_{R}}=1$. Thus both
sides of the system are in a temperature regime, well above their respective
$T_{{}_{\mathrm{cr}}}$ values, wherein they behave as harmonic lattices with a
phonon band of $0<\Omega<\sqrt{4k_{{}_{L,R}}}$ composed mainly of
noninteracting phonons which gives $0<\Omega/2\pi\lesssim 0.32$ for the left
oscillator and $0<\Omega/2\pi\lesssim 0.14$ for the right one Li04a . In the
weak-coupling limit ($k_{{}_{0}}\ll 1$) thermal resonance can occur only for
frequency values in the overlapping region of these phonon bands, which
implies that $\omega_{\mathrm{m}}/2\pi<0.14$; our result
$\omega_{\mathrm{m}}/2\pi\sim 0.09$ is in good agreement with the above
estimate. Therefore, the net energy flow from the external source into the
thermal reservoirs is accomplished by the external driving due to its
interaction with —and ensuing alteration of— the phonon bands activated by the
thermal bias imposed at the boundaries.

Figure 3: (Color online) (a) Power spectra of the two oscillators in the left
(red) and right (blue) sides of the contact for a driving frequency of
$\omega=0.001$. (b) Same as (b) but for $\omega_{\mathrm{m}}=0.6$. Vertical
dashed and dot-dashed lines correspond to the cut-off frequencies of the left
and right phonon bands respectively. Vertical solid line denotes the
$\omega_{\mathrm{m}}/2\pi$ value. Same $V_{{}_{L}}$, $k_{{}_{L}}$, $\lambda$,
$k_{{}_{0}}$, $N$, and $n_{{}_{c}}$ values as in Fig. 2.

In Fig. 4 we plot the corresponding temperature profiles for selected $\omega$
values reported in Fig. 2. It is clear that, near and at
$\omega_{\mathrm{m}}$, there is a change in the sign of the slope
corresponding to the left side of the system, which signals that there is an
energy flux into the hot reservoir. The slope of the temperature profile in
the right side is also increased with respect to its value in the adiabatic
regime because of the additional energy flux afforded by the external driving.
Due to the temperature jump at the right boundary the slope of the temperature
profile is higher than that in the left side, which correlates well with the
fact that $J_{{}_{R}}>J_{{}_{L}}$ for $\omega=\omega_{\mathrm{m}}$. Now,
contrary to Ref. Beraha15 , all of the temperature profiles reported in Fig. 4
exhibit a discontinuity at the interface. This result is consistent with the
existence of the contact between both lattices, as shown in numerous studies
Li04a ; Hu05 ; Luo07 . We corroborated that this temperature jump is still
present when simulations were performed with the exact values of the employed
parameter set reported in the aforementioned reference and different
production time intervals, as can be readily seen in the Fig. 4 inset.

Figure 4: (Color online) Temperature profiles for selected $\omega$ values,
from low to high ones, increasing from bottom to top. Same $V_{{}_{L}}$,
$k_{{}_{L}}$, $\lambda$, $k_{{}_{0}}$, $N$, and $n_{{}_{c}}$ values as in Fig.
2. Inset is the temperature profile for the same values of the parameter set
as in Ref. Beraha15  with a resonant frequency of $\omega_{\mathrm{m}}=0.3$
for production time intervals of $5\times 10^{6}$ (circle), $2\times 10^{7}$
(triangle), and $10^{8}$ (square) time units.

As already explained in Sec. II, both the magnitude of the elastic constant
and the strength of the on-site potential can be simultaneously controlled by
means of the $\lambda$ parameter. Therefore, it is sufficient to study the
properties of the considered system as a function of the aforementioned
parameter. Figure 5 reports the dependence of $J_{\mathrm{m}}$ versus
$\lambda$. This figure clearly shows that, for $\lambda<1$, we have
$J_{{}_{R}}>J_{{}_{L}}$, i.e., more heat flows into the cold reservoir. Next,
as $\lambda$ increases, and thus the asymmetry in the system decreases,
$J_{{}_{R}}$ decreases. For $\lambda>1$ values there is an almost complete
suppression of the heat flux in the right side of the lattice. In this
situation the hight of the valley and the harmonic constant stiffness become
large enough to confine the oscillators in the potential valley, thus
preventing any significant heat flux altogether into the right segment; only a
steady heat flow into the left (hot) reservoir through the corresponding
segment is present. In the inset we plot $\omega_{\mathrm{m}}$ versus
$\lambda$, and it can be readily observed that, as $\lambda$ increases, the
resonant frequency has a sharp increase in its value when $\lambda\sim 0.6$.

Figure 5: (Color online) Heat flux $J_{\mathrm{m}}$ vs scaling parameter
$\lambda$. Same $V_{{}_{L}}$, $k_{{}_{L}}$, $k_{{}_{0}}$, $N$, and
$n_{{}_{c}}$ values as in Fig. 2. Inset is $\omega_{\mathrm{m}}$ vs $\lambda$.
Continuous lines are a guide to the eye.

For the instances in Fig. 5 wherein $\lambda<0.6$ we have
$\omega_{\mathrm{m}}<1$, as can be seen in the inset, with a dependence of the
heat flux on the driving frequency very similar to that already depicted in
Fig. 2. In Fig. 6 we display the data of the heat flux as a function of the
driving frequency $\omega$ for the particular $\lambda=0.8$ value. It is to be
noted that a peak associated with a $\omega<1$ frequency is still present at
$\omega_{\mathrm{l}}\sim 0.8$, but is now accompanied by a second peak, that
now becomes the absolute maximum and thus the resonant frequency, at
$\omega_{\mathrm{m}}\sim 1.6$, as can be readily noticed, with
$J_{{}_{R}}>J_{{}_{L}}$.

Figure 6: (Color online) (a) Heat flux vs driving frequency $\omega$ for
$\lambda=0.8$. Same $V_{{}_{L}}$, $k_{{}_{L}}$, $k_{{}_{0}}$, $N$, and
$n_{{}_{c}}$ as in Fig. 2. The two peaks of the heat flux correspond to
frequencies $\omega_{\mathrm{l}}=0.8$ and $\omega_{\mathrm{m}}=1.6$.
Continuous lines are a guide to the eye.

The power spectra of the oscillators to the left and right of the contact for
the $\lambda=0.8$ value of Fig. 5 are displayed in Fig. 7(a). The increased
structural symmetry of the system is reflected in the qualitative similarity
of the two spectra, and the ensuing overlap over the entire frequency range
accounts for the widening of the range wherein subresonant behavior can be
observed. The left peak of the $J$ vs $\omega$ plot in Fig. 6 can be
associated with the leftmost ones of the displayed spectra since they are
located at $\Omega/2\pi\sim 0.13$, a value almost identical to
$\omega_{\mathrm{l}}/2\pi$. However, the value for which thermal resonance
appears, $\omega_{\mathrm{m}}/2\pi\sim 0.254$, coincides with a spike-like
value in both spectra, as can be readily appreciated. Just as in the instance
depicted in Fig. 3, $\omega_{\mathrm{m}}$ lies in the frequency range wherein
the phonon bands overlap. These are given by $0<\Omega/2\pi\lesssim 0.32$ for
the left side and $0<\Omega/2\pi\lesssim 0.28$ for the right one; thus
$\omega_{\mathrm{m}}/2\pi<0.28$, as is indeed observed. The $\lambda=2$
instance is displayed in Fig. 7(b). As previously mentioned, now $V_{{}_{R}}$
becomes relevant, and thus so does the influence of the anharmonic FK
potential. In this case the lower bound of the phonon band is raised by
$\sqrt{V}_{{}_{R}}$ and the phonon band is shifted to
$\sqrt{V_{{}_{R}}}<\Omega<\sqrt{V_{{}_{R}}+4k_{{}_{R}}}$ Li04a . For the
considered conditions the latter is $0.5<\Omega/2\pi<0.68$, which has no
possible overlap with the left phonon band $0<\Omega/2\pi\lesssim 0.32$.
Therefore it is clear that the only open channel available for heat carrying
phonons in the low-frequency region is now afforded by the latter phonon band.
And indeed the resonant frequency $\omega_{\mathrm{m}}/2\pi\sim 0.127$ lies
within this last frequency range. Thus the heat flux into the left (hot)
reservoir is dominated by low-frequency acoustic phonons. Finally,
$\omega_{\mathrm{m}}$ coincides approximately with the location of the highest
(leftmost) peak of the spectrum and exactly with a small spike-like
perturbation at that same frequency value.

Figure 7: (Color online) (a) Power spectra of the two oscillators in the left
(red) and right (blue) sides of the contact for a $\lambda=0.8$ value. (b)
Same as (b), but for $\lambda=2$. Vertical dashed and dot-dashed lines
correspond to the cutoff frequencies of the left and right phonon bands,
respectively. Vertical solid line denotes the corresponding
$\omega_{\mathrm{m}}/2\pi$ values in each case. Same $V_{{}_{L}}$,
$k_{{}_{L}}$, $k_{{}_{0}}$, $N$, and $n_{{}_{c}}$ values as in Fig. 2.

The temperature profiles for various driving frequency values corresponding to
the $\lambda=0.8$ case depicted in Fig. 6 are plotted in Fig. 8, and an
asymmetry between the two sides of the system is evident, as expected. Now for
the $\omega$ values considered the slopes in the left side are very similar,
indicating that the magnitude of the heat flux towards the left reservoir is
almost the same in all cases, including $\omega=0.8$, i.e. the left most peak.
However, in the right side the slope corresponding to $\omega=1.6$, i.e., the
rightmost peak and absolute maximum, is clearly greater than the other ones,
indicating a heat flux towards the right (colder) reservoir of larger
magnitude than that in the left side. We also notice that the jump in the part
of the temperature profile in contact with the cold reservoir that was
observed for the corresponding $\omega_{\mathrm{m}}$ case in Fig. 4 is now
absent. This effect is most certainly a consequence of the increased symmetry
of both segments, since it is well known that the temperature profile of the
homogeneous FK lattice has no discontinuities in the boundaries Hu98 .

Figure 8: (Color online) Temperature profiles for selected $\omega$ values,
from low to high ones, increasing from bottom to top in the case of
$\lambda=0.8$. Same $V_{{}_{L}}$, $k_{{}_{L}}$, $k_{{}_{0}}$, $N$, and
$n_{{}_{c}}$ values as in Fig. 2.

The interface elastic constant $k_{{}_{0}}$ is a very important parameter as
it plays the role of coupling the two lattices. By adjusting this parameter
one can control the heat flow through the system, as previous work has shown
Hu06 . Indeed, once its value is fixed, then the smaller the coupling is, the
smaller the heat current is through the system. In Fig. 9 we present the
variation of $J_{\mathrm{m}}$ as a function of $k_{{}_{0}}$. It is clear that
the average value of $J_{\mathrm{m}}$ presents a monotonic increase as the
magnitude of the elastic constant grows. However, for $k_{{}_{0}}<0.4$ values
we have $J_{{}_{L}}<J_{{}_{R}}$, i.e., a higher energy flux towards the colder
reservoir. On the other hand, if $k_{{}_{0}}>0.4$ the growth of $J_{{}_{R}}$
diminishes, whereas $J_{{}_{L}}$ keeps increasing its magnitude.

Figure 9: (Color online) Heat flux $J_{\mathrm{m}}$ vs $k_{{}_{0}}$. Same
$V_{{}_{L}}$, $k_{{}_{L}}$, $\lambda$, $N$, and $n_{{}_{c}}$ values as in Fig.
2. Continuous lines are a guide to the eye.

The change in relative magnitude of $J_{{}_{L}}$ and $J_{{}_{R}}$ as
$k_{{}_{0}}$ increases can be explained by examining the power spectra
corresponding to low and high values of the interface elastic constant that
are plotted in Fig. 10; in both instances the left and right phonon bands are
$0<\Omega/2\pi\lesssim 0.32$ and $0<\Omega/2\pi\lesssim 0.14$, respectively.
For $k_{{}_{0}}=0.4$ more power is available in the right spectrum at high
frequencies within the region wherein the phonon bands overlap. Thus the
external driving interacts more closely with the right segment, resulting in a
well defined thermal resonance with $\omega_{\mathrm{m}}/2\pi\sim 0.15$ and,
as a consequence, $J_{{}_{R}}>J_{{}_{L}}$. Next, by examining the spectra for
the $k_{{}_{0}}=0.8$ case depicted in Fig. 10(b) it can be observed that they
become entangled and begin to form a whole. This behavior can be explained by
noticing that, for this high $k_{{}_{0}}$ value, both segments interact more
strongly and act as a single system. Therefore a single phonon band, i.e.
$0<\Omega/2\pi\lesssim 0.32$, determines the frequency values of the
noninteracting phonons within the two segments. Thus the ensuing resonant
frequency $\omega_{\mathrm{m}}/2\pi\sim 0.22$ lies within the aforementioned
phonon band and in a frequency range wherein the overlap of the two spectra is
largest. Now, as these can still be distinguished from one another, it is
clear that more power is provided by the left spectrum. This in turn entails
$J_{{}_{L}}>J_{{}_{R}}$, in agreement with the results of Fig. 5.

Figure 10: (Color online) (a) Power spectra of the two oscillators in the left
(red) and right (blue) sides of the contact for a value of the interface
elastic constant $k_{{}_{0}}=0.4$. (b) Same as (b) but for $k_{{}_{0}}=0.8$.
Vertical dashed and dot-dashed lines correspond to the cut-off frequencies of
the left and right phonon bands, respectively. Vertical solid line denotes the
corresponding $\omega_{\mathrm{m}}/2\pi$ values in each case. Same
$V_{{}_{L}}$, $k_{{}_{L}}$, $\lambda$, $N$, and $n_{{}_{c}}$ values as in Fig.
2.

Finally, we would like to discuss the possibility of an experimental
realization of this device. For typical atom lattices, room temperature,
measured in kelvins, corresponds to a dimensionless temperature
$T_{{}_{0}}\in[0.1,1]$ Li12 ; thus the employed $T_{{}_{0}}$ value is within
the adequate range considering future technological applications. For a
lattice constant of 1 Å, a lattice of the size herein used is about 32
nanometers long, a size scale within the reach of current technology in order
to be built. Furthermore, at molecular levels a modulation of the coupling
between two molecules can be achieved experimentally in molecular junctions
by, for example, harmonically varying the distance among them, therefore
modulating the coupling between the molecules.

##  IV Concluding remarks

To summarize: we have studied energy transport control in a one-dimensional
segmented system composed of two FK lattices connected by a time modulated
coupling. This model affords a convenient way to study dynamical control of
heat transport and obtain results that might be relevant for nanoscale
devices. Our analysis reveals that, as far as the resonance property is
concerned, there is much similarity between the FK and the harmonic models and
is consistent with previous findings on the subject Zhang11 . By an
appropriate scaling we have reduced the number of parameters involved in the
description of the system and determined that the phonon heat transport
properties in the harmonic limit of the undriven lattice are crucial to
control the heat flux through the segments when thermal resonance is present.
If values of the amplitude of the on-site potential and harmonic constant in
the left segment are greater than those in the right one, energy flow into the
colder reservoir is greater than the flow into the hot one, whereas the
reverse is true in the opposite case. Also, heat flow into the colder
reservoir is higher than into the hotter one for low values of the harmonic
coupling constant between segments; the opposite being true for higher values
of that same constant. We expect that our results hold for larger system sizes
since the herein studied system is essentially the same as that in Ref.
Beraha15 . Therefore, the shift to lower frequencies of the resonant frequency
value reported in the latter work can be expected as the system size
increases.

Previously the thermal transport properties of a harmonic lattice system
consisting of two semi-infinite leads at different temperature and connected
by a time-modulated coupling were studied by means of the nonequilibrium
Green’s function formalism Cuansing10 . There is a net energy flow out of the
warmer lead, but for the colder lead energy flow direction depends on the
values of both the driving frequency and temperature. Later, a model similar
to that herein employed and in Ref. Beraha15 , but with the thermal reservoirs
modeled as infinite harmonic lattices and the on-site potential being harmonic
as well Beraha16 , was proposed and studied analytically by the same technique
as in the last reference. One of the transport regimes studied corresponds to
a heat pump against the imposed thermal gradient. It would be interesting to
explore the possibility of obtaining the same effect with the anharmonic model
herein employed.

###### Acknowledgements.

The author thanks CONACYT, México for financial support, Brandon Armando
Matínez Torres for his help in rendering the figures, and Juan M. López and
Susana Sanchez Merodio for insightful comments and discussions.

## References

  * (1) M. Tsutsui, M. Taniguchi, and T. Kwai, Nano Lett. 8, 3293 (2008).
  * (2) Z. Juang et al., Nature Nanotech. 2, 698 (2007).
  * (3) M. Tsutsui and Y.-C. Chen, Sci. Rep. 9, 18677 (2019).
  * (4) Y. T. Yang et al., Nano Lett. 6, 583 (2006).
  * (5) D. Rugar, R. Budakian, H. J. Mamin, and B. W. Chui, Nature 430, 329 (2004).
  * (6) A. N. Cleland and M. L. Roukes, Nature 392, 160 (1998).
  * (7) S. Volz et al., Eur. Phys. J. B 89, 15 (2016).
  * (8) N. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi, and B. Li, Rev. Mod. Phys. 84, 1045 (2012).
  * (9) N. A. Roberts and D. G. Walker, Int. J. Therm. Sci. 50, 648 (2011).
  * (10) M. Terraneo, M. Peyrard, and G. Casati, Phys. Rev. Lett. 88, 094302 (2002).
  * (11) L. Wang and B. Li, Phys. Rev. Lett. 101, 267203 (2008).
  * (12) C. W. Chang, D. Okawa, A. Majumdar, and A. Zettl, Science 314, 1121 (2006).
  * (13) R. Xie et al., Adv. Funct. Mater. 21, 1602 (2011).
  * (14) F. Zhan, N. Li, S. Kohler, and P. Hänggi, Phys. Rev. E 80, 061115 (2009).
  * (15) N. Li, F. Zhan, P. Hänggi, and B. Li, Phys. Rev. E 80, 011125 (2009).
  * (16) J. Ren and B. Li, Phys. Rev. E 81, 021111 (2010).
  * (17) C. Van den Broeck and R. Kawai, Phys. Rev. Lett. 96, 210601 (2006).
  * (18) B. Q. Ai, D. He, and B. Hu, Phys. Rev. E 81, 031124 (2010).
  * (19) S. Zhang, J. Ren, and B. Li, Phys. Rev. E 84, 031122 (2011).
  * (20) L. Arrachea, E. R. Mucciolo, C. Chamon, and R. B. Capaz, Phys. Rev. B 86, 125424 (2012).
  * (21) D. Segal and A. Nitzan, Phys. Rev. E 73, 026109 (2006).
  * (22) M. Rey, M. Strass, S. Kohler, P. Hänggi, and F. Sols, Phys. Rev. B 76, 085337 (2007).
  * (23) N. Beraha, A. Soba, R. Barreto, and M. F. Carusela, Physica A 433, 9 (2015).
  * (24) B. Li, L. Wang, and G. Casati, Phys. Rev. Lett. 93, 184301 (2004).
  * (25) B. Hu, L. Yang, and Y. Zhang, Phys. Rev. Lett. 97, 124302 (2006).
  * (26) B. Hu and L. Yang, Chaos 15, 015119 (2005).
  * (27) C. Luo and J. Dong, Phys. Rev. B 75, 132301 (2007).
  * (28) B. Hu, B. Li, and H. Zhao, Phys. Rev. E 57, 2992 (1998).
  * (29) E. C. Cuansing and J.-S. Wang, Phys. Rev. E 82, 021116 (2010).
  * (30) N. Beraha, A. Soba, and M. F. Carusela, Physica A 464, 64 (2016).