# condLSTM-Q: A novel deep learning model for predicting Covid-19 mortality in
fine geographical Scale

HyeongChan Jo  Division of Biology and Biological Engineering, Caltech  These
authors contributed equally to this work. Juhyun Kim  The Division of Physics,
Mathematics and Astronomy, Caltech  These authors contributed equally to this
work. Tzu-Chen Huang  Walter Burke Institute for Theoretical Physics, Caltech
Yu-Li Ni  Division of Biology and Biological Engineering, Caltech
Corresponding author<EMAIL_ADDRESS>M.D  Division of Biology and Biological
Engineering, Caltech  Corresponding author<EMAIL_ADDRESS>
###### Abstract

Predictive models with a focus on different spatial-temporal scales benefit
governments and healthcare systems to combat the COVID-19 pandemic. Here we
present the conditional Long Short-Term Memory networks with Quantile output
(condLSTM-Q), a well-performing model for making quantile predictions on
COVID-19 death tolls at the county level with a two-week forecast window. This
fine geographical scale is a rare but useful feature in publicly available
predictive models, which would especially benefit state-level officials to
coordinate resources within the state. The quantile predictions from
condLSTM-Q inform people about the distribution of the predicted death tolls,
allowing better evaluation of possible trajectories of the severity. Given the
scalability and generalizability of neural network models, this model could
incorporate additional data sources with ease, and could be further developed
to generate other useful predictions such as new cases or hospitalizations
intuitively.

##  1 Introduction

The coronavirus disease 2019 (COVID-19) pandemic has taken more than 1.1
million lives worldwide and more than 224,000 lives in the United States as of
October 2020. Predicting the trend of the pandemic in terms of deaths and
positive cases precisely has been crucial, as it allows the governments and
healthcare systems to effectively distribute and prioritize resource
allocations [1].

Here we present the conditional Long Short-Term Memory networks with Quantile
output (condLSTM-Q) model, a novel deep learning model predicting the spatial
and temporal distribution of COVID-19 outbreak. The condLSTM-Q was developed
during the Caltech COVID-19 Initiative [2], a campaign aimed to develop and
explore novel models to complement the classic epidemiological models. The
optimization goal for participating models was to predict daily death tolls
attributed to COVID-19 in each county across the United States with a two-week
future prediction window. In addition, the models were required to provide
estimations of 10-quantiles (0.1 to 0.9 quantile with intervals of 0.1) as the
outputs, as quantile prediction is useful when forecasting the extremes. [3,
4]. The county-level spatial resolution was selected to inform state officials
and local governments so that they can make timely decisions on resource
allocations and strengthen their pandemic preparedness. The condLSTM-Q was
among one of the best performers by the end of the campaign in early June. To
test the robustness of this model, the model was then deposited without
architectural modifications since mid-June, 2020, and has continued to be
trained and output predictions using updated data. Since its deposit, the
model’s predictions were comparable to the well-known, publicly available
predictions by the Institute for Health Metrics and Evaluation (IHME) [5]
throughout May-October when this manuscript was written, showcasing the
robustness of condLSTM-Q.

The foundation of the condLSTM-Q was based on the Long Short-Term Memory
networks (LSTM), a common neural network model well-suited for time series
predictions [6]. There are previous works of using LSTM and other neural
networks for COVID-19 forecasting, such as risk assessment of countries [7],
and predicting national COVID-19 mortality rates [8]. Nevertheless, to our
knowledge no county-level, quantile prediction model existed when our
prototype was developed. Lack of county-based predictions could hugely impact
the prevention and control of COVID-19 pandemic at local governments, which
plays pivotal roles in the pandemic preparedness in the United States. Making
quantile predictions rather than single-point predictions also gives a sense
of the distribution of the predicted values, which can help governments to
understand the situation better. Thus we decided to explore if the LSTM-based
model could provide good quantile predictions at the county level. We improved
the classical LSTM model by adapting and building on a "conditional" LSTM
architecture which takes in static data more naturally along with time series
data [9, 10]. In addition, we utilized the flexibility of neural networks to
output a distribution of 10-quantiles that was required by the initiative.

Given the robust performance of condLSTM-Q and the fact that we have not
exploited all possible variations of the model, this architecture from our
pilot experiment is worthy of further exploration. Interdisciplinary
collaboration between experts in machine learning and epidemiology would
greatly facilitate building better variants that will provide both longer
prediction windows and clearer interpretability.

##  2 Methods

###  2.1 Data Sources and Data Preprocessing

Our model includes a variety of data that is expected to have correlations
with the death toll. The original source of each data is given in Section 5.
The data used in our model fall into one of the following seven categories:
(1) COVID-19 mortality and confirmed cases provided by New York Times; (2)
demographics and local health resources such as age composition, mortality
rate by diseases, and the number of hospitals, provided by the Yu group at
University of California, Berkeley; (3) county-wise gross domestic product
(GDP) from the Bureau of Economic Analysis; (4) population density and
geographical data from 2010 Census; (5) mobility changes in response to
COVID-19 provided by Descartes Labs; (6) policy actions in response to
COVID-19 such as state of emergency declaration, safe-at-home order, and
business closure, from Covidvis team at University of California, Berkeley
[11] and the U.S. Department of Health and Human Services; (7) and the U.S.
pneumonia and influenza mortality report from the National Center for Health
Statistics Mortality Surveillance System in Centers for Disease Control and
Prevention (CDC). All the data were county-level, and the latest start date of
the time series data was March 1, 2020.

Preprocessing specifics for each dataset are listed in the following:

  * •

Mobility data from the Descartes Labs was provided for 2,721 counties out of
3,114 counties in total. The missing 393 counties’ mobility data were filled
in with their corresponding state-level data. Missing dates in mobility data
of Descartes Labs were interpolated with the data on the closest existing date
of the same day in a week to reflect the weekly pattern inherent in the data.
As of November 8, 2020, 14,049 values were missing, out of total 688,413 data
points (about 2%). They were filled in using either the aforementioned same-
day interpolation or spline interpolation, and the data from the first and the
last day of recording was repeated to fill the missing data before and after
the existing data, if necessary. Spline interpolation has an advantage over
the same-day interpolation in handling a large number of missing data over a
long period of time, whereas the latter can capture the weekly pattern that
cannot be maintained in the former method.

  * •

Seasonality feature was extracted from the U.S. pneumonia and influenza data,
under the premise that the COVID-19 will follow a seasonality of virus that
becomes more dormant during the summer and severe in the winter. More
specifically, the multiplicative seasonality was extracted from the state-wide
pneumonia and influenza mortality rate during flu seasons from 2013 to 2020
provided by CDC [12].

  * •

Among 64 features of the demographics and local health data, 43 features
including the population estimate and the mortality rate of various underlying
diseases were selected as the static features.

  * •

For the policy actions features, declaration of the state of emergency, safe-
at-home action, and the closure of inessential business were selected as the
static features.

All numerical values were standardized before the model training step.

###  2.2 Implementation of condLSTM-Q

The classical LSTM which the condLSTM-Q’s backbone was based on is shown in
Fig 1A. LSTM is a type of Recurrent Neural Network (RNN) that can be trained
to learn the mapping from the time series features to the time series of
interest. Compared to previous architectures of RNN, LSTM has the advantage of
remembering long-term dependencies by passing the information through its cell
state and hidden states recurrently. However, the classical LSTM architecture
cannot take into account non-time series data (hereinafter referred to as
categorical data) in a natural way; such features have to be stacked into the
same dimension of the time series data during data processing as shown in Fig
1 A, in order to take both types of data into the input stream. This way of
data processing undermines the optimal performance of LSTM.

One approach to overcome this limitation is to initialize the initial states
of the LSTM units in response to the categorical data. Intuitively, the
categorical data should provide a priori information to the prediction, rather
than real-time information. Thus, a "conditional" module can be added to get
these categorical data and feed in "priors" as the hidden states, as
implemented in [9], in contrast to the usual LSTM where the hidden states are
initialized to zeros or random noise [13].

This is the way we treat the categorical data in our model, conditional LSTM
with Quantile output model (condLSTM-Q; Fig 1 B). In this model, categorical
inputs are passed through a fully connected layer into the model as the hidden
states to the initial step of the LSTM layer. Through this design, we could
input 50 categorical data such as Income, Age and Gender Distribution,
Medicare coverage, and Population in their original forms of one scalar per
feature into our model. After the initialization step, the cell states and
hidden states in LSTM were then updated as in the classical LSTM, based on 8
time series features including mobility of the population, new cases per day,
and new deaths per day.

Another feature of condLSTM-Q is that it generates predictions on 10-quantiles
(“Q"). To get the sense of the distribution of the predicted values, which
help us understand the situation and the quality of the prediction better (a
forecast with higher variance results in huge uncertainty in whatever
derivatives from the forecast), our model forecasts the death counts in each
quantile from 0.1 to 0.9 with 0.1 increments by a variant of the multi-output
LSTM. In this model, the outputs from the LSTM layer are fed into 9 parallel
dense layers to generate a vector of dimension 9 for each day in the
prediction, where each component of the output vector is a forecast of each
quantile. A standard measure for the accuracy of such a quantile forecasting
is the pinball loss defined as

$\displaystyle loss=\max(q\cdot\epsilon,(q-1)\cdot\epsilon),$ (1)

where $0<q<1$ is the quantile and $\epsilon=y-\hat{y}$ is the difference
between the true target value $y$ and the predicted target value $\hat{y}$. As
such, our model is trained to minimize the average pinball loss

$\displaystyle\frac{1}{9}\sum_{i=1}^{9}\mathrm{max}(q_{i}\cdot\boldsymbol{\epsilon_{i}},(q_{i}-1)\cdot\boldsymbol{\epsilon_{i}}),$
(2)

where $q_{i}=i/10,1\leq i\leq 9$ and
$\boldsymbol{\epsilon_{i}}=\boldsymbol{y}-\boldsymbol{\hat{y_{i}}}$ is the
difference between the true target vector and the predicted target vector for
quantile $q_{i}$.

###  2.3 Hyperparameters and the Training

The time series data after preprocessing is an array of shape
$(\texttt{\\#}\textrm{counties},\texttt{\\#}\textrm{dates},\texttt{\\#}\textrm{features})$,
and the categorical data is an array of shape
$(\texttt{\\#}\textrm{counties},\texttt{\\#}\textrm{features})$. For each date
in the time series data, the data is further split into a history window of
size $7$ and a target window of size $14$, so that condLSTM-Q can learn to
predict the mortality rate for the next 14 days based on the history of prior
7 days.

For hyperparameter tuning, we held out the last $21$ days’ worth of data for
validation. As of August 4, 2020, we used 136 days’ worth of data until July
14, 2020, for the training set, and the remaining 21 days’ worth of data for
the validation set. Since the data was from 3,114 counties, and the input and
the target for the model should be 7 and 14 days long, this led to 316,224
samples for training and 3,114 samples for validation.

The optimal hyperparameters selected were 128 units in LSTM, a learning rate
of 0.001, and a dropout rate of 0.2. In this setting, overfitting was observed
after $~{}30$ epochs, so the model was trained with the ADAM optimizer [14]
for 20 epochs to avoid overfitting.

##  3 Result

Figure 1: Classical LSTM vs condLSTM-Q Architecture. Classical LSTM
architecture (Top). The classical LSTM can predict time series better than
traditional RNN by remembering and passing cell states and hidden states over
the timesteps. The hidden states are usually initialized to zero or random
values. Categorical features have to be stacked and replicated to the same
dimension to match the time series features. The classical structure outputs
one time series. The condLSTM-Q architecture (Bottom). In condLSTM-Q,
categorical inputs are passed through a fully connected layer as hidden states
to the initial step ("Conditional" based on the static information) of the
LSTM layer. A dropout layer is also applied to the dense layer to prevent
overfitting to the categorical inputs. The output from the LSTM layer is fed
into 9 parallel dense layers to provide 9 outputs for 10-quantiles ("Q"). Each
of the 9 parallel outputs is given to the pinball loss function during
training, thus the parallel dense layer gets updated to match the
corresponding quantiles by minimizing the summed loss.

###  3.1 Predictions and Performance

The condLSTM-Q provides a 14-day, county-level prediction with 10-quantiles.
To illustrate this, we show predictions for several representative counties,
identified by their Federal Information Processing Standards- or FIPS-based
county code, including Cook County, Los Angeles, New York, Wayne,
Philadelphia, Hennepin, Maricopa and Montgomery from June 12 to June 25, 2020
(See Fig 3). This prediction interval was immediately after we finalized and
froze the model architecture. The model was trained with data up to June 11,
2020, and was agnostic to the data afterward (i.e. post-June 12). The
condLSTM-Q was able to predict the death trends in different phases in the
pandemic with different dynamic ranges. For example, the number of daily
deaths in Los Angeles fluctuated between 20 and 60 per day whereas the counts
in Philadelphia were roughly a third of that.

To validate the performance of the condLSTM-Q, we aggregated the sum of
reported mortality cases of the counties and plotted the national trend with
New York Times’ statistics (Fig 2). At first, the precision of predictive
values was subjected to limited data for model training. The predictions
gradually improved around mid-April, and successfully predicted the descending
trend of COVID-19 mortality from May to July 2020 and the “second wave” of
COVID-19 since late July 2020. To observe the changes in prediction accuracy
from day 1 to day 14 of the prediction (i.e. across the two-week prediction
window), we aligned the day 1, 3, 7, 10, and 14 of each of the prediction
windows to their corresponding dates. The overall prediction accuracy was
improved with more available training data. Before mid-April 2020, distal
prediction (e.g. day 14) was much more varied compared to proximal
predictions; however, after the model became stable with enough training data,
the performance of both distal and proximal prediction converged with a
smaller spread, and the precision of distal prediction was not inferior to
proximal predictions.

Visualization of the nationwide, county-based prediction from August 6 to
August 19, 2020, is shown in Fig 4. This interval approximately covered
periods of the COVID-19 “second wave” with the highest number of daily death
attributed to COVID-19 in the United States, which impacted the Southwest and
Southeast regions the most. With county-level resolution, condLSTM-Q
demonstrated its ability to pick up inhomogeneous hot spots within specific
states, whereas other coarser geographical models could only observe an
average trend of the whole state. In addition, at this resolution, one can
observe state-state interactions of spread in the state border. For instance,
quite a few counties in Arizona and Nevada had trends more similar to their
neighboring California hot spots than to other counties within the states.

###  3.2 Comparison with other models

To evaluate the performance of condLSTM-Q, we compared its predictions with
the IHME model [5]. As condLSTM-Q provides county-wise predictions while the
IHME model generates state-wise predictions, we aggregated predictions of
condLSTM-Q over each county into corresponding states and measured root mean
square error (RMSE) of the state-wise forecasts on two-week intervals. We also
averaged across the 9 quantile predictions from condLSTM-Q to estimate the
mean of the distribution, in order to allow direct comparison against single
value predictions from IHME. Note also that due to the format of the IHME
model available, predictions from the two models are matched on a bimonthly
basis.

As shown in Fig 5, the state-wise predictions of condLSTM-Q were consistently
comparable to the IHME model. The condLSTM-Q initially had higher RMSE in
early May 2020, but showed better performance over most intervals since then
after training data became ample. Despite the fact that condLSTM-Q model was
trained using county-level data with the pinball loss function, it
demonstrated a robust performance on a different geographical scale under a
different metric (i.e. RMSE). We also measured the RMSE by setting zeros as
control(i.e. placing zeros on the full 14-days interval across all 50 states
in the United States). Our results revealed that both IHME and our condLSTM-Q
model presented a similar pattern with the control, indicating that there
could be an uncaptured variance in both models. We also noticed that there was
an abrupt peak of the epidemic curve (approximately 20,000 deaths) on June 24,
2020, when New York Times made bulk adjustments in its data due to changes in
tallying criteria.

Figure 2: Nationwide prediction on the total number of deaths, aggregated from
counties. At each time point, the actual data from that day is shown in blue,
and the predicted values returned from models trained until 1, 3, 7, 10, and
14 days ago are shown in different colors. The model’s overall accuracy is
relatively low in the beginning when there was not enough data for the model
to be trained on, but the prediction started to follow the trend well since
late April Figure 3: Predictions in representative locations Representative
counties including Cook, Los Angeles, New York, Wayne, Philadelphia, Hennepin,
Maricopa, and Montgomery. The condLSTM-Q was able to keep track of the death
trends in different phases in the pandemic with different dynamic ranges. For
instance, Los Angeles had daily deaths fluctuating between 20 - 60 per day
whereas Philadelphia counts were roughly a third of that.

Figure 4: County-wise prediction by condLSTM-Q. County-wise absolute death
counts (Top) and the differences with the ground truth (Bottom), summed from
August 6 to August 19, 2020. This interval was roughly the peak of the “second
wave" of deaths in the United States. Note that the top and the bottom figures
have different color scales.

Dates | 14-day prediction RMSE
---|---
condLSTM-Q |  IHME |  control
2020-05-04 | 54.814 | 38.189 | 67.881
2020-05-19 | 18.396 | 20.812 | 38.046
2020-06-03 | 13.709 | 14.152 | 27.152
2020-06-24 | 74.154 | 74.295 | 77.597
2020-07-04 | 16.291 | 14.790 | 28.945
2020-07-18 | 48.690 | 47.098 | 64.760
2020-08-06 | 32.927 | 33.749 | 55.354
2020-08-21 | 14.276 | 21.912 | 37.628
2020-09-02 | 16.325 | 18.915 | 34.328
2020-09-18 | 13.123 | 14.421 | 28.974

Figure 5: Performance analysis. State-level prediction (Top). State-wise, two-
week prediction RMSE of condLSTM-Q and IHME model. We matched the starting
dates of our predictions to the dates when the IHME model was updated, where
each model would have access to the training data up to the day before the
onset of the two-week prediction. We also tried placing zeros in all the
predictions and calculated RMSE as a control. Both IHME and condLSTM-Q showed
a similar pattern with the control, indicating that there is an uncaptured
variance for both models. County(FIPS)-level performance measured by pinball
loss and RMSE (Bottom). The error for both metrics showed a steady decline
since April with several jumps from June to August, which is from the
abnormality in New York Times’ data due to their bulk adjustments following
the changes in tallying criteria.

###  3.3 Usefulness of conditional layer for categorical data

With a classical LSTM network which does not have a conditional layer for
initializing the hidden states based on categorical data as in condLSTM-Q, the
categorical data should be stacked to the same dimension of the time series
data to be fed into the model. As mentioned in the Introduction, this may lead
to suboptimal performance of the model because it introduces constant values
in time series features which actually is not sequential data. To test this,
we compared the prediction accuracy of two models: the model with the
aforementioned stacking method, which we refer to as the "pseudo-categorical"
LSTM model, and condLSTM-Q. The pseudo-categorical LSTM model had a pinball
loss of $0.115$ during the prediction period of May 22 - June 4, 2020, and
$0.0994$ during June 5 - June 18, 2020. The pinball loss of condLSTM-Q, on the
other hand, was $0.0959$ and $0.0744$, respectively, over the same prediction
periods.

To investigate whether such a decrease in pinball loss was observed in every
county, we looked at the difference in pinball loss between condLSTM-Q and
pseudo-categorical LSTM model in each county. Fig 6 A shows a distribution of
such differences in counties with a high number of deaths ($>50$). The
difference was obtained by subtracting the pinball loss of a pseudo-
categorical model from condLSTM-Q’s loss, so its bias to the left means
condLSTM-Q’s loss is significantly lower than the pseudo-categorical model
($p<0.0005$, one-sided Wilcoxon signed-rank test)

The difference between these two models is well-demonstrated in Fig 6 B and C,
which shows example predictions of May 22 - June 4 and June 5 - June 18, 2020,
in New York - the county with the largest difference in pinball loss between
the models. When the condLSTM-Q was already capturing a down-turned trend, the
pseudo-categorical LSTM was still predicting an upward trend. We also found
that the condLSTM-Q had a much tighter spread of prediction, allowing planning
based on the predictions much more possible.

Figure 6: Effectiveness of conditional architecture. (A) A histogram of
differences in pinball loss between condLSTM-Q and the pseudo-categorical
model, in counties with a high number of total death (>50). The loss from the
pseudo-categorical model was subtracted from condLSTM-Q’s loss, so negative
values mean condLSTM-Q has lower pinball loss. (B, C) Representative case
study of New York, in two different time frames. The condLSTM-Q not only
captures the trend better, but also has a tighter distribution when compared
to the pseudo-categorical model trained on the same training data with
categorical data stacked to match the time series.

###  3.4 Explainability of the model

Figure 7: Importance of each categorical feature. The permutation feature
importance of the 10 most important categorical features in three prediction
windows (May 15-28, July 15-28, and October 8-21, 2020) is shown. The
categorical features were more important during the earlier phase rather than
in the later phase. The list of important features also changed over time: GDP
and heart disease mortality for the earlier phase; population, mortality, and
the number of eligible people for medicare for the later phase. Asterisks
indicate statistical significance based on one-sample sign test to original
validation loss (*p<0.05, **p<0.01) Figure 8: Importance of each time series
feature. The permutation feature importance of the time series features in
three prediction windows (May 15-28, July 15-28, and October 8-21, 2020) is
shown. The number of cases, deaths, and flu seasonality are the three most
important time series features in the model. The importance of seasonality,
however, went down nearly to zero in the later phase. Asterisks indicate
statistical significance from one-sample sign test to original validation loss
(*p<0.05, **p<0.01)

To quantify the relative importance of individual features and see which has
the most information for the prediction, permutation feature importance was
measured for each time series and categorical feature. As suggested in
previous studies [15, 16], we measured the feature importance based on how the
model’s validation loss has increased after permuting each feature in a
validation set. Increased validation loss after shuffling a feature shows that
the model is relying on the feature for the prediction, whereas an unchanged
validation loss means the model is largely ignoring the feature during
prediction.

For categorical data, each feature was shuffled across counties, as it does
not have a temporal component. Time series features, on the other hand, were
shuffled across counties first, and then permuted across time. As each county
can have drastically different values in time series features, shuffling
across counties and time separately ensures more realistic input data for
LSTM, compared to the features shuffled across counties and time altogether.
Such permutation was done 10 times for each feature, and the resulting 10
input data sets were given to the trained model to generate predictions. The
validation loss was calculated from each prediction and compared to the
original validation loss to see how the shuffled features affected the model’s
performance. As the importance of each feature can change over time,
permutation feature importance was measured at three different time frames
(May 15-28, July 15-28, and October 8-21, 2020) with the same procedure.

The result is shown in Fig 7 and Fig 8. For categorical data (Fig 7), its
importance to the prediction was generally higher during the earlier phase,
when there were insufficient time series data for the model to be trained on.
The top three features with the greatest importance during the earlier phase
in May were GDP in 2015, GDP in 2016, and heart disease mortality. GDP was
also one of the most important features in the later phase in July, but it was
not as important as the population estimate (2018) and estimated mortality
(2015-2017). These two features were the most important categorical features
in October as well, followed by the number of eligible people for Medicare
(2018).

In contrast, the ranking of the importance of the time series features (see
Fig 8) was relatively stable compared to categorical features. The number of
deaths and cases were the two most important features in all three time
frames. In the beginning, the number of cases was more important, but was
surpassed by the number of deaths later on. Seasonality also had high
importance, especially during the earlier phase in May and July 2020 when
there were fewer data for the model to be trained on.

##  4 Discussion

###  4.1 Contribution and Novelty

This is the first use of conditional RNN in the prediction of COVID-19 trend
that we know of. The model has the merit of providing fine-scale, county-level
predictions of daily deaths with a two-week window. The performance of the
model was robust and comparable to, if not slightly better than, the well-
known IHME predictions. Given that the application of this model type is very
new in the field, and that we are far from testing all possible architectures
with different input features, one could foresee that the model’s performance
in terms of stability, precision, and length of the forecast window could
still be enhanced with further optimization. We also see that this type of
model could be generalized to trend forecasting of other infectious diseases
in the future and is not limited to COVID-19 itself.

County-level spatial resolution is a feature not available in other famous and
publicly available forecasts such as IHME, DELPHI, and Los Alamos National
Laboratory [17, 18] which provide state-wise predictions. This provides unique
opportunities for investigating trends within the states, and interactions
between counties along state borders. In addition, through the analysis of the
importance of the categorical data which will be discussed below, one could
learn which features are risk factors that affect the death trend and provide
handles for officials to ameliorate the risks.

###  4.2 Effect of categorical features on predictions

Figure 9: Changes in nationwide prediction after modifying categorical
features. For each model making predictions over 2 weeks from May 15th (top),
July 15th (middle), and October 8th, 2020 (bottom), new predictions after
increasing and decreasing a categorical feature are shown in blue and orange,
along with the original prediction in green. Each altered feature has been
selected based on permutation feature importance from each of the models, and
features were either increased or decreased by 3 standard deviations. The
values shown in this figure are the predicted numbers of deaths in the United
States.

One advantage of using conditional LSTM is that it provides insight into the
effect of each categorical and time series feature on the predictions. Fig 9
shows the total number of deaths predicted by each model after altering
categorical features with the highest permutation feature importance. As
mentioned in the Result, GDP in 2015 and 2016 and heart disease mortality were
the most important features in the beginning, but later on, the dominant
factors shifted to the population estimate, estimated mortality, and the
number of eligible people for Medicare; an increase in any of these features
led to an increase in predicted counts and vice versa.

Such a pattern in the earlier phase shows the possible vulnerability of people
with heart disease to COVID-19, which is consistent with previous studies
reporting elevated risk of poorer outcomes with COVID-19 in people with
congenital heart disease [19] and coronary heart disease [20, 21, 22]. The
positive correlation between GDP and the predicted number of deaths in the
earlier phase is also consistent with a study claiming a positive correlation
between GDP and the number of confirmed cases of COVID-19 in China [23] \-
nevertheless, the GDP estimates the value of goods and services produced by
each county [24], which inevitably correlates with its population and other
variables. Therefore, it would require further investigation to remove the
dependencies of each factor.

In the later phase, on the other hand, the important features have shifted to
population estimates, mortality rate, and the number of eligible people for
Medicare. This may indicate that the number of deaths due to COVID-19 in the
later phase is more related to the overall healthcare resources and general
health condition of the population in each county.

###  4.3 Additional Benefits

The LSTM-based structure is scalable and flexible. We were limited to a few
available data sources during the initiative and have not expanded actively to
other data sources that possibly provide additional information for
predictions, as testing the robustness of this architecture was the original
goal of this pilot study. Once there are additional data sources, the
condLSTM-Q could adapt and be retrained to extract information from the new
data source without adjusting the architecture. For instance, if one
hypothesizes that the temperature affects the death trend, the model can be
easily retrained and tested with an additional time series feature of
temperature added to the input data stream. One could also foresee that, when
the vaccine is eventually distributed, the immunization rate in regions could
be one additional important feature.

##  5 Material Availability

Code and Model for our study are available in:
https://github.com/cjackal/COVID-SKTW

Demographics and Health Resource dataset from the Yu Group at UC Berkeley:
https://raw.githubusercontent.com/Yu-Group/covid19-severity-
prediction/data/county_data_abridged.csv

COVID-19 mortality dataset from New York Times:
https://raw.githubusercontent.com/nytimes/covid-19-data/master/us-counties.csv

Local area GDP dataset from U.S. Bureau of Economic Analysis (BEA):
https://www.bea.gov/news/2019/local-area-gross-domestic-product-2018

Land area and Population density dataset from 2010 Census:
https://factfinder.census.gov/faces/tableservices/jsf/pages/productview.xhtml?src=bkmk

Mobility dataset from DescartesLabs:
https://raw.githubusercontent.com/descarteslabs/DL-COVID-19/master/DL-us-
mobility-daterow.csv

COVID-19 related policies dataset from HHS:
https://healthdata.gov/dataset/covid-19-state-and-county-policy-orders

##  6 Acknowledgement

The authors thank Prof Yaser Abu-Mostafa, and the Teaching Assistants of CS156
in Caltech for organizing the COVID19 prediction initiative and for providing
the data pipeline for parsing data sources. We thank Isaac Yen-Hao Chu M.D.
for reading the manuscript. Y.L.Ni was supported by Taipei Veterans General
Hospital-National Yang-Ming University Excellent Physician Scientists
Cultivation Program, No.103-Y-A-003.

## References

  * [1] IHME COVID-19 health service utilization forecasting team and Christopher JL Murray.  Forecasting covid-19 impact on hospital bed-days, icu-days, ventilator-days and deaths by us state in the next 4 months.  medRxiv, 2020.
  * [2] Yaser S. Abu-Mostafa.  Caltech cs156 covid-19 model.  https://cs156.caltech.edu/, September 2020.
  * [3] Najeebullah Khan, Shamsuddin Shahid, Liew Juneng, Kamal Ahmed, Tarmizi Ismail, and Nadeem Nawaz.  Prediction of heat waves in pakistan using quantile regression forests.  Atmospheric Research, 221:1 – 11, 2019.
  * [4] Angela Brennan, Paul C. Cross, and Scott Creel.  Managing more than the mean: using quantile regression to identify factors related to large elk groups.  Journal of Applied Ecology, 52(6):1656–1664, 2015.
  * [5] IHME.  Institute for health metrics and evaluation covid-19 estimate.  http://www.healthdata.org/covid/data-downloads, October 2020.
  * [6] Hochreiter and Schmidhuber.  Long Short-Term Memory.  Neural Computation, 9(8):1735–1780, November 1997.  Publisher: MIT Press.
  * [7] Ratnabali Pal, Arif Ahmed Sekh, Samarjit Kar, and Dilip K. Prasad.  Neural network based country wise risk prediction of COVID-19.  Applied Sciences, 10(18):6448, September 2020.  arXiv: 2004.00959.
  * [8] Patricia Melin, Julio Cesar Monica, Daniela Sanchez, and Oscar Castillo.  Multiple Ensemble Neural Network Models with Fuzzy Response Aggregation for Predicting COVID-19 Time Series: The Case of Mexico.  Healthcare, 8(2):181, June 2020.  Number: 2 Publisher: Multidisciplinary Digital Publishing Institute.
  * [9] Philippe Rémy.  Conditional rnns made easy with tensorflow and keras.  https://github.com/philipperemy/cond_rnn, March 2020.
  * [10] Andrej Karpathy and Fei-Fei Li.  Deep visual-semantic alignments for generating image descriptions.  CoRR, abs/1412.2306, 2014.
  * [11] The Covidvis team.  Covidvis.  https://covidvis.berkeley.edu/, May 2020.
  * [12] CDC.  Pneumonia and influenza mortality surveillance from the national center for health statistics mortality surveillance system.  https://gis.cdc.gov/grasp/fluview/mortality.html, October 2020.
  * [13] Hans-Georg Zimmermann, Christoph Tietz, and Ralph Grothmann.  Forecasting with recurrent neural networks: 12 tricks.  In Neural Networks: Tricks of the Trade, pages 687–707. Springer, 2012.
  * [14] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * [15] Leo Breiman.  Random forests.  Machine learning, 45(1):5–32, 2001.
  * [16] Aaron Fisher, Cynthia Rudin, and Francesca Dominici.  All models are wrong, but many are useful: Learning a variable’s importance by studying an entire class of prediction models simultaneously, 2019\.
  * [17] COVID Analytics.  Delphi model.  https://www.covidanalytics.io/projections/, October 2020.
  * [18] Los Alamos National Laboratory.  Lanl model.  https://covid-19.bsvgateway.org/, October 2020.
  * [19] Weiyi Tan and Jamil Aboulhosn.  The cardiovascular burden of coronavirus disease 2019 (covid-19) with a focus on congenital heart disease.  International Journal of Cardiology, 309:70 – 77, 2020.
  * [20] Ruchong Chen, Wenhua Liang, Mei Jiang, Weijie Guan, Chen Zhan, Tao Wang, Chunli Tang, Ling Sang, Jiaxing Liu, Zhengyi Ni, Yu Hu, Lei Liu, Hong Shan, Chunliang Lei, Yixiang Peng, Li Wei, Yong Liu, Yahua Hu, Peng Peng, Jianming Wang, Jiyang Liu, Zhong Chen, Gang Li, Zhijian Zheng, Shaoqin Qiu, Jie Luo, Changjiang Ye, Shaoyong Zhu, Xiaoqing Liu, Linling Cheng, Feng Ye, Jinping Zheng, Nuofu Zhang, Yimin Li, Jianxing He, Shiyue Li, and Nanshan Zhong.  Risk factors of fatal outcome in hospitalized subjects with coronavirus disease 2019 from a nationwide analysis in china.  Chest, 158(1):97 – 105, 2020.
  * [21] Yingyu Chen, Xiao Gong, Lexun Wang, and Jiao Guo.  Effects of hypertension, diabetes and coronary heart disease on covid-19 diseases severity: a systematic review and meta-analysis.  medRxiv, 2020.
  * [22] C Chen, C Chen, JT Yan, N Zhou, JP Zhao, and DW Wang.  Analysis of myocardial injury in patients with covid-19 and association between concomitant cardiovascular diseases and severity of covid-19.  Zhonghua xin xue guan bing za zhi, 48(7):567—571, July 2020.
  * [23] Yi Zhang, Hanwen Tian, Yinglong Zhang, and Yiping Chen.  Is the epidemic spread related to gdp? visualizing the distribution of covid-19 in chinese mainland, 2020.
  * [24] U.S. Bureau of Economic Analysis (BEA).  Local area gross domestic product, 2018.  https://www.bea.gov/system/files/2019-12/lagdp1219.pdf, October 2020\.


$\displaystyle\hskip
28.45274pt+\big{(}z^{2}-1+y^{2}\big{)}\Big{(}\frac{1}{2}-iM\Big{)}F\Big{(}-\frac{1}{2}+iM,\frac{1}{2}+iM;1;\frac{(z-1)^{2}-y^{2}}{(z+1)^{2}-y^{2}}\Big{)}\Bigg{|}\,dy$
$\displaystyle\leq$ $\displaystyle
Cz^{-\frac{1}{2}}(z-1)^{1+2s-n(\frac{1}{p}-\frac{1}{q})}\quad\mbox{\rm for
all}\quad z\in[1,\infty)\,.$

In the second zone we use (7.10), (7.11), and (7.14). Thus, we have to
estimate the next two integrals:

$\displaystyle A_{11}$ $\displaystyle:=$ $\displaystyle z^{2}\int_{(z,r)\in
Z_{2}(\varepsilon,z)}r^{2s-n(\frac{1}{p}-\frac{1}{q})}\frac{1}{((z-1)^{2}-r^{2})\sqrt{(z+1)^{2}-r^{2}}}\,dr\,,$
$\displaystyle A_{12}$ $\displaystyle:=$ $\displaystyle z^{2}\left(1+\ln
z\right)\int_{(z,r)\in
Z_{2}(\varepsilon,z)}r^{2s-n(\frac{1}{p}-\frac{1}{q})}\frac{1}{((z-1)^{2}-r^{2})\sqrt{(z+1)^{2}-r^{2}}}\,dr\,.$

We apply (7.10) to $A_{11}$ and obtain

$\displaystyle A_{11}\leq C_{\varepsilon}z^{2}\int_{(z,r)\in
Z_{2}(\varepsilon,z)}r^{2s-n(\frac{1}{p}-\frac{1}{q})}\frac{1}{[(z+1)^{2}-r^{2}]}\frac{1}{\sqrt{(z+1)^{2}-r^{2}}}\,dr\leq
C_{\varepsilon}z^{-\frac{1}{2}}(z-1)^{1+2s-n(\frac{1}{p}-\frac{1}{q})}$

for all $z\in[1,\infty)$, while

$\displaystyle A_{12}$ $\displaystyle\leq$ $\displaystyle
C_{\varepsilon}z^{-\frac{1}{2}}(z-1)^{1+2s-n(\frac{1}{p}-\frac{1}{q})}\left(1+\ln
z\right)\quad\mbox{\rm for all}\quad z\in[1,\infty)\,.$

Proposition is proven. $\square$

To complete the proof of the theorem we write

$\displaystyle\int_{0}^{1-e^{-t}}r^{2s-n(\frac{1}{p}-\frac{1}{q})}|K_{0}(r,t)|\,dr$
$\displaystyle\leq$ $\displaystyle
e^{-t[2s-n(\frac{1}{p}-\frac{1}{q})]}\int_{0}^{e^{t}-1}y^{2s-n(\frac{1}{p}-\frac{1}{q})}\frac{1}{[(e^{t}-1)^{2}-y^{2}]\sqrt{(e^{t}+1)^{2}-y^{2}}}$
$\displaystyle\times\Bigg{|}\big{(}e^{t}-e^{2t}-iM(1-e^{2t}-y^{2})\big{)}F\Big{(}\frac{1}{2}+iM,\frac{1}{2}+iM;1;\frac{(e^{t}-1)^{2}-y^{2}}{(e^{t}+1)^{2}-y^{2}}\Big{)}$
$\displaystyle\hskip
28.45274pt+\big{(}e^{2t}-1+y^{2}\big{)}\Big{(}\frac{1}{2}-iM\Big{)}F\Big{(}-\frac{1}{2}+iM,\frac{1}{2}+iM;1;\frac{(e^{t}-1)^{2}-y^{2}}{(e^{t}+1)^{2}-y^{2}}\Big{)}\Bigg{|}\,dy$
$\displaystyle\leq$ $\displaystyle
Ce^{-t[\frac{1}{2}+2s-n(\frac{1}{p}-\frac{1}{q})]}(e^{t}-1)^{1+2s-n(\frac{1}{p}-\frac{1}{q})}(1+t)\,.$

Thus,

$\displaystyle\|(-\bigtriangleup)^{-s}u(x,t)\|_{{L}^{q}({\mathbb{R}}^{n})}$
$\displaystyle\leq$ $\displaystyle
C\Big{(}e^{\frac{t}{2}}(1-e^{-t})^{2s-n(\frac{1}{p}-\frac{1}{q})}+e^{-t[\frac{1}{2}+2s-n(\frac{1}{p}-\frac{1}{q})]}(e^{t}-1)^{1+2s-n(\frac{1}{p}-\frac{1}{q})}(1+t)\Big{)}\|\varphi_{0}(x)\|_{{L}^{p}({\mathbb{R}}^{n})}$
$\displaystyle\leq$ $\displaystyle
C(1+t)e^{\frac{t}{2}}(1-e^{-t})^{2s-n(\frac{1}{p}-\frac{1}{q})}\|\varphi_{0}(x)\|_{{L}^{p}({\mathbb{R}}^{n})}.$

Theorem is proven. $\square$

## References

  * [1] Abramowitz, M., Stegun, I.A.: Handbook of mathematical functions with formulas, graphs, and mathematical tables. National Bureau of Standards Applied Mathematics Series, 55, Washington, DC, 1964
  * [2] Bateman, H., Erdelyi, A.: Higher Transcendental Functions. v.1,2, McGraw-Hill, New York, 1953.
  * [3] Birrell, N. D., Davies, P.C.W. : Quantum fields in curved space, Cambridge ; New York : Cambridge University Press, 1984.
  * [4] Bony, J.-F., Hafner, D.: Decay and non-decay of the local energy for the wave equation in the De Sitter - Schwarzschild metric, arXiv:0706.0350v1
  * [5] Brenner, P.: On $\,L^{p}-L^{q}\,$ estimates for the wave-equation. Math. Zeitschrift 145 (1975) 251-254.
  * [6] Brevik, I., Simonsen, B.: The scalar field equation in Schwarzschild-de Sitter space. Gen. Relativity Gravitation 33 (2001), no. 10, 1839–1861.
  * [7] Brozos-Vázquez, M., García-Río, E., Vázquez-Lorenzo, R.: Locally conformally flat multidimensional cosmological models and generalized Friedmann-Robertson-Walker spacetimes, J. Cosmol. Astropart. Phys. JCAP12(2004)008 doi:10.1088/1475-7516/2004/12/008.
  * [8] Dafermos, M., Rodnianski, I.: The wave equation on Schwarzschild-de Sitter spacetimes,arXiv:0709.2766.
  * [9] Chandrasekhar, S.: The Mathematical Theory of Black Holes, Oxford,Clarendon Press ; New York : Oxford University Press, 1998.
  * [10] Chrusciel, P.T., Pollack,D.: Singular Yamabe metrics and initial data with exactly Kottler-Schwarzschild-de Sitter ends, arXiv:0710.3365v1
  * [11] De Sitter, W.: On Einstein’s Theory of Gravitation, and its astronomical consequences.II,III. Royal Astronimcal Society. 77 (1917) 155-184; 78 (1917) 3-28.
  * [12] Einstein, A.: Kosmologische Betrachtungen zur allgemeinen Relativitätstheorie, Sitzungsber Preuss. Akad. Wiss. Berlin (1917) 142-152.
  * [13] Finster, F., Kamran, N., Smoller, J., Yau, S.-T.: Decay of solutions of the wave equation in the Kerr geometry. Comm. Math. Phys. 264 (2006), no. 2, 465–503.
  * [14] Friedrich, H., Rendall, A.: The Cauchy problem for the Einstein equations. Einstein’s field equations and their physical implications. Lecture Notes in Phys., 540, Springer, Berlin (2000) 127–223.
  * [15] Galstian, A.: $L_{p}$-$L_{q}$ decay estimates for the wave equations with exponentially growing speed of propagation. Appl. Anal. 82 (3) (2003) 197–214.
  * [16] Heinzle, J. M., Rendall, A., Power-law inflation in spacetimes without symmetry. Commun. Math. Phys. 269 (2007) 1-15.
  * [17] Hörmander, L.: The analysis of linear partial differential operators. IV. Fourier integral operators.  Grundlehren der Mathematischen Wissenschaften, 275. Springer-Verlag, Berlin, 1994.
  * [18] Kronthaler, J.: The Cauchy problem for the wave equation in the Schwarzschild geometry. J. Math. Phys. 47 (2006), no. 4, 042501, 29 pp.
  * [19] Littman, W.: The wave operator and $L_{p}$ norms. J. Math. Mech. 12 (1963) 55–68.
  * [20] Littman, W., McCarthy, C., Rivière, N.: The non-existence of $L^{p}$ estimates for certain translation-invariant operators. Studia Math. 30 (1968) 219–229.
  * [21] López-Ortega, A.: Quasinormal modes of $D$-dimensional de Sitter spacetime. Gen. Relativity Gravitation 38 (2006), no. 11, 1565–1591.
  * [22] M$\o$ller, C.: The theory of relativity. Oxford, Clarendon Press, 1952.
  * [23] Pecher, H.: $L^{p}$-Abschätzungen und klassische Lösungen für nichtlineare Wellengleichungen.I, Math. Zeitschrift. 150 (1976) 159-183.
  * [24] Peral, J.C.: $L^{p}$ estimates for the wave equation. J. Funct. Anal. 36 (1) (1980) 114-145.
  * [25] Racke, R.: Lectures on Nonlinear Evolution Equations. Aspects of Mathematics, Vieweg, Braunschweig/Wiesbaden, 1992.
  * [26] Rendall, A.: Asymptotics of solutions of the Einstein equations with positive cosmological constant. Ann. Henri Poincaré 5 (6) (2004) 1041-1064.
  * [27] Shatah, J., Struwe, M.: Geometric wave equations. Courant Lecture Notes in Mathematics, 2. New York University, Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 1998.
  * [28] Sonego, S., Faraoni, V.: Huygens’ principle and characteristic propagation property for waves in curved space-times, J. Math. Phys. 33(2) (1992) 625-632.
  * [29] Yagdjian, K.: A note on the fundamental solution for the Tricomi-type equation in the hyperbolic domain, J. Differential Equations 206 (2004) 227-252.
  * [30] Yagdjian, K.: Global existence in the Cauchy problem for nonlinear wave equations with variable speed of propagation, New trends in the theory of hyperbolic equations, 301-385, Oper. Theory Adv. Appl., 159, Birkh$\ddot{\rm a}$user, Basel, 2005.
  * [31] Yagdjian, K.: Global existence for the $n$-dimensional semilinear Tricomi-type equations, Comm. Partial Diff. Equations 31 (2006) 907-944.
  * [32] Yagdjian, K.: Self-similar solutions of semilinear wave equation with variable speed of propagation. J. Math. Anal. Appl. 336 (2007) 1259-1286..
  * [33] Yagdjian, K. Galstian, A.: Fundamental Solutions for Wave Equation in de Sitter Model of Universe, ISSN 1437-739X, University of Potsdam, August, Preprint 2007/06.


# Coupling trapped ions to a nanomechanical oscillator

Moritz Weegen  Department of Chemistry, University of Basel,
Klingelbergstrasse 80, 4056 Basel, Switzerland  Swiss Nanoscience Institute,
University of Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland  Martino
Poggio  Department of Physics, University of Basel, Klingelbergstrasse 82,
4056 Basel  Swiss Nanoscience Institute, University of Basel,
Klingelbergstrasse 82, 4056 Basel, Switzerland  Stefan Willitsch
<EMAIL_ADDRESS>Department of Chemistry, University of Basel,
Klingelbergstrasse 80, 4056 Basel, Switzerland  Swiss Nanoscience Institute,
University of Basel, Klingelbergstrasse 82, 4056 Basel, Switzerland

###### Abstract

Cold ions in traps are well-established, highly controllable quantum systems
with a wide variety of applications in quantum information, precision
spectroscopy, clocks and chemistry. Nanomechanical oscillators are used in
advanced sensing applications and for exploring the border between classical
and quantum physics. Here, we report on the implementation of a hybrid system
combining a metallic nanowire with laser-cooled ions in a miniaturised ion
trap. We demonstrate resonant and off-resonant coupling of the two systems and
the coherent motional excitation of the ion by the mechanical drive of the
nanowire. The present results open up avenues for mechanically manipulating
the quantum motion of trapped ions, for the development of ion-mechanical
hybrid quantum systems and for the sympathetic cooling of mechanical systems
by trapped ions and vice versa.

ion traps, nanowires, nanomechanical oscillators, hybrid quantum systems,
quantum technologies

††preprint: APS/123-QED

Laser-cooled ions stored in harmonic traps [1] count among the best controlled
quantum systems [2]. Single ions or strings of ions can be cooled to the
quantum regime of motion in the trap, thus realising quantum harmonic
oscillators which find a variety of applications across different branches of
sciences and technology. The ability to prepare tailored quantum states of
motion and entangle them with the ions’ internal degrees of freedom, as well
as their excellent coherence properties [3], render them a leading platform
for quantum computing [4, 5, 6] and quantum simulation [7, 8, 9]. Similarly,
trapped ions form the basis of some of today’s most accurate clocks [10] and
have enabled new types of ion-neutral collision and chemical experiments at
very low energies [11, 12, 13]. Their unique properties have also been
harnessed in hybrid experiments with ultracold atoms [14, 15].

Figure 1: Ion-nanowire hybrid system. a, Schematic of the experimental setup.
An Ag2Ga nanowire was positioned close to trapped 40Ca+ ions in a miniaturised
RF trap. The ions were cooled with laser beams of wavelengths 397 nm and 866
nm. Oscillations of the nanowire were imaged by projecting its shadow onto a
dual-element photodiode with a laser beam at 866 nm. The resonance
fluorescence of the ions generated during laser cooling was imaged by an EMCCD
camera coupled to a microscope above the experimental chamber. A beam splitter
directed half of the collected ion fluorescence to a photomultiplier tube
(PMT) for the detection of driven ion motion using photon correlation. See
text for further details. b, Fluorescence image of a single trapped ion at
rest (top) and driven by the nanowire (bottom). c, Optical microscope image of
the nanowire and its support at $20\times$ magnification.

Traditionally, trapped ions have been manipulated using laser [2], electric
[16, 17], magnetic and microwave [18, 19] fields. An intriguing prospect is
their manipulation by other types of oscillators, e.g., other trapped ions
[20, 21, 22] or nanomechanical systems. As explored theoretically in Refs.
[23, 24, 25, 26, 27, 28], ion-mechanical hybrid systems offer possibilities
for extending the available techniques for the cooling, manipulation, control
and readout of both constituents up to realising novel hybrid quantum systems
[29, 30] in which entanglement between ions and mechanical objects can be
realised.

In this context, the research on mechanical oscillators on the nanometer scale
has progressed rapidly in recent years. Their properties as objects on the
border between classical and quantum physics make them excellent candidates
for the realisation of ion-mechanical hybrid systems [31]. Cryogenic cooling
of nanomechanical oscillators in combination with other cooling techniques has
enabled the leap from the classical to the quantum regime by preparing them
close to their motional ground state [32, 33, 34, 35, 36]. They have also been
successfully coupled to ultracold atoms [37, 38, 39, 40, 41].

Here, we present a hybrid system consisting of trapped, laser-cooled 40Ca+
ions in a miniaturised linear radiofrequency (RF) ion trap coupled to a
conductive nanowire [42] as previously explored theoretically in Ref. [28].
The present hybrid system was realised with single trapped ions, strings of
several ions or larger ion ensembles forming three-dimensional Coulomb
crystals [11]. A bias voltage applied to the nanowire coupled the two systems
by their mutual electrostatic interaction. Mechanically driving the
oscillation of the nanowire modulated the electric field experienced by the
ions leading to an effective transfer of energy when the drive was resonant
with the motion of the ions in the trap.

Fig. 1 shows a schematic of the experimental setup with the most important
components (see Methods for details). Ions were confined and laser-cooled in a
microstructured linear-quadrupole radiofrequency trap [11] consisting of
electrodes in a stacked-wafer configuration. The ions were observed by imaging
their resonance fluorescence generated during laser cooling onto an EMCCD
camera coupled to a microscope (Fig. 1b). Motional excitation of the ions was
quantified using the photon correlation method [43, 44, 45].

The nanomechanical oscillator used was a Ag2Ga cantilever positioned at the
tip of a tungsten support with $250$ µm base diameter and $\approx 4.8$ mm
length. The length of the nanowire was measured with a microscope to be
$\approx 16$ µm with a diameter of $\approx 150$ nm as specified by the
manufacturer. Fig. 1c shows a microscope image of the nanowire on the tungsten
tip. The oscillator was mounted on three nanopositioners allowing its free
deployment in all spatial directions inside the trap. The wire was
mechanically driven in the direction of the longitudinal trap axis ($z$) by a
piezoelectric actuator attached to the assembly.

Figure 2: Oscillations of the nanowire and the ions. a, Mechanical excitation
spectrum of the nanowire obtained by imaging its vibration with a dual-element
photodiode. The features represent excitations of the eigenmodes $n$ of the
nanowire on its support. The strong $n=4$ mode at 422 kHz was used for the
excitation of the ion. Black dashed lines correspond to eigenfrequencies
obtained from theoretical modelling. b, Histograms of arrival times of
resonance-fluorescence photons of the ions correlated with the periodic drive
signal of the nanowire under near-resonant (top) and far off-resonant (bottom)
motional excitation by the mechanical oscillator. The black dashed lines show
fits of the data from which the velocity amplitude $v_{\text{max}}$ of the
ions during their oscillation was obtained (Methods).

The eigenmodes and, therefore, also the vibration frequencies of mechanical
oscillators are determined by their mass, material and geometry. To
characterise the vibrational modes of the present nanowire, the drive
frequency $f_{\text{drive}}$ applied to the piezo was scanned. Its
oscillations were measured by imaging the tip of the oscillator placed in the
waist of a focused laser beam onto a dual-element Si PIN photodiode [46, 47]
(Methods). Fig. 2a shows a mechanical excitation spectrum thus obtained in the
frequency interval from 100 to 500 kHz.

In this range, the strongest mechanical response was found around
$f_{\text{drive}}=422$ kHz. Structural mechanics simulations performed with
the COMSOL Multiphysics program [48] indicate that this resonance corresponds
to the $n=4$ eigenmode [49] which represents a combined oscillation of the
nanowire and its support (Methods).

To demonstrate the coupling of the ion-nanowire system, the motional
excitation of the ions by a resonant drive of the mechanical oscillator was
explored. For this purpose, the piezo frequency was set to the mechanical
resonance of the nanowire at 422 kHz. The oscillation frequency of the ion
along the longitudinal trap axis was matched by suitably adjusting the static
potentials applied to the electrodes (Methods). Close to resonance between the
nanowire and ion oscillations, both systems coupled efficiently resulting in a
transfer of energy from the mechanical oscillator to the ion (bottom image in
Fig. 1b) [28].

During driven oscillations in the trap, the ions experienced periodically
varying Doppler shifts with respect to the cooling lasers modulating their
resonance fluorescence. The arrival times of fluorescence photons on a
photomultiplier tube (PMT) were correlated with the periodic piezo drive at
$f_{\text{drive}}$ using lock-in detection. Ion velocities were then extracted
from the time-resolved fluorescence profiles (Methods). Fig. 2b shows typical
time-resolved fluorescence traces for a near-resonant
($f_{\text{drive}}=421.0$ kHz, top panel) and and far off-resonant
($f_{\text{drive}}=366.2$ kHz, bottom panel) drive of the ion motion. In these
measurements, a voltage $V_{\text{nw}}=1.2$ V was applied to the nanowire
which was sinusoidally driven with a voltage amplitude $V_{\text{piezo}}=5$ V
applied to the piezo. The equilibrium distance between ion and nanowire was
determined from a fit to experimental data (see below) as $d\approx 300$ µm. A
clear modulation of the fluorescence yield can be seen at strong excitation
while the signal appears largely unmodulated with an off-resonant drive. The
maximum ion velocities $v_{\text{max}}$ during oscillation were determined to
be $v_{\text{max}}~{}=~{}55.3\pm 2.0$ m/s and $1.1\pm 2.2$ m/s for the
examples shown at the top and bottom in Fig. 2b, respectively.

To confirm that the motion of the ion was indeed a consequence of the
mechanical action of the nanowire rather than a spurious excitation from
oscillating stray electric fields in the experiment, the frequency of the
voltage applied to the piezo to drive the nanowire was scanned across the
strong mechanical resonance at 422 kHz observed in Fig. 2a. In each scan step,
the oscillation frequency of the ion along the $z$ axis was matched to the
piezo drive frequency with an accuracy of $\pm 0.2$ kHz by adjusting the
static trapping potentials. The motional excitation of the trapped ions was
quantified as above by measuring their maximum velocities $v_{\text{max}}$
using the photon correlation method. Fig. 3 shows the squared maximum
velocities of the ions $v_{\text{max}}^{2}$ (red crosses), which are
proportional to their maximum kinetic energy in the trap, in relation to the
piezo drive frequency and superimposed onto the mechanical drive spectrum of
the nanowire (orange trace). The measurement was performed with a string of
two ions.

Figure 3: Mechanical excitation of the ion by the nanowire. Squared maximum
velocities $v^{2}_{\text{max}}$ (red crosses) of two trapped ions over a range
of drive frequencies $f_{\text{drive}}$ applied to the nanowire. The orange
trace corresponds to the mechanical excitation spectrum of the nanowire shown
in Fig. 2a. Motional excitation of the ion could only be observed around the
mechanical resonance of the nanowire demonstrating the resonant coupling of
the two systems. Classical dynamics simulations (purple dashed line) reproduce
the frequency response of the ion excitation at the specific experimental
parameters. All traces were normalised to the experimental data for
comparison. Error bars represent a combination of the errors of the fits for
determining $v_{\text{max}}$ from time-resolved ion-fluorescence profiles and
the statistical standard errors of 3 measurements.

Fig. 3 shows that the motion of the ions is only excited when the piezo drives
the motion of the mechanical oscillator. This direct correlation of the ion
motion with the vibration of the nanowire proves the mechanical nature of the
ion excitation and thus the successful coupling between the two systems. The
dashed purple curve in Fig. 3 represents the results of a classical simulation
of the motional excitation of a single ion by the action of the nanowire. The
simulation treats the excitation of the ion by the nanooscillator by
approximating the system as two oscillating point charges corresponding to the
leading term in a multipole-expansion of their mutual electrostatic
interactions (Ref. [28] and Methods). While this simplified treatment cannot
capture the asymmetric lineshape of the nanowire’s mechanical resonance (which
we attribute to overlapping mechanical modes of the entire assembly), it
reproduces the general features of the frequency response of the excitation of
the ion by the nanooscillator. The modelling also shows that the degree of
motional excitation of the ion results from the balance of its continuous
cooling by the lasers with the transfer of energy from the nanowire.

Figure 4: Variation of coupling parameters. Maximum ion velocities
$v_{\text{max}}$ as a function of a. the voltage $V_{\text{piezo}}$ of the
piezo actuator driving the nanowire, b. the ion-nanowire distance $d$, c. the
voltage $V_{\text{nw}}$ applied to the nanowire, and d. the number of ions in
the trap. Dashed lines show the results of a classical simulation of the
system, see text for details. Error bars represent a combination of the errors
of the fits for determining $v_{\text{max}}$ from time-resolved ion-
fluorescence profiles and the statistical standard errors of 3 measurements.

To further explore the properties of the coupled system, its response to
variations of salient experimental parameters was studied. Fig. 4a and b show
the maximum velocities of the ions excited by the nanooscillator on the 422
kHz resonance as a function of the piezo drive voltage and the ion-nanowire
distance, respectively. As expected, the excitation of the ions is increased
with increasing the oscillation amplitude of the nanowire, i.e., with
increasing piezo voltage, and decreased with increasing distance between the
two subsystems. The dashed blue lines show least-squares fits of an
approximated classical simulation model to the combined data. In the fit, the
amplitude $A$ of the nanowire oscillation and the ion-nanowire distance $d$
were adjusted as free parameters (Methods). The good agreement illustrates
that the dynamics of the system and its dependence on these specific
parameters is well captured by the theory.

Fig. 4c shows the variation of the ion’s maximum velocity subject to the
resonant drive at varying bias voltages applied to the nanowire. The coupling,
and hence the motional excitation of the ion, seem to increase linearly with
the charge on the nanowire in the voltage range studied. The dashed blue line
shows the results of simulations assuming the same parameters as used in Figs.
4a,b. While the simulation captures the trend correctly, it predicts a weaker
dependence of the ion velocity on the nanowire voltage than observed
experimentally. We attribute this discrepancy to the simplifications of the
model, in particular to the neglect of asymmetries in the charge distribution
of the nanowire on its holder, i.e., to higher-order terms in the multipole
expansion of the electrostatic interaction between the two systems, and to the
possible presence of stray electric fields in the setup.

Finally, Fig. 4d shows the dependence of the maximum velocity on the number of
trapped ions subjected to the drive by the nanooscillator. It can be observed
that the kinetic energy imparted to each ion is independent of the other ions
in the trap, consistent with the notion that the center-of-mass-motion of the
entire ion crystal is excited by the nanowire. This is also supported by the
ion images which show a uniform blurring of the fluorescence and hence uniform
excitation of multi-ion crystals in this scenario.

The present study demonstrated the coupling of trapped ions to a
nanomechanical oscillator in a novel hybrid system consisting of a
miniaturised linear radiofrequency ion trap and a nanowire assembly. In the
present experiments, the nanowire effectively acted as a miniaturised,
movable, oscillating trap electrode which resonantly couples to the ions in a
highly controllable fashion. Resonant excitation of the motion of trapped
ions, ranging from single ions to small Coulomb crystals, was achieved in the
classical regime and control over the ion motion was demonstrated by the
variation of different coupling parameters.

The present results lay the basis for further experiments. As explored
theoretically in Ref. [28], cooling the ion into the quantum regime of motion
in the trap by, e.g., resolved-sideband methods [2] paves the way for
generating different quantum states of ion motion through the action of the
mechanical oscillator. These include, among others, large coherent states and
multi-component Schrödinger-cat states which are challenging to engineer by
conventional optical means [50]. Additionally, using nanooscillators in the
quantum regime, e.g., low-mass oscillators like carbon nanotubes in a
dilution-refrigerator environment, offers perspectives for realising a
completely quantum hybrid system in which the exchange of discrete phonons,
quantum entanglement and mutual sympathetic cooling can be realised [28].

In this context, the key advantages of using nanomechanical oscillators are
the capability for laserless state manipulation of the ions, the easy control
over the mechanical drive strength as well as the possibility to quickly
adjust the coupling strength on-the-fly by the change of parameters such as
the nanowire’s position and applied voltages. Moreover, their small size
renders them attractive for integration in scalable miniaturised devices as
required for large ion quantum networks [6]. In addition to these prospects,
the present experiment can be further developed in a variety of directions.
Driving the nanowire optically may offer more localised control over its
oscillation without driving other parts of the assembly such as its holder
[46]. Integrated optics could improve the optical in-situ readout of the
nanowire resonances. Lastly, the use of other types of nanomechanical
oscillators such as membranes as candidates for the mechanical quantum-state
manipulation of ions could be investigated.

###### Acknowledgements.

We thank Dr. Anatoly Johnson, Georg Holderied, Philipp Knöpfel and Grischa
Martin for technical assistance. We thank Dr. Panagiotis Fountas for his
contributions to the development of the experimental setup and Dr. Adrien
Poindron for useful discussions. This work is supported by the Swiss
Nanoscience Institute (grant nr. P1808), the Swiss National Science Foundation
(grant nr. 200021_204123), and the University of Basel.

## Methods

### Ion trapping and cooling

The linear-quadrupole RF ion trap consisted of four wafers at a separation of
400 µm fabricated from laser-cut aluminium oxide sputtered with gold following
the layout proposed in Ref. [28]. Two diagonally opposed wafers formed the RF
electrodes. The remaining two wafers were sectioned into seven individually
addressable segments of 400 µm length to which static voltages were applied
for confining the ions along the longitudinal axis of the trap and
manipulating the position of the ions. Typical static voltages applied were on
the order of 0.5 V to 1 V for stable ion confinement along the longitudinal
trap axis at frequencies $f_{z}\approx 400$ kHz. The secular oscillation
frequencies of the ions were measured using resonant motional excitation by
additional RF fields applied to one of the trap electrodes [51]. The
radiofrequency used was $f_{\text{RF}}=21.629$ MHz at an amplitude
$V_{\text{RF}}\approx 50$ V.
40Ca+ ions were loaded by photoionisation of Ca atoms emanating from an oven
source using two diode laser beams at 423 nm and 375 nm inside the trap. The
ions were Doppler laser-cooled on the
$(4\text{s})~{}^{2}\text{S}_{1/2}\leftrightarrow(4\text{p})~{}^{2}\text{P}_{1/2}\leftrightarrow(3\text{d})~{}^{2}\text{D}_{3/2}$
system of optical cycling transitions using diode laser beams at 397 nm and
866 nm [11]. The laser beams were inserted at an angle of $\approx 45^{\circ}$
to all principal trap axes (see Fig. 1) to ensure cooling in all spatial
directions.

### Vibrational spectrum of the nanowire

Oscillations of the nanowire (NaugaNeedles) were measured by placing its tip
in the focus of a laser beam (here, a part of the 866 nm ion cooling beam was
used) and imaging its shadow onto a dual-element Si PIN photodiode [46]. The
amplitude of the vibrations was strongest at the oscillator’s resonance
frequencies. The driven motion of the oscillator lead to modulations of the
difference signal from the two elements of the photodiode. The signal was
demodulated at the drive frequency $f_{\text{drive}}$ using a lock-in
amplifier (Zurich Instruments HF2LI) yielding the mechanical excitation
spectrum displayed in Fig. 2a. The observed features correspond to combined
excitations of the nanowire on its tungsten holder. The measured
eigenfrequencies were found to be consistent with simulations performed with
the structural mechanics module of COMSOL Multiphysics [48]. In these
simulations, the tungsten holder geometry was modeled as a cylindrical base
(diameter $d_{\text{holder}}=250$ µm and length $l_{\text{holder}}=4.15$ mm)
with a conical tip of length $l_{\text{tip}}=715$ µm. The Ag2Ga nanowire was
modeled as a cylinder of diameter $d_{\text{nw}}=150$ nm and length
$l_{\text{nw}}=16$ µm placed on the tip of the holder. The $n=4$ mode used in
the present experiments corresponds to a vibration with a predominant
oscillation amplitude of the nanowire and minor contributions of the support.

### Photon correlation method

A photon correlation method [43, 44, 45] was employed for the detection and
quantification of ion motion. This method was originally introduced for the
detection of excess ion micromotion in radiofrequency traps [43] and is used
here for the measurement of secular ion motion driven by the charged
mechanical oscillator. Briefly, laser-cooled ions undergoing oscillatory
motion in a trap exhibit periodic variations of their fluorescence because of
a continuously varying detuning of the laser frequency with respect to the
cooling resonance in the frame of the ion due to the Doppler effect. The
photon scattering rate $R(t)$ is thus modulated according to

$R(t)=R_{0}\frac{(\Gamma_{12}/2)^{2}}{(\Gamma_{12}/2)^{2}+(\delta_{0}-kv_{\text{max}}\cos{(\omega
t-\varphi)})^{2}}$ (1)

where $R_{0}$ is the scattering rate on resonance, $\Gamma_{12}$ the natural
linewidth of the transition, $\delta_{0}$ the detuning of the cooling-laser
frequency from resonance for an ion at rest, $k$ is the projection of the
laser’s wave vector on the direction of motion of the ion, $v_{\text{max}}$ is
the maximum velocity of the ion during vibration, $\omega$ is the oscillation
frequency, $t$ is the time and $\varphi$ is a phase determined by the initial
conditions of the experiment. The arrival times of the fluorescence photons on
the PMT were measured with a time-to-amplitude converter (TAC, Stanford
Research Systems SR620) in correlation with the drive of the nanowire yielding
the time-resolved fluorescence curves displayed in Fig. 2b. A fit of the data
to Eq. 1 yields the velocity amplitudes $v_{\text{max}}$ imparted to the ion
by the mechanical drive of the nanowire. Uncertainties of $v_{\text{max}}$
quoted include both fit and statistical errors of the measurement.

### Theoretical modelling

The ion-nanowire hybrid system was modelled using classical-dynamics
simulations following the treatment of Ref. [28]. The nanowire was modeled as
a single point charge at position
$\vec{x}_{\text{nw}}=(x_{\text{nw}},y_{\text{nw}},z_{\text{nw}})$,
corresponding to the leading term of the multipole expansion of the
interaction of its asymmetric charge distribution with the ion. The total
electric potential $\Phi_{\text{tot}}$ experienced by an ion at position
$\vec{x}=(x,y,z)$ was thus described by the sum of the trapping potential
$\Phi_{\text{trap}}$ and a Coulomb interaction:

$\Phi_{\text{tot}}=\Phi_{\text{trap}}+k_{\text{c}}\frac{q_{\text{nw}}q_{\text{ion}}}{\lvert\vec{r}\rvert}$
(2)

Here,
$\lvert\vec{r}\rvert=\sqrt{(x-x_{\text{nw}})^{2}+(y-y_{\text{nw}})^{2}+(z-z_{\text{nw}})^{2}}$
is the effective ion-nanowire distance, $q_{\text{nw}}$ and $q_{\text{ion}}$
are the charges of the nanowire and the ion and $k_{\text{c}}$ is the Coulomb
constant. The mechanical drive of the nanowire was added by substituting
$z_{\text{nw}}\rightarrow z_{\text{nw}}+A\cos{(\omega_{\text{nw}}t)}$, where
$A$ is the amplitude and $\omega_{\text{nw}}$ the angular frequency of the
nanowire oscillation. The nanowire was assumed to vibrate in the $z$-direction
in line with the direction of the piezo drive in the experiment.
$\Phi_{\text{trap}}$ was expressed by a time-independent harmonic
pseudopotential [1].
Taylor expansion of the interaction potential to second order yielded an
approximate simulation model with equations of motions corresponding to a
periodically driven harmonic oscillator (see also [52]). From analysis of the
Taylor-expanded potential, $v_{\text{max}}$ behaves as a function of $A$ and
$d$ as

$v_{\text{max}}\propto(\frac{1}{d^{3}}-3\frac{z_{\text{nw}}^{2}}{d^{5}})A$ (3)

where $d$ is the distance between the ion-nanowire equilibrium positions. This
expression was used to fit the data in Figs. 4a,b.

The equations of motion were solved numerically using the Velocity Verlet
algorithm [53]. Typical integration time steps were $\Delta t=1$ ns with
$1\times 10^{6}$ steps leading to a total simulation time of 1 ms. The
nanowire oscillation frequency was matched to the axial trap frequency
$\omega_{\text{nw}}=\omega_{\text{z}}$ as in the experiments. The nanowire’s
amplitude response $A(f_{\text{drive}})$ to the drive frequency
$f_{\text{drive}}$ was modeled by a Lorentzian function [42]:

$A(f_{\text{drive}})=A_{0}\frac{f_{\text{R}}^{2}/Q}{\sqrt{(f_{\text{R}}^{2}-f_{\text{drive}}^{2})^{2}+f_{\text{R}}^{2}f_{\text{drive}}^{2}/Q^{2}}}$
(4)

where $f_{\text{R}}=422$ kHz is the $n=4$ resonance frequency and $A_{0}$ is
the resonant amplitude response at $f_{\text{drive}}=f_{\text{R}}$. The
quality factor $Q=f_{\text{R}}/\Delta f=117$ was estimated from the full-
width-at-half-maximum (FWHM) $\Delta f=3.6$ kHz of the 422 kHz resonance peak
(see fig. 3).

## References

  * Major  _et al._ [2005] F. G. Major, V. N. Gheorghe, and G. Werth,  _Charged Particle Traps_ (Springer, Berlin & Heidelberg, 2005).
  * Leibfried  _et al._ [2003] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Quantum dynamics of single trapped ions, Rev. Mod. Phys. 75, 281 (2003).
  * Wang  _et al._ [2021] P. Wang, C. Luan, M. Qiao, M. Um, J. Zhang, Y. Wang, X. Yuan, M. Gu, J. Zhang, and K. Kim, Single ion qubit with estimated coherence time exceeding one hour, Nat. Comm. 12, 233 (2021).
  * Häffner  _et al._ [2008] H. Häffner, C. F. Roos, and R. Blatt, Quantum computing with trapped ions, Phys. Rep. 469, 155 (2008).
  * Flühmann  _et al._ [2019] C. Flühmann, T. L. Nguyen, M. Marinelli, V. Negnevitsky, K. Mehta, and J. P. Home, Encoding a qubit in a trapped-ion mechanical oscillator, Nature 566, 513 (2019).
  * Pino  _et al._ [2021] J. M. Pino, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, M. S. Allman, C. H. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, C. Ryan-Anderson, and B. Neyenhuis, Demonstration of the trapped-ion quantum CCD computer architecture, Nature 592, 209 (2021).
  * Lanyon  _et al._ [2011] B. P. Lanyon, C. Hempel, D. Nigg, M. Müller, R. Gerritsma, F. Zähringer, P. Schindler, J. T. Barreiro, M. Rambach, G. Kirchmair, M. Hennrich, P. Zoller, R. Blatt, and C. F. Roos, Universal digital quantum simulation with trapped ions, Science 334, 57 (2011).
  * Blatt and Roos [2012] R. Blatt and C. F. Roos, Quantum simulations with trapped ions, Nat. Phys. 8, 277 (2012).
  * Whitlow  _et al._ [2023] J. Whitlow, Z. Jia, Y. Wang, C. Fang, J. Kim, and K. R. Brown, Quantum simulation of conical intersections using trapped ions, Nature Chemistry 15, 1509 (2023).
  * Ludlow  _et al._ [2015] A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, Optical atomic clocks, Rev. Mod. Phys. 87, 637 (2015).
  * Willitsch [2012] S. Willitsch, Coulomb-crystallised molecular ions in traps: methods; applications; prospects, Int. Rev. Phys. Chem. 31, 175 (2012).
  * Meir  _et al._ [2016] Z. Meir, T. Sikorsky, R. Ben-shlomi, N. Akerman, Y. Dallal, and R. Ozeri, Dynamics of a Ground-State Cooled Ion Colliding with Ultracold Atoms, Phys. Rev. Lett. 117, 243401 (2016).
  * Willitsch [2017] S. Willitsch, Chemistry with controlled ions, Adv. Chem. Phys. 162, 307 (2017).
  * Willitsch [2015] S. Willitsch, Ion-atom hybrid systems, Proc. Int. Sch. Phys. Enrico Fermi 189, 255 (2015).
  * Tomza  _et al._ [2019] M. Tomza, K. Jachymski, R. Gerritsma, A. Negretti, T. Calarco, Z. Idziaszek, and P. S. Julienne, Cold hybrid ion-atom systems, Rev. Mod. Phys. 91, 035001 (2019).
  * Meekhof  _et al._ [1996] D. M. Meekhof, C. Monroe, B. E. King, W. M. Itano, and D. J. Wineland, Generation of nonclassical motional states of a trapped atom, Phys. Rev. Lett. 76, 1796 (1996).
  * Srinivas  _et al._ [2023] R. Srinivas, C. M. Löschnauer, M. Malinowski, A. C. Hughes, R. Nourshargh, V. Negnevitsky, D. T. C. Allcock, S. A. King, C. Matthiesen, T. P. Harty, and C. J. Ballance, Coherent control of trapped-ion qubits with localized electric fields, Phys. Rev. Lett. 131, 020601 (2023).
  * Johanning  _et al._ [2009] M. Johanning, A. Braun, N. Timoney, V. Elman, W. Neuhauser, and C. Wunderlich, Individual addressing of trapped ions and coupling of motional and spin states using rf radiation, Phys. Rev. Lett. 102, 073004 (2009).
  * Ospelkaus  _et al._ [2011] C. Ospelkaus, U. Warring, Y. Colombe, K. R. Brown, J. M. Amini, D. Leibfried, and D. J. Wineland, Microwave quantum logic gates for trapped ions, Nature 476, 181 (2011).
  * Brown  _et al._ [2011] K. R. Brown, C. Ospelkaus, Y. Colombe, A. C. Wilson, D. Leibfried, and D. J. Wineland, Coupled quantized mechanical oscillators, Nature 471, 196 (2011).
  * Harlander  _et al._ [2011] M. Harlander, R. Lechner, M. Brownnutt, R. Blatt, and W. Hänsel, Trapped-ion antennae for the transmission of quantum information, Nature 471, 200 (2011).
  * An  _et al._ [2022] D. An, A. M. Alonso, C. Matthiesen, and H. Häffner, Coupling two laser-cooled ions via a room-temperature conductor, Phys. Rev. Lett. 128, 063201 (2022).
  * Tian and Zoller [2004] L. Tian and P. Zoller, Coupled ion-nanomechanical systems, Phys. Rev. Lett. 93, 266403 (2004).
  * Hensinger  _et al._ [2005] W. K. Hensinger, D. W. Utami, H.-S. Goan, K. Schwab, C. Monroe, and G. J. Milburn, Ion trap transducers for quantum electromechanical oscillators, Phys. Rev. A 72, 041405 (2005).
  * Liu  _et al._ [2007] T. Liu, M. Feng, and K.-L. Wang, Exact solution of quantum dynamics of a cantilever coupling to a single trapped ultracold ion, Commun. Theor. Phys. 47, 561 (2007).
  * Nicacio  _et al._ [2013] F. Nicacio, K. Furuya, and F. L. Semião, Motional entanglement with trapped ions and a nanomechanical resonator, Phys. Rev. A 88, 022330 (2013).
  * Kotler  _et al._ [2017] S. Kotler, R. W. Simmonds, D. Leibfied, and D. J. Wineland, Hybrid quantum systems with trapped charged particles, Phys. Rev. A 95, 022327 (2017).
  * Fountas  _et al._ [2019] P. Fountas, M. Poggio, and S. Willitsch, Classical and quantum dynamics of a trapped ion coupled to a charged nanowire, New J. Phys 21, 013030 (2019).
  * Daniilidis and Häffner [2013] N. Daniilidis and H. Häffner, Quantum interfaces between atomic and solid-state systems, Annu. Rev. Condens. Matter Phys. 4, 83 (2013).
  * Kurizki  _et al._ [2015] G. Kurizki, P. Bertet, Y. Kubo, K. Molmer, D. Petrosyan, P. Rabl, and J. Schmiedmayer, Quantum technologies with hybrid systems, Proc. Natl. Acad. Sci. 112, 3866 (2015).
  * Bachtold  _et al._ [2022] A. Bachtold, J. Moser, and M. I. Dykman, Mesoscopic physics of nanomechanical systems, Rev. Mod. Phys. 94, 045005 (2022).
  * O’Connell  _et al._ [2010] A. D. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bialczak, M. Lenander, E. Lucero, M. Neeley, D. Sank, H. Wang, M. Weides, J. Wenner, J. M. Martinis, and A. N. Cleland, Quantum ground state and single-phonon control of a mechanical resonator, Nature 464, 697 (2010).
  * Chan  _et al._ [2011] J. Chan, T. P. Mayer Alegre, A. H. Sfavi-Naeini, J. T. Hill, A. Krause, S. Gröblacher, M. Aspelmeyer, and O. Painter, Laser cooling of a nanomechanical oscillator into its quantum ground state, Nature 10, 1038 (2011).
  * Teufel  _et al._ [2011] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W. Simmonds, Sideband cooling of micromechanical motion to the quantum ground state, Nature 475, 359 (2011).
  * Rossi  _et al._ [2018] M. Rossi, D. Mason, J. Chen, Y. Tsaturyan, and A. Schliesser, Measurement-based quantum control of mechanical motion, Nature 563, 353 (2018).
  * Tebbenjohanns  _et al._ [2021] F. Tebbenjohanns, M. L. Mattana, M. Rossi, M. Frimmer, and L. Novotny, Quantum control of a nanoparticle optically levitated in cryogenic free space, Nature 595, 378 (2021).
  * Vochezer  _et al._ [2018] A. Vochezer, T. Kampschulte, K. Hammerer, and P. Treutlein, Light-mediated collective atomic motion in an optical lattice coupled to a membrane, Phys. Rev. Lett. 120, 073602 (2018).
  * Jöckel  _et al._ [2015] A. Jöckel, A. Faber, T. Kampschulte, M. Korppi, M. T. Rakher, and P. Treutlein, Sympathetic cooling of a membrane oscillator in a hybrid mechanical-atomic system, Nat. Nanotech. 10, 1038 (2015).
  * Wang and Lekavicius [2020] H. Wang and I. Lekavicius, Coupling spins to nanomechanical resonators: Toward quantum spin-mechanics, Appl. Phys. Lett 117, 230501 (2020).
  * Schmid  _et al._ [2022] G. Schmid, C. T. Ngai, M. Ernzer, M. Bosch Aguilera, T. M. Karg, and P. Treutlein, Coherent feedback cooling of a nanomechanical membrane with atomic spins, Phys. Rev. X 12, 011020 (2022).
  * Verhagen  _et al._ [2012] E. Verhagen, S. Deléglise, S. Weis, A. Schliesser, and T. J. Kippenberg, Quantum-coherent coupling of a mechanical oscillator to an optical cavity mode, Nature 482, 63 (2012).
  * Biedermann  _et al._ [2010] L. B. Biedermann, R. C. Tung, A. Raman, R. G. Reifenberger, M. M. Yazdanpanah, and C. R. W., Characterization of silver-gallium nanowires for force and mass sensing applications, Nanotechnology 21, 305701 (2010).
  * Berkeland  _et al._ [1998] D. J. Berkeland, J. D. Miller, J. C. Bergquist, W. M. Itano, and D. J. Wineland, Minimization of ion micromotion in a paul trap, J. Appl. Phys. 83, 5025 (1998).
  * Wang  _et al._ [2010] B. Wang, J. W. Zhang, Z. H. Lu, and L. J. Wang, Direct measurement of micromotion speed in a linear quadrupole trap, J. Appl. Phys. 108, 013108 (2010).
  * Keller  _et al._ [2015] J. Keller, H. L. Partner, T. Burgermeister, and T. E. Mehlstäubler, Precise determination of micromotion for trapped-ion optical clocks, J. Appl. Phys. 118, 104501 (2015).
  * Gloppe  _et al._ [2014] A. Gloppe, P. Verlot, E. Dupont-Ferrier, A. Siria, P. Poncharal, G. Bachelier, P. Vincent, and O. Arcizet, Bidimensional nano-optomechanics and topological backaction in a non-conservative radiation force field, Nat. Nanotech. 9, 920 (2014).
  * de Lépinay  _et al._ [2017] L. M. de Lépinay, B. Pigeau, B. Besga, P. Vincent, P. Poncharal, and O. Arcizet, A universal and ultrasensitive vectorial nanomechanical sensor for imaging 2d force fields, Nat. Nanotech. 12, 156 (2017).
  * [48] www.comsol.com, Comsol multiphysics version 3.5a and 4.3a, www.comsol.com.
  * Hauer  _et al._ [2013] B. D. Hauer, C. Doolin, K. S. D. Beach, and J. P. Davis, A general procedure for thermomechanical calibration of nano/micro-mechanical resonators, Ann. Phys. 339, 181 (2013).
  * Franke  _et al._ [2023] J. Franke, R. Muleady, S. R. Kaubruegger, F. Kranzl, R. Blatt, A. M. Rey, M. K. Joshi, and C. F. Roos, Quantum-enhanced sensing on optical transitions through finite-range interactions, Nature 621, 740 (2023).
  * Drewsen  _et al._ [2004] M. Drewsen, A. Mortensen, R. Martinussen, P. Staanum, and J. L. Sørensen, Nondestructive identification of cold and extremely localized single molecular ions, Phys. Rev. Lett. 93, 243201 (2004).
  * Hunger  _et al._ [2011] D. Hunger, S. Camerer, M. Korppi, A. Jöckel, T. W. Hänsch, and P. Treutlein, Coupling ultracold atoms to mechanical oscillators, C. R. Phys. 12, 871 (2011).
  * Spreiter and Walter [1999] Q. Spreiter and M. Walter, Classical molecular dynamics simulation with the velocity verlet algorithm at strong external magnetic fields, J. Comp. Phys. 152, 102 (1999).


# Endowments, Patience Types, and Uniqueness in Two-Good HARA Utility
Economies

Andrea Loi  Andrea Loi, Dipartimento di Matematica e Informatica
Università di Cagliari, Italy<EMAIL_ADDRESS>and  Stefano Matta  Stefano
Matta, Dipartimento di Scienze economiche e Aziendali
Università di Cagliari, Italy<EMAIL_ADDRESS>
###### Abstract.

This paper establishes a link between endowments, patience types, and the
parameters of the HARA Bernoulli utility function that ensure equilibrium
uniqueness in an economy with two goods and two impatience types with additive
separable preferences. We provide sufficient conditions that guarantee
uniqueness of equilibrium for any possible value of $\gamma$ in the HARA
utility function
$\frac{\gamma}{1-\gamma}\left(b+\frac{a}{\gamma}x\right)^{1-\gamma}$. The
analysis contributes to the literature on uniqueness in pure exchange
economies with two-goods and two agent types and extends the result in [4].

The first author was supported by INdAM. GNSAGA - Gruppo Nazionale per le
Strutture Algebriche, Geometriche e le loro Applicazioni and by KASBA, funded
by Regione Autonoma della Sardegna. Both authors were supported by STAGE,
funded by Fondazione di Sardegna.

Keywords: Uniqueness; Excess demand function; Additive separable preferences;
HARA utility; Polynomial approximation.

JEL Classification: C62, D51, D58.

##  1\. Introduction

The relationship between risk aversion, the number of consumer types $I$ and
the uniqueness of the equilibrium price has been analysed in a recent article
[4]. More precisely, in an economy with two goods and an arbitrary number $I$
of impatience types, where each type has additive separable preferences with a
HARA Bernoulli utility function
$u_{H}(x):=\frac{\gamma}{1-\gamma}\left(b+\frac{a}{\gamma}x\right)^{1-\gamma}$,
it has been shown that the equilibrium is unique if the parameter $\gamma$ is
in the range $\left(1,\frac{I}{I-1}\right]$.

While it is well known [2, 5, 6] the effect on uniqueness when $\gamma$ takes
a value between $0$ and $1$ (in the case $\gamma=1$ the function becomes
logarithmic), it is not known what conditions ensure uniqueness when $\gamma$
is greater than $2$. In this perspective, [4] analysed only the particular
case when $\gamma=3$ and $I=2$, and found sufficient conditions that ensure
uniqueness.

It is a natural question to ask whether a similar result can be found outside
the above interval. More specifically, this would mean finding sufficient
conditions that guarantee uniqueness of the equilibrium for any value of the
$\gamma$ parameter. This is related to [1]’s remarks on the difficulty of
finding a sufficient condition, expressible in closed form, that would allow,
for DARA-type utilities, the introduction of more heterogeneity among agents
in order to overcome the restrictive assumption of identical endowments to
ensure uniqueness (see [1, Proposition 2]).

This paper provides a positive answer to the question above for HARA
utilities, an important subclass of the DARA type. More precisely, our main
result, Theorem 2, shows the connection between endowments, patience types and
the parameters of the HARA utility function that ensure the uniqueness of the
equilibrium. To obtain this result, we will follow the approach of [4], where
the excess demand function is approximated by a polynomial whose variable, the
price, is raised to a power dependent on $\gamma$. An algebraic result, Lemma
4, which links the existence of a double root of a polynomial to an inequality
involving its coefficients, allows us to use a topological argument to prove
our main result.

For an overview of the literature on uniqueness, in addition to the well-known
contributions by [3] and [5], we also refer the reader to the two recent
contributions by [1] and [7] for two-good, two-agent pure exchange economies.

This short note is organised as follows. Section 2 analyses the economic
setting using the polynomial approach. Section 3 proves our main result.

##  2\. Preliminaries

Consider an economy with two goods and $I=2$ impatience types, where type $i$
has preferences represented by the utility function

(1) $u_{i}(x,y)=u_{H}(x)+\beta_{i}u_{H}(y),$

where $u_{H}$ is HARA, i.e.

(2)
$u_{H}(x):=\frac{\gamma}{1-\gamma}\left(b+\frac{a}{\gamma}x\right)^{1-\gamma},\
\gamma>0,\gamma\neq 1,a>0,b\geqslant 0.$

Let $\varepsilon$ be a rational number $\frac{m}{n}$, $m,n\in\mathbb{N}$
sufficiently close to $\frac{1}{\gamma}$. Suppose that $\gamma>2$ and, hence,
$n>2m$. Denoting by $(e_{i},f_{i})$ consumer i’s endowments, the standard
maximisation problem over the budget constraint $px_{i}+y_{i}\leqslant
pe_{i}+f_{i}$ gives (see [4, formula (14)]) the aggregate excess demand
function for good $x$:

(3)
$\sum_{i=1}^{2}\frac{b-bp^{\varepsilon}\sigma_{i}+a\varepsilon\left(pe_{i}+f_{i}\right)}{a\varepsilon\left(p+\sigma_{i}p^{\varepsilon}\right)}-(e_{1}+e_{2}),$

where

$\varepsilon\approx\frac{1}{\gamma},\ \sigma_{i}:=\beta_{i}^{\varepsilon},\
i=1,2.$

Following [4], we combine terms over a common denominator and take the
numerator, then we collect terms in $p$, divide by $p^{\varepsilon}$, and we
get:

$\begin{split}&p(-ae_{1}\sigma_{1}\varepsilon-
ae_{2}\sigma_{2}\varepsilon-b\sigma_{1}-b\sigma_{2})+p^{1-\varepsilon}(af_{1}\varepsilon+af_{2}\varepsilon+2b)+\\\
&p^{\varepsilon}(-ae_{1}\sigma_{1}\sigma_{2}\varepsilon-
ae_{2}\sigma_{1}\sigma_{2}\varepsilon-2b\sigma_{1}\sigma_{2})+af_{1}\sigma_{2}\varepsilon+af_{2}\sigma_{1}\varepsilon+b\sigma_{1}+b\sigma_{2}\end{split}$

Recalling that $\varepsilon=\frac{m}{n}$ and by letting, with a slight abuse
of notation, $x:=p^{1/n}$, we rewrite the previous expression in decreasing
order as follows:

(4)
$A(e,\sigma,a,b)x^{n}+B(f,\sigma,a,b)x^{n-m}+C(e,\sigma,a,b)x^{m}+D(f,\sigma,a,b),$

where

(5) $\displaystyle A(e,\sigma,a,b):=$
$\displaystyle-(e_{1}\sigma_{1}+e_{2}\sigma_{2})-\frac{b}{a\varepsilon}(\sigma_{1}+\sigma_{2})<0,$
$\displaystyle B(f,\sigma,a,b):=$
$\displaystyle(f_{1}+f_{2})+\frac{2b}{a\varepsilon}>0,$ $\displaystyle
C(e,\sigma,a,b):=$
$\displaystyle(e_{1}+e_{2})\sigma_{1}\sigma_{2}-\frac{2b}{a\varepsilon}\sigma_{1}\sigma_{2}<0,$
$\displaystyle D(f,\sigma,a,b):=$
$\displaystyle(f_{1}\sigma_{2}+f_{2}\sigma_{1})+\frac{b}{a\varepsilon}(\sigma_{1}+\sigma_{2})>0.$

######  Lemma 1.

If the following conditions hold

(6) $\beta_{1}<\beta_{2},\,e_{1}\leqslant e_{2},\,f_{1}\geqslant f_{2},$ (7)
$b\geqslant\frac{a}{\gamma}\left(\frac{\beta_{2}}{\beta_{1}}\right)^{\frac{2}{\gamma}}(e_{2}+f_{1}),$

then the polynomial (4) satisfies the inequality

(8) $A(e,\sigma,a,b)D(f,\sigma,a,b)-B(f,\sigma,a,b)C(e,\sigma,a,b)<0$

[Proof]

We will follow, mutatis mutandis, the same line of reasoning of the proof of
[4, Theorem 2], the only difference here is that we deal with an arbitrary
value of $\gamma$. The formula
$A(e,\sigma,a,b)D(f,\sigma,a,b)-B(f,\sigma,a,b)C(e,\sigma,a,b)$ can be written
as

$(\sigma_{2}-\sigma_{1})(e_{1}f_{2}\sigma_{1}-e_{2}f_{1}\sigma_{2})+E(e,f,\sigma,a,b),$

where

$E(e,f,\sigma,a,b):=-\frac{b^{2}}{a^{2}\varepsilon^{2}}\left(\sigma_{1}-\sigma_{2}\right)^{2}+\frac{b}{a\varepsilon}\left[\left(e_{1}+e_{2}+f_{1}+f_{2}\right)\sigma_{1}\sigma_{2}-(e_{1}+f_{2})\sigma_{1}^{2}-(e_{2}+f_{1})\sigma_{2}^{2}\right].$

Observe that, by condition (6),

$(\sigma_{2}-\sigma_{1})(e_{1}f_{2}\sigma_{1}-e_{2}f_{1}\sigma_{2})\leqslant(\sigma_{2}-\sigma_{1})f_{1}(e_{1}\sigma_{1}-e_{2}\sigma_{2})<(\sigma_{2}-\sigma_{1})f_{1}\sigma_{2}(e_{1}-e_{2})\leqslant
0.$

Moreover, $E(e,f,\sigma,a,b)\leqslant 0$ if and only if

$b\geqslant
a\varepsilon\frac{\left[(e_{1}+e_{2}+f_{1}+f_{2})\sigma_{1}\sigma_{2}-(e_{1}+f_{2})\sigma_{1}^{2}-(e_{2}+f_{1})\sigma_{2}^{2}\right]}{\left(\sigma_{1}-\sigma_{2}\right)^{2}}.$

Again, by (6), we can write

$\displaystyle
a\varepsilon\frac{\left[(e_{1}+e_{2}+f_{1}+f_{2})\sigma_{1}\sigma_{2}-(e_{1}+f_{2})\sigma_{1}^{2}-(e_{2}+f_{1})\sigma_{2}^{2}\right]}{\left(\sigma_{1}-\sigma_{2}\right)^{2}}$
$\displaystyle<a\varepsilon\frac{\left[(e_{1}+e_{2}+f_{1}+f_{2})\sigma_{1}\sigma_{2}\right]}{\left(\sigma_{1}-\sigma_{2}\right)^{2}}$
$\displaystyle<a\varepsilon\left(\frac{\sigma_{2}}{\sigma_{1}}\right)^{2}(e_{2}+f_{1}).$

Thus, since $\sigma_{i}=\beta_{i}^{\varepsilon}$, $i=1,2$ and
$\gamma=\frac{1}{\varepsilon}$, the proof of the lemma follows.

##  3\. Main result

In this section we present our main result, Theorem 2. As far as uniqueness is
concerned, we will assume an arbitrary $\gamma>2$. In fact, the case
$\gamma\in(1,2]$ is a particular case of [4, Theorem 1], while the case
$\gamma\leqslant 1$ is a well known result in the literature [2, 5, 6].

Observe that the zero set of aggregate demand function amounts to studying the
zeros of polynomial (4). In fact, according to [4]’s approach it is possible
to approximate $\gamma$ with a rational number, since $\mathbb{Q}$ is dense in
$\mathbb{R}$, in such a way that the cardinality of the set of regular
equilibria does not decrease [4, Lemma 9]. To provide a geometric insight,
this corresponds to small perturbations of the aggregate demand function that
do not allow a decrease in the number of the equilibria.

######  Theorem 2.

In an economy with two goods and two impatient types with Hara preferences
(1), if the conditions (6) and (7) hold, then the equilibrium price is unique.

######  Remark 3.

For the general type of DARA, [1] observe there is not a closed-form
expression that ensures uniqueness, but conditions (6) and (7) represent a
closed-form expression for HARA utilities, an important subclass of utilities
of type DARA. They are the same as those presented in [4], here suitably
generalised.

[Proof]By [4, Theorem 1]) there exists uniqueness of equilibrium if and only
if the polynomial (4), $P(x)$, has a unique positive root. We will prove that
the inequality (8), which holds by Lemma 1, implies that $P(x)$ has a unique
positive root. Assume by contradiction that $P(x)$ has more than one positive
root. Since $P(x)$ belongs to the path-connected space of polynomials
$Ax^{n}+Bx^{n-m}+Cx^{m}+D$, with non zero coefficients such that $AD-BC<0$, it
follows by the continuous dependence of the roots of a polynomial on its
coefficients that $P(x)$ has indeed a double positive root. Hence one can
achieve the conclusion of Theorem 2 by the following algebraic lemma.

######  Lemma 4.

If the polynomial (4), $P(x)={Ax^{n}+Bx^{n-m}+Cx^{m}+D}$, $ABCD\neq 0$, has a
double positive root, then $AD-BC\geqslant 0$.

[Proof]By contradiction, let $\alpha>0$ be a double root of $P(x)$, that is,
$(x-\alpha)^{2}$ divides $P(x)$. The following table shows the pattern of the
remainders after the first $k$ steps of the division.

Step |  Remainder
---|---
$1$ |  $2\alpha Ax^{n-1}-\alpha^{2}Ax^{n-2}+Bx^{n-m}+Cx^{m}+D$
$2$ |  $3\alpha^{2}Ax^{n-2}-2\alpha^{3}Ax^{n-3}+Bx^{n-m}+Cx^{m}+D$
$3$ |  $4\alpha^{3}Ax^{n-3}-3\alpha^{4}Ax^{n-4}+Bx^{n-m}+Cx^{m}+D$
$\vdots$ |  $\vdots$
$k$ |  $(k+1)\alpha^{k}Ax^{n-k}-k\alpha^{k+1}Ax^{n-k-1}+Bx^{n-m}+Cx^{m}+D$

From $n-k-1=n-m$, we get $m=k+1$ and then we can rewrite the reminder
accordingly:

$m\alpha^{m-1}Ax^{n-m+1}+[B-(m-1)\alpha^{m}A]x^{n-m}+Cx^{m}+D.$

Continuing the division with this new reminder, the next table reveals again
the following pattern:

Step |  Remainder
---|---
$1$ |  $[B+(m+1)\alpha^{m}A]x^{n-m}-m\alpha^{m+1}Ax^{n-m-1}+Cx^{m}+D$
$2$ |  $[2\alpha B+(m+2)\alpha^{m+1}A]x^{n-m-1}-[\alpha^{2}B+(m+1)\alpha^{m+2}A]x^{n-m-2}+Cx^{m}+D$
$3$ |  $[3\alpha^{2}B+(m+3)\alpha^{m+2}A]x^{n-m-2}-[2\alpha^{3}B+(m+2)\alpha^{m+3}A]x^{n-m-3}+Cx^{m}+D$
$\vdots$ |  $\vdots$
$k$ |  $k\alpha^{k-1}B+(m+k)\alpha^{m+k-1}A]x^{n-m-k+1}-[(k-1)\alpha^{k}B+(m+k-1)\alpha^{m+k}A)]x^{n-m-k}+Cx^{m}+D$

From $n-m-k=m$, we get $k=n-2m$. We can rewrite the reminder as follows:

$\begin{split}&[(n-2m)\alpha^{n-2m-1}B+(n-m)\alpha^{n-m-1}A]x^{m+1}-\\\
&-[(n-2m-1)\alpha^{n-2m}B+(n-m-1)\alpha^{n-m}A]x^{m}+Cx^{m}+D,\end{split}$

that, reordering terms, becomes

$\begin{split}&[(n-2m)\alpha^{n-2m-1}B+(n-m)\alpha^{n-m-1}A]x^{m+1}+\\\
&+[C-(n-2m-1)\alpha^{n-2m}B-(n-m-1)\alpha^{n-m}A]x^{m}+D.\end{split}$

Starting with this new reminder, the last pattern is suggested by the
following table:

Step |  Remainder
---|---
$1$ |  $[C+(n-2m+1)\alpha^{n-2m}B+(n-m+1)\alpha^{n-m}A]x^{m}-[(n-2m)\alpha^{n-2m+1}B+(n-m)\alpha^{n-m+1}A]x^{m-1}+D$
$2$ |  $[2\alpha C+(n-2m+2)]\alpha^{n-2m+1}B+(n-m+2)\alpha^{n-m+1}A]x^{m-1}-[\alpha^{2}C+(n-2m+1)\alpha^{n-2m+2}B+(n-m+1)\alpha^{n-m+2}A]x^{m-2}+D$
$3$ |  $[3\alpha^{2}C+(n-2m+3)]\alpha^{n-2m+2}B+(n-m+3)\alpha^{n-m+2}A]x^{m-2}-[2\alpha^{3}C+(n-2m+2)\alpha^{n-2m+3}B+(n-m+2)\alpha^{n-m+3}A]x^{m-2}+D$
$\vdots$ |  $\vdots$
$k$ |  $[k\alpha^{k-1}C+(n-2m+k)]\alpha^{n-2m+k-1}B+(n-m+k)\alpha^{n-m+k-1}A]x^{m-k+1}-[(k-1)\alpha^{k}C+(n-2m+k-1)\alpha^{n-2m+k}B+(n-m+k-1)\alpha^{n-m+k}A]x^{m-k}+D$

After $k=m$ divisions, the remainder reduces to a first degree polynomial:

$m\alpha^{m-1}C+(n-m)\alpha^{n-m-1}B+n\alpha^{n-1}A]x-[(m-1)\alpha^{m}C+(n-m-1)\alpha^{n-m}B+(n-1)\alpha^{n}A]+D.$

Under the hypothesis that $(x-\alpha)^{2}$ divides $P(x)$, the coefficients of
this last remainder must vanish, that is:

$\begin{cases}m\alpha^{m-1}C=-(n-m)\alpha^{n-m-1}B-n\alpha^{n-1}A\\\
D=(m-1)\alpha^{m}C+(n-m-1)\alpha^{n-m}B+(n-1)\alpha^{n}A.\\\ \end{cases}$

Multiplying the second equation by $m\alpha^{m-1}$, we get

$m\alpha^{m-1}D=(m-1)\alpha^{m}m\alpha^{m-1}C+m(n-m-1)\alpha^{n-1}B+m(n-1)\alpha^{n+m-1}A,$

where, substituting $m\alpha^{m-1}C$ with the RHS of the first equation and
multiplying by $A$, we obtain

$m\alpha^{m-1}AD=(n-2m)\alpha^{n-1}AB+(n-m)\alpha^{n+m-1}A^{2}.$

Moreover, we observe that

$m\alpha^{m-1}BC=-(n-m)\alpha^{n-m-1}B^{2}-n\alpha^{n-1}AB.$

We can then write

$m\alpha^{m-1}(AD-
BC)=(n-m)\alpha^{n-m-1}(\alpha^{2m}A^{2}+B^{2})+2(n-m)\alpha^{n-1}AB,$

as

$(n-m)\alpha^{n-m-1}(\alpha^{2m}A^{2}+B^{2}+2\alpha^{m}AB),$

or, equivalently,

$(n-m)\alpha^{n-m-1}(\alpha^{m}A+B)^{2}.$

Hence, we have

$AD-BC=\frac{n-m}{m}\alpha^{n}(\alpha^{n}A+B)^{2}\geqslant 0,$

yielding the desired contradiction.

######  Remark 5.

It should be possible to give a more elegant proof of Lemma 4 by using an
approach based on the discriminant of a polynomial instead of a division
algorithm as in our proof. However, this alternative approach seems to lead to
very complicated calculations that the authors were unable to handle.

## References

  * Geanakoplos and Walsh [2018] Geanakoplos, J., Walsh, K. J., 2018, Uniqueness and stability of equilibrium in economies with two goods, Journal of Economic Theory, 174, 261-272.
  * Hens and Loeffler [1995] Hens, T., Loeffler, A., 1995, Gross substitution in financial markets, Econ. Lett. 49(1), 39- 43.
  * Kehoe [1998] Kehoe, Timothy J., 1998, Uniqueness and stability. In: Kirman, Alan (Ed.), Elements of General Equilibrium Analysis. Wiley Blackwell, pp. 38-87 (Chapter 3).
  * Loi and Matta [2022] Loi, A., Matta, S., 2022, Risk Aversion and Uniqueness of Equilibrium in Economies with Two Goods and Arbitrary Endowments, The B.E. Journal of Theoretical Economics, https://doi.org/10.1515/bejte-2021-0150.
  * Mas-Colell, Andreu, [1991] Mas-Colell, Andreu, 1991, On the uniqueness of equilibrium once again. In: Barnett, William A., Cornet, Bernard, D’Aspermont, Claude, Gabszewicz, Jean, Mas-Colell, Andreu (Eds.), Equilibrium Theory and Applications. Cambridge University Press, pp. 275-296 (Chapter 12).
  * Mas-Colell, A. et al., [1995] Mas-Colell, Andreu, Whinston, Michael D., Green, Jerry R., 1995, Microeconomic Theory. Oxford University Press.
  * Toda and Walsh [2017] Toda, A. A., Walsh, K. J., 2017, Edgeworth box economies with multiple equilibria, Econ. Theory Bull., 5, 65-80.