\section{Introduction}
Ever since the emergence of modern financial markets in the 17th century, specialists have made an effort to maximize profits from trading goods, services, stocks and securities. Consequently, consistent alpha generation \cite{gross2005consistent}, a term to describe the edge over the market or the ability to make profit, has remained a hot topic \cite{fung2002news, saad1998comparative}. With the advent of the digital age and the ability to process large amounts of data, many strategies have revolved around new techniques and computational breakthroughs in an effort to predict the upward and downward trends of assets (e.g., \cite{lee2009using, ni2011stock}). The reliability of these predictions often remains questionable \cite{fama1995random, vui2013review}. This difficulty in predicting traditional asset data is all the more true for the recent digital asset (e.g. Bitcoin) markets, which have proven to be extraordinarily volatile \cite{ali2014economics, dyhrberg2016bitcoin, yermack2015bitcoin}.

However, the recent rise of AI and ML has opened up new opportunities. AI techniques can analyze enormous amounts of data and observe decades of conventional experience at a rate that is only limited by the complexity of the technique and the available computational power \cite{bottou2010large, bottou2018optimization}. Furthermore, modern machine learning algorithms are capable of discovering patterns in data that are hidden to the human eye due to the volume and complexity of the data \cite{witten2016data}. Unsurprisingly, investment and asset management firms have started to invest billions into developing AI-based asset trading and investment applications \cite{bloomberg2019quants}.

In this paper, we present our findings and experiences in designing and implementing \emph{RCURRENCY}, an automated trading engine tasked with managing a portfolio of digital assets (Section \ref{sec:model}). The primary goal of our research is to provide insights into the capabilities of AI-powered digital asset trading in a highly volatile and therefore challenging live market environment. Effectively managing a portfolio requires the AI-based training engine to estimate when to sell off assets, buy new ones, or hold the current position. Decisions are made based on data gathered from a digital asset exchange at a regular interval.
For this purpose, we combine the prediction with previous (known) values and apply a trading strategy to derive and execute a final buy/hold/sell decision. In practice, achieving a stable portfolio is considered a high bar for automatic trading engines to master.

We experimentally evaluate different trading strategies based on actual historical market data (Section \ref{sec:experiments}). Our results show that the prediction accuracy is highly dependent on the configuration of the neural network. The Root Mean Squared Error over a prediction can increase significantly when incorrect parameters are chosen.
We thoroughly evaluate the most important parameter spaces and their impact on the prediction accuracy. When an adequate set of parameters is chosen, the trading engine is capable of estimating future asset values with an error of approximately 0.4\%. In practice, this makes it possible to reliably predict up and down motions, as well as approximate the absolute change in value. This, in turn, allows the system to maintain and increase portfolio value over time.

\section{Related Work}

The suggestion of predicting stock prices utilizing historic data has been an active research topic for decades. One of the early prediction methods based on an inter-temporal general equilibrium model was proposed by \citet{balvers1990predicting}.  These early methods used a purely analytical approach applying basic mathematics and statistics. However, some schools of thought even then disputed the merit of the approach and results of predicting the stock prices. For instance, it has been empirically shown that if a random walk strategy accurately reflects the reality, then all models concerning (stock) market predictions are meaningless~\cite{fama1995random}.

Nowadays, more variables are added to the prediction model, such as the market volume \cite{brooks1998predicting}, twitter sentiment \cite{zhang2011predicting}, investor sentiment, or media news content \cite{tetlock2007giving}.
Whereas fundamental analysis describes the examination of data focused on macro-economic variables \cite{abarbanell1997fundamental}, technical analysis puts more emphasis on extracting trends and patterns of an investment instrument's price, volume, breadth and trading activities \cite{thomsett1999mastering,pring1991technical,meyers1994technical}.
The methodology proposed and analyzed in this paper focuses on technical analysis while it could be expanded to support fundamental analysis in future works. Most recently, multiple studies have used AI techniques for prediction. Examples of concrete techniques include particle swarm optimization \cite{nenortaite2004stocks}, nearest neighbor classification \cite{teixeira2010method}, and neural networks \cite{lam2004neural}. An empirical study by \citet{gu2018empirical} compared multiple AI techniques in combination with trading techniques in the conventional stock market. Their fundamental result is that overall AI can be successfully utilized and the predictive capability of the AI techniques carries a significant advantage. However, using AI in highly volatile markets such as digital assets \cite{dyhrberg2016bitcoin} has gone largely unexplored to this date.

Some initial efforts to predict prices and trends of Bitcoin and other digital currencies have been undertaken by \citet{mcnally2018predicting} and \cite{alessandretti2018anticipating}. \citet{mcnally2018predicting} explore prediction of Bitcoin prices based on the Simple Moving Average (SMA) technical indicator, through RNN, LSTM and ARIMA methods.

\citet{alessandretti2018anticipating} investigate price prediction of a large number of cryptocurrencies using tree boosting techniques, as well as an LSTM-based method. The researchers simulate trading of currencies, taking into account transaction fees, but fail to explore the myriad of strategies traders use to determine their investments. In either case the number of input variables is limited; resulting the loss of potentially valuable network input data.

\section{The RCURRENCY Live Trading System} \label{sec:model}

In order to address the challenge of making reliable predictions of cryptocurrency markets for live trading, we designed the RCURRENCY system and structured it into three separate components. The responsibility of the first component is to collect and preprocess the timely-ordered series of data. The second component feeds this data into the neural network and, by doing so, trains it to recognize patterns in the processed data.
These two components are sufficient to grant the system predictive capabilities for single steps in the time series data. The third component allows the system to make trading decisions depending not only on a single prediction, but also on an elaborate analysis of trends and patterns. It implements and applies various trading strategies, each of which calculates a trading signal based on its respective predefined rules. The flow diagram of the RCURRENCY system is visualized in Figure~\ref{fig:rcurrency_flow}.

\begin{figure*}
    \centering
    \includegraphics[width=0.82\textwidth]{images/flow.png}
    \caption{Structure and data flow in RCURRENCY.}
    \label{fig:rcurrency_flow}
\end{figure*}

\subsection{Training Data}

A neural network-based prediction method is only as reliable as the data that can be leveraged to train the network. As a result, an essential step in the processing pipeline is to collect the raw data and turn it into preprocessed training data for the learning algorithm to digest. Our system obtains raw Bitcoin training data by calling an API provided by CryptoCompare\footnote{https://min-api.cryptocompare.com/}. This company offers services concerning over 3.000 digital assets, including the ability to retrieve historical price data. Such data comprises candlestick data points in one-minute time frames, where each data point contains four values: the \emph{open}, \emph{high}, \emph{low}, and \emph{close} price of the time frame. The system runs a simple script to recursively and sequentially call this API to fetch historical Bitcoin price data, which consists of approximately 80.000 hourly data points, spanning back to August 2010. In contrast, live data on the current market condition is instead fetched directly from a cryptocurrency exchange to minimize acquisition delays.

The raw data then undergoes preprocessing. Unfortunately, approximately the first 30 percent of the data set contains many near-zero values that tend to deteriorate the prediction capabilities of the network. These values do not contain sufficient variance for the gradient descent of the error function to converge. As a consequence, all data up until May 2013 has been manually removed, roughly up until the time Bitcoin reached values around one hundred USD.

Next, the system applies some technical analysis on the resulting data set, which involves extracting trends and patterns from the price data by constructing technical trading indicators. Each indicator is calculated according to their respective formula, where each formula uses some (or all) of the data set as input. RCURRENCY implements the financial technical indicators.
The calculated indicator values are added as additional rows to the training data set. This extra training data provides the neural network with more information about the market and should thus lead to higher predictive accuracy.

The data is subsequently converted to what is called a \emph{stationary data set}. Such a data set is created by taking the difference of each data point from the previous data point in a column-based manner.

Finally, row-based normalization is applied to the data set to fit the resulting values in the range of -1 and 1 in order to obtain better results and reduce the training time \cite{sola1997importance}.
The new row-normalized matrix is mathematically defined as follows:

Given a matrix $x$, with dimensions $m \times n$, for the new matrix $y$ holds:
\begin{equation}
y_{i,j} = 2 * (\dfrac{x_{i,j} - y^{min}}{y^{max} - y^{min}}) - 1, 1 \leq i \leq m, \& 1 \leq j \leq n.
\end{equation}

\subsection{Recurrent Neural Network}

There are various options available regarding networks capable of consuming a series of data points in order to make predictions on future values, such as the time series analysis technique ARIMA \cite{mills1991time}, or neural networks. As evaluated by \citet{kohzadi1996comparison}, a neural network outperforms the ARIMA network over longer periods of time, justifying its preference for our system. Specifically a recurrent neural network functions as RCURRENCY's predictor, which is a particular type of neural network capable of exhibiting temporal dynamic behavior. As first introduced by \cite{mikolov2010recurrent}, these networks model the connections between layer nodes as a directed graph along a temporal sequence. The input data for RCURRENCY's recurrent neural network consists of so-called candlesticks, data structures containing numerical stock market information on assets. The values in a candlestick are dependent on previous candlestick values and do therefore not follow a random walk \cite{lo1988stock}.

A conventional recurrent neural network experiences problems with the gradient descent when learning long-time dependencies, as depicted in the research by \citet{bengio1994learning}. When an LSTM layer is used in the architecture of the recurrent neural network, the problems of exploding and vanishing gradients are solved \cite{Hochreiter1997}. An example of a long-term dependency in the current context is given by \citet{tsinaslanidis2014prediction}. Their paper shows that a previously encountered subsequence of the time series data could be similar to the current trend. The consequent actions of the price following the previously encountered subsequence could carry value in the predictions of the actions following the current subsequence.

The recurrent neural network encompasses three layers:

\begin{enumerate}
    \item \textit{Linear layer}\\
    The linear layer is required to map the nodes of the input layer to the number of nodes in the hidden layer. It is a fully connected layer, of which every edge has an individual weight.

    \item \textit{Fast LSTM layer}\\
    The FastLSTM layer is a special, faster version of the conventional LSTM layer. It combines the calculation of the input, forget and output gates and the hidden state in a single step.

The composite functions of the standard LSTM cell are based on the paper from \cite{graves2013speech}.
and are defined as follows:
    \begin{equation}
        \begin{split}
            i_t = \sigma(W_{xi}x_t + W_{hi}h_{t-1} + W_{ci}c_{t-1}+b_i)\\
            f_t = \sigma(W_{xf}x_t + W_{hf}h_{t-1} + W_{cf}c_{t-1}+b_f)\\
            c_t = f_tc_{t-1} + i_t\text{tanh}(W_{xc}x_t + W_{hc}h_{t-1} + b_c)\\
            o_t = \sigma(W_{xo}x_t + W_{ho}h_{t-1} + W_{co}c_{t-1}+b_o)\\
            h_t = o_t\text{tanh}(c_t).
        \end{split}
    \end{equation}
    However, as the network implements a FastLSTM layer rather than a conventional LSTM layer, the composite functions are slightly altered \cite{mlpack_fastlstm}.
    \begin{equation}
        \begin{split}
            i_t = \sigma(W_{xi}x_t + W_{hi}h_{t-1} + b_i)\\
            f_t = \sigma(W_{xf}x_t + W_{hf}h_{t-1} + b_f)\\
            c_t = f_tc_{t-1} + \text{tanh}(W_{xc}x_t + W_{hc}h_{t-1} + b_c)\\
            o_t = \sigma(W_{xo}x_t + W_{ho}h_{t-1} + b_o)\\
            h_t = \text{tanh}(c_t).
        \end{split}
    \end{equation}
    \item \textit{Output layer with Mean Squared Error (MSE) performance function}\\

    The output layer evaluates the neural network, in case it is used for training purposes. If the neural network is used to predict a values, the value that is returned reflects the output of the specified output layer. In this case, the MSE function is used to evaluate the network.
    It does that as follows:
    \begin{equation}
     MSE = \dfrac{1}{m * n} \sum\limits_{i=1}^m \sum\limits_{j=1}^n (x_{ij} - \hat{x}_{ij})^2
    \end{equation}
    where $x_{ij}$ is the actual value, $\hat{x}_{ij}$ is the predicted value in row $i$ and column $j$, $m$ is the total amount of rows, and $n$ is the total amount of columns.
\end{enumerate}

\noindent The network also uses the Adam optimizer, introduced by \citet{kingma2014adam}, to speed up the training process. The Adam optimizer is computationally efficient, requires little memory, is invariant to diagonal re-scaling of the gradients, and is well suitable for large quantities of data, or problems with many parameters. Therefore, it is a good fit for the solution.
The resulting formula for updating parameters is:
\begin{equation}
    \begin{split}
        \hat{m}_t = \dfrac{m_t}{1-\beta_1^t},\\
        \hat{v_t} = \dfrac{v_t}{1-\beta_2^t},\\
        a_t = \dfrac{a*\sqrt{1-\beta_2^t}}{1-\beta_1^t},\\
        \theta_{t+1} = \theta_t - \dfrac{a_t * m_t}{\sqrt{\hat{v}_t} + \epsilon}\hat{m}_t,
    \end{split}
\end{equation}
where $\alpha$ is the step size, $m_t$ is the mean at time $t$, and $v_t$ is the uncentered variance of the gradients at time $t$.

\subsection{Trading Strategies and Prediction}

Technical analysis of digital assets is indeed sensitive to the trading strategy, i.e, the methods of buying and selling based on predefined rules to make trading decisions.
There have been longstanding discourses on how useful trading strategies are in practice \cite{knez1996estimating}. Common trading strategies based on time series pattern analysis, such as momentum and contrarian, as well as other technical analysis strategies have had varying degrees of success \cite{conrad1998anatomy, bessembinder1998market, teixeira2010method}. As such, it was decided to add this element to the experimental AI, giving it the capability of applying four different trading strategies based on technical indicators. Due to the vast number of strategies available and conceivable, priority was given to those that are most commonly adopted by existing trading platforms, such as TradingView\footnote{https://www.tradingview.com/} or Cryptowatch\footnote{https://cryptowat.ch/}. Each trading strategy is briefly discussed below.

Every time period (minutely, hourly, daily) the neural network makes a prediction of the next period's asset price. This value is digested by the chosen trading strategy or strategies to determine a decision signal, which is most commonly facilitated by observing the intersection between two or more moving averages. Volatility of the portfolio value can be contained by setting limits on the minimum and maximum amount of portfolio value that can be traded each time period. Altering these limits and changing the trading strategies determines the level of risk management in trading.

\begin{itemize}

\item Rate-of-Change (ROC)\\
The ROC strategy considers the change of the predicted value relative to the last actual value. When the change is below, in between, or above a lower and upper threshold, a decision is made whether to sell, buy or hold.

\item Relative Strength Indicator (RSI)\\
This strategy tracks upward and downward trends in asset value. A sufficiently long series of positive value changes result in an asset being marked as overbought (indicating the holder should sell) and vice versa.

\item Double Exponential Moving Averages (DEMA)\\
The DEMA strategy utilizes two moving averages to determine whether to buy, hold or sell an asset. Crossovers of these averages indicate buy and sell moments as the momentum changes between long-term and short-term asset value.

\item Moving Average Convergence / Divergence (MACD)\\
MACD is a well-known strategy based on the difference between two moving averages with different lengths, called \emph{slow-} and \emph{fast moving average}. This is then compared to the moving average of the difference itself, called the \emph{signal line}. A decision to buy, hold or sell is made whenever the difference crosses the signal line.

\item Random Walk\\
A controversial yet widely regarded investment theory known as the efficient market hypothesis states that all available information is immediately reflected in the value of an asset, and therefore no consistently profitable trading strategies should exist \cite{basu1977investment}. Consequently, a series of random actions should be able to approximate the results of a series of calculated trades.
\end{itemize}

\subsection{Implementation}
RCURRENCY is implemented in C++ with various libraries making up its core.
First, the historical data is fetched using the API provided by CryptoCompare. Since the API is currently not rate-limited or metered, this step takes only minutes. Live data is fetched in real-time by creating a WebSocket connection with the exchange (Binance) and fetching the candles after they have formed.

After obtaining the historical data, it is preprocessed and fed to the neural network. This process utilizes two well-known libraries: Armadillo\footnote{http://arma.sourceforge.net/}
and TA-Lib\footnote{http://ta-lib.org/}.
The Armadillo library provides highly optimized matrix computation functionality desired for regular data set manipulation. TA-Lib provides all necessary functionality for efficiently calculating technical indicators from stock market data. As such, preprocessing and preparation of the training dataset is completed in under a minute.

The recurrent neural network, as described in Section \ref{sec:model}, is implemented using MLPack \cite{mlpack2018}. MLPack is a fast and flexible C++ machine learning library, implementing various types of configurable neural networks and machine learning methods. It utilizes the Armadillo library for linear algebra computations, allowing for relatively low network training times. Training the network on all of the historical data between May 15th, 2013 and January 31st, 2019 takes approximately 10.5 seconds on a single commodity machine (Apple MacBook Pro model A1502, 2.4GHz Intel Core i5-4258U, 8GB RAM). The prediction of a data point as well as an associated trading decision can be obtained in under a millisecond.

\section{Experimental Evaluation} \label{sec:experiments}
Considering the fact that the training data set is time-dependent, Time Series $k$-fold Cross Validation was applied to assess the accuracy of the AI's predictive capabilities and express them numerically in a meaningful manner. The training data set is divided into a training part (80\%) of size $p$ samples and a validation part (20\%) of size $q$ samples, where each sample represents the values of one time unit.

\begin{figure}[t!]
    \centering
    \includegraphics[width=0.41\textwidth]{images/plot-folds.png}
    \caption{Performance of the network under k-fold cross-validation.}
    \label{fig:k-fold}
\end{figure}

The value $k$ represents how many samples are used in the validation step, or how many time units are predicted in each training-validation cycle. Next, a software environment was developed to facilitate this validation method. The software continuously trained the model on samples $1,\dots,p$, produced a prediction for the next hour $p + 1$, validated the prediction using validation sample $p + 1$, trained the model anew on samples $1,\dots,p + 1$, and so on. Setting the value of $k$ to $1$ means $q$ cycles are performed to complete one series of cross validation. An increasing value of $k$ means a decreasing amount of training-validation cycles. After each cycle, the error is calculated and stored. Finally, an error function is used to measure the performance. A range of values for $k$ was tested in order to obtain the lowest margin of error, as shown in Figure \ref{fig:k-fold} where the Root Mean Squared Error is expressed for the number of chosen folds.

As expected, the AI makes fewer predictive errors for a lower value of $k$, as a larger number of folds reduces the amount of data that is used to train the network before the prediction cycle(s). The figure also shows that the prediction error on all four elements stays well below fifty dollars for lower values of $k$. This effectively allows the system to accurately predict up- and downward trends, as the absolute change between each interval is on average higher than the prediction error. Since the initial plan for the trading engine was to make hourly predictions in the first place, the value of $k$ was fixed to one. A new validation cycle is automatically performed each hour after the JavaScript API pulls new data.

\subsection{Parameter Optimization}
The performance of the neural network is highly dependent on several parameters that fine-tune different aspects of the network. These parameters relate to the topology of the network itself, as well as mathematical functions used to increase the speed of learning and overcome local optima.

To optimize the hyperparameters, the validation technique of the previous section was applied while conducting an exhaustive search for the minimal RMSE by modifying a single parameter at a time. While related works using neural networks often remains silent on how and why specific hyper-parameters have been chosen, in practice they are crucial to the accuracy of the RNN and therefore to the success of any RNN-based agent.

\subsubsection{Edge Weight Initialization}
\begin{figure}[t]
    \centering
    \includegraphics[width=0.41\textwidth]{images/plot-edge-weight.png}
    \caption{Impact of the different weight initialization values on the RMSE.}
\label{fig:edge-weight}
\end{figure}

The initialization values of these edges will notably influence the training time of the neural network. A more optimal setting of initial weights will lead to a faster and better convergence of the final weight values, and this will reduce the training time significantly \cite{yam2001feedforward}. As Figure~\ref{fig:edge-weight} shows, random numbers in the range of -0.75 and 0.75 proved to have the lowest RMSE across the different folds of the Bitcoin data set.

\begin{figure}[t]
   \subfloat{\label{fig:optimization}
      \includegraphics[width=0.40\textwidth]{images/plot-optimizers.png}}
    \newline
   \subfloat{\label{fig:rho}
      \includegraphics[width=0.40\textwidth]{images/plot-rho.png}}
    \caption{Impact of parameter optimization on the RMSE.}\label{bs1}
\end{figure}

\subsubsection{Optimization Algorithm}

The role of the optimization algorithm is to adjust the weights on each training iteration of the network, impacting both the speed of convergence as well as the accuracy of the prediction \cite{YAM2000219}. Stochastic Gradient Descent (SGD) methods are the most commonly used optimization algorithms. After testing it with different step sizes against the Adam optimizer \cite{kingma2014adam}, the latter proved to give the best balance of speed and accuracy at a learning rate of 0.01. The results are shown in Figure \ref{fig:optimization}.

\subsubsection{Rho Value (Look-back)}

The rho value sets the length of the back-propagation through time (BPTT) used to train the network. Effectively, it determines how many recent values are used to predict the output at the next step in the time series. As shown in Figure \ref{fig:rho}, the optimal value was found to be around 150.

\subsubsection{Hidden Layer Complexity}

The complexity of the hidden layer refers to the number of hidden nodes used in the FastLSTM layer. Increasing the number of nodes brings the potential to learn more complex rules, but at the same time requires more data to train. Furthermore, it increases the risk of over-fitting and thereby losing the ability to generalize, which results in inaccurate results when the network is used with slightly different data sets. Figure~\ref{fig:hidden-layer} shows that the RMSE was the lowest at approximately 36 hidden layer nodes.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.41\textwidth]{images/plot-hidden-layer.png}
    \caption{Performance of the different hidden layer sizes, measured against the RMSE.}
\label{fig:hidden-layer}
\end{figure}

\subsection{Evaluating Trading Strategy Results}

This section discusses the results obtained by validating the trading engine and each of its five previously described trading strategies on historical data; a technique popularly known as \emph{backtesting}. The set of historical data points acts as a simulation of live input and updates the portfolio by executing hypothetical ('dummy') trades. The performance of each trading strategy in minimizing losses to the portfolio is dependent on the quality of the prediction used to determine a course of action, as well as the quality of the strategy itself. The recurrent neural network in these experiments utilizes the optimized parameter configuration as mentioned in the previous section to maximize prediction accuracy.

Each experiment concerns trading the BTC/USDT pair, with BTC as the base currency and USDT as the quote currency. The goal is to maximize the portfolio's total monetary value in the quote currency. The total value is calculated by converting the portfolio to quote currency at each period.

However, merely achieving a higher monetary value does not imply successful trading if the trader is outperformed by the market itself. If the BTC/USDT value was only ever increasing during a certain time period, then even taking no action whatsoever would have turned a profit. Incidentally, this is known as the \emph{Buy-And-Hold} strategy. With this strategy, a designated asset volume is simply bought and held; no further trading actions are taken. This strategy therefore provides an appropriate baseline for comparison, also known as "trading against Alpha". The results of this comparison and the experiments in general are shown in Figure \ref{fig:portfolio-development}.

\begin{figure} [t]
    \centering
    \includegraphics[width=0.41\textwidth]{images/backtesting/plot-portfolio-growth.png}
    \caption{Portfolio growth of 10.000 USDT between January 31th and September 13th, 2019.}
    \label{fig:portfolio-development}
\end{figure}

\noindent The AI was trained on approximately 90\% of the data set and used to predict the remaining 10\% of the data. Each strategy was tested independently, each receiving an initial portfolio containing 10.000 USDT and subsequently used for trades using the stock market data spanning January 31st to September 13th 2019, approximately 5.400 hours. The Alpha baseline was generated by simply converting the 10.000 USDT into BTC at $t = 0$ and taking no further trading actions. Some limitations have been imposed on the trading agent in order to achieve more stable results. Primarily, a maximum of 25\% of the total portfolio value may be traded at a time. Secondarily, no trades will be made if the expected value increase across the traded volume is smaller than the trading fee (set at 0.1\%).

Figure \ref{fig:portfolio-development} shows that during this part of historical Bitcoin data the RSI strategy managed to achieve the highest portfolio value, with the MACD strategy coming in as a very close second and the DEMA strategy performing only slightly worse. However, all three strategies managed to maintain a stable portfolio value, keeping up with and even outperforming the Alpha baseline. The ROC strategy performs almost as badly as the Random Walk, both barely managing to retain a third of the Alpha baseline value.

Although the AI is able to outperform the Alpha baseline with certain strategies, the question remains whether the profit merits the risks involved. One way of answering this question is by calculating the Sharpe Ratio \cite{sharpe1966mutual}, which gives an indication of how well an investor is compensated by the returns of an asset compared to a risk-free asset with respect to the risk. The Sharpe Ratio is calculated by subtracting some risk-free asset return rate from the considered asset return rate, divided by the standard-deviation of the considered asset return rate. Figure~\ref{fig:sharpe-ratios} shows the monthly Sharpe Ratios of each of the six trading strategies measured over the final seven (virtual) months of the backtesting period, ranging from February 13th to September 13th, 2019. The return rate $R$ is calculated as the percentage of excess returns relative to the start of each month. As for determining the risk-free rate of return $R_f$, monthly U.S. Treasury Bills (T-Bills) made for an appropriate risk-free asset, as the portfolio is measured in USDT. The T-Bills during the backtesting period had a risk-free rate of roughly 2\%, as indicated on the website of the United States Treasury Department \cite{usgovtreasury2019}. To put these numbers into perspective, a Sharpe Ratio above is 1 is generally considered a good investment, above 2 is considered excellent \cite{goetzmann2002sharpening}.

The diagram reveals that even though the RSI and MACD strategies attained the highest portfolio values respectively, it is the DEMA strategy that holds the highest Sharpe Ratio. Since the returns ($R$) of the RSI and MACD strategies were higher and the risk-free rate of return ($R_f$) of each month was the same for every strategy, the explanation is that the RSI and MACD strategies displayed a higher volatility (expressed as $\sigma_R$) than the DEMA strategy during this backtesting period. The Random Walk strategy suffered both suffered major losses and a high volatility, thus its Sharpe Ratio fell below negative. Although the ROC strategy suffered similar losses, it managed to keep a positive Sharpe Ratio due to lower volatility levels.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.41\textwidth]{images/backtesting/plot-sharpe-ratio.png}
    \caption{Sharpe Ratios of the six trading strategies between February 13th and September 13th, 2019.}
\label{fig:sharpe-ratios}
\end{figure}

\section{Conclusion and Future Work}

With RCURRENCY, we have built a viable AI-based trading agent that is able to operate within the highly volatile cryptocurrency market and even generate profit under real market conditions. RCURRENCY uses a FastLSTM-backed recurrent neural network to predict future stock prices. Network inputs are a combination of direct data (four values representing asset value) and seven values derived from this data (technical indicators), while the four outputs represent the predictions for the four asset values. Performance is measured using the Root Mean Squared Error, which emphasizes the gravity of large prediction errors.
Results obtained from the TSKCV suggest that the network is capable of predicting new asset values with approximately 0.4\% error across the four output channels. The strategies have varying rates of success in maintaining a steady portfolio value throughout the historical dataset but every single one performs better than the randomized strategy. Despite the fact that RCURRENCY is tested for Bitcoin, its design allows the system to support other digital assets as well.

No strategy was able to predict large and sudden drops in asset value correctly. This, however, is a common problem of existing AI-based solutions and future research could shed light on the question whether utilizing a broader set of input data can improve the ability to predict such events effectively.

Neural networks generally perform better given more training data. One specific input readily available is volume data, which measures what quantities of an asset have been traded in a given period of time. Unfortunately, volume data had to be excluded as it was unknown from which exchanges the volume data originated, a problem that would likely lead to poor predictions. Finding out the origins of the volume data would provide the system with additional training data.

Despite the fact that RCURRENCY is tested for Bitcoin, the design allows the system to support other digital assets as well. The training data has to include the historical data of the alternative digital asset. Subsequently, the hyper-parameters have to be verified or re-optimized for the newly added data. Another addition, is combining the historical data of multiple assets in different slices per asset, supplying the model with the possibility to utilize the correlations between the slices.

Lastly, it is a common trader's tactic to combine multiple trading strategies to maximize profits. Further simulations could be executed in order to test the performance disparity of the trading engine when combining two, three or even more trading strategies.

\footnotesize{
\bibliographystyle{IEEEtranN}


\section{Introduction}\label{sec:intro}

Microcavity-polaritons hold much promise for the next generation of optoelectronic devices \cite{Sanvitto2016}; including low-threshold lasers and light-emitters \cite{Das2011}, optimisation machines \cite{Berloff2017}, optical circuitry and logic gates \cite{Liew2008}, as well as a platform for quantum simulations and quantum computing \cite{Angelakis2017}. They are quasiparticles that arise from the strong light-matter coupling between excitons in quantum wells or wires, and confined photonic modes in optical microcavities \cite{Hopfield1958, Kavokin2017}. They can form polariton condensates \cite{Byrnes2014} --- a coherent state in which a single mode is macroscopically occupied --- in analogy to atomic Bose-Einstein condensates \cite{Anderson1995}. In addition, modern fabrication techniques permit fine control over the lateral confinement and potential landscape of polaritons. This has enabled the exploration of fascinating phenomena in lattices of polaritons \cite{Schneider2016}, such as exciton-polariton topological insulators \cite{Klembt2018}, Dirac cones in polaritonic graphene \cite{Jacqmin2014}, and Josephson junctions \cite{Lagoudakis2010}.

Much of the pioneering work in polaritonics employed inorganic materials and cryogenic temperatures \cite{Byrnes2014}; this is because the temperature must be lower than the binding energy of the Wannier excitons in these systems. Organic materials are an attractive alternative as they allow room temperature condensation \cite{Cohen2010} due to the much larger binding energy of their Frenkel excitons \cite{Bassani2003}; in addition their processing and fabrication can also be much simpler \cite{Keeling2020}. Organic materials feature high-energy vibrational modes, which enables the efficient condensation of polaritons when matching a vibrational resonance \cite{Plumhof2014, Scafirimuto2017}. This opens the interesting perspective for selectively condensing into different eigenmodes of polariton lattices \cite{Scafirimuto2021}. Whilst there are theoretical treatments of vibron-mediated condensation in planar microcavities \cite{Mazza2009, Mazza2013, Strashko2018, Martinez2019, Justyna2014}, a model that treats this phenomenon in lattices of polaritons is lacking.

In this article, we introduce a general model for simulating organic exciton-polariton condensates in lattices. Specifically we consider lattices of zero-dimensional microcavities, each of which hosts a localised condensate mode, that hybridise due to evanescent coupling. A nonresonant laser creates a reservoir of excitons that can relax to the condensate by emitting a vibron; this is a very efficient process so long as the vibronic energy is similar to the energy difference between the reservoir and a condensate mode. This introduces an energy-selective mechanism, in contrast to nonresonant excitation in inorganic systems. This can be exploited to selectively populate one or more eigenstates of the lattice and should allow the observation of novel dynamics. We highlight the key features of this new model with a few examples. After demonstrating the basics with a single zero-dimensional cavity, we treat the simplest coupled system, a dimer, and observe how either the symmetric or anti-symmetric mode can be selectively populated by modulating the detuning. Then, we examine the dimerised chain and show that a homogeneous pumping of the chain can result in exciting the topological mid gap states, or alternatively the bulk states, by controlling the detuning. Finally, we show a signature of nonreciprocal transport that emerges when driving just one site in the middle of the chain. Taken together, these results exemplify some of the new interesting physics to be found in lattices of organic polaritons that feature different dynamics to their inorganic counterparts.

The rest of this article is organised as follows. In section \ref{sec:analytical_model} we introduce the analytical model; we describe the Hamiltonian and derive the mean-field equations. Then in section \ref{sec:results} we demonstrate the interesting features of these equations for simulating lattices of organic polaritons, by exploring a few examples; a single cavity, a dimer, and a dimerised chain where we show the formation of edge states and a signature of nonreciprocal transport. We discuss the results in section \ref{sec:discussion}, and explain more details of our methods in section \ref{sec:methods}.

\section{Analytical model}\label{sec:analytical_model}

Zero-dimensional microcavities can host localised polariton condensates. For example, one can carve a defect microcavity into a planar Distributed-Bragg-Reflector by Focused-Ion Beam milling; this leads to a discretisation of the lower polariton branch into localised modes within the defect microcavity, with the ground state well separated from higher energy states and featuring a flat dispersion \cite{Scafirimuto2017}. For such systems it is reasonable for us to effectively describe the condensate of the microcavity as a single bosonic mode. This allows us to treat the dynamics of coupled arrays of such cavities with a lattice model. Thus we introduce a Hamiltonian for a general array as

\begin{equation}\label{eqn:Hamiltonian}
    H = \sum_m \left( \omega_\s{C} \psi_m^\dag \psi_m + \omega_\s{R} a_m^\dag a_m + \omega_\s{V} b_m^\dag b_m \right) + H_\s{I} + H_\s{F},
\end{equation}

\noindent where the energy and corresponding creation operator for the condensate are denoted by $\omega_\s{C}$ and $\psi_m^\dag$, for the excitonic reservoir by $\omega_\s{R}$ and $a_m^\dag$, and for the vibrons by $\omega_\s{V}$ and $b_m$; in the case of the $m$th microcavity. By fabricating two (or more) microcavities with a small overlap, there will be evanescent coupling between the condensates, which can be described by

\begin{equation}\label{eqn:condensate_interaction}
    H_\s{I} = \frac{1}{2} \sum_{mn} \Omega_{mn} \left( \psi_m^\dag \psi_n + \psi_n^\dag \psi_m \right),
\end{equation}

\noindent where $\Omega_{mn}=\Omega_{nm}$ is the coupling energy between the $m$th and $n$th cavity. In organic microcavities, scattering from the excitonic reservoir to the polaritonic states are often facilitated by an electronic transition that is not strongly coupled \cite{Keeling2020, Mazza2009, Plumhof2014, Scafirimuto2017}. Here we consider a vibron-mediated relaxation which can be modelled with a Fr{\"o}lich-type Hamiltonian

\begin{equation}
    H_\s{F} = g \sum_m \left( a_m b_m^\dag \psi_m^\dag + a_m^\dag b_m \psi_m \right),
\end{equation}

\noindent where $g$ is the interaction constant between excitons and vibrons. We describe particle decay and a nonresonant pump via a master equation (see Methods), from which we can derive mean-field equations for the number of reservoir excitons $N_m = \langle a_m^\dag a_m \rangle$, the energy flow $J_{mn} = \langle a_m b_m^\dag \psi_n^\dag \rangle$, the flow of condensate polaritons $\rho_{mn} = \langle \psi_m^\dag \psi_n \rangle$, as well as the number of condensate polaritons $\rho_m = \rho_{mm}$. The equations read

\begin{align}
\hbar \frac{\text{d}\rho_{mn}}{\text{d}t} &= -\gamma_\s{C} \rho_{mn} + \text{i} g \left( J_{mn}^* - J_{nm} \right) + \text{i}\sum_{\ell=1}^M \left( \Omega_{m \ell} \rho_{n \ell}^* - \Omega_{n \ell} \rho_{m \ell} \right), \nonumber \\
\hbar \frac{\text{d}N_m}{\text{d}t} &= -\gamma_R N_m + P_m - \text{i} g \left( J_{mm}^* - J_{mm} \right), \label{eqn:mean_field}\\
\hbar \frac{\text{d}J_{mn}}{\text{d}t} &= -\left( \text{i} \delta +  \frac{\Gamma}{2} \right)J_{mn} + \text{i} g N_m^R \left( \delta_{mn} + \rho_{nm} \right) + \text{i} \sum_{\ell=1}^M \Omega_{n \ell} J_{m \ell} \nonumber,
\end{align}

\noindent where $\gamma_\s{C}$, $\gamma_\s{R}$, $\gamma_\s{V}$, and $\Gamma = \gamma_\s{C} + \gamma_\s{R} + \gamma_\s{V}$ are the inverse decay rates of the condensate, reservoir, vibrons, and their sum respectively; $P_m$ is the pumping rate to the reservoir of the $m$th cavity, and $\delta = \omega_\s{R} - (\omega_\s{C} + \omega_\s{V})$ is the detuning between the reservoir and the sum of the condensate and vibron energies. Thus we have a system of $M(2M+1)$ equations and variables (where $M$ is the number of lattice sites) which we can solve for the dynamics. We now gain some insight and intuition with these equations by exploring some simple models, and also observe a few interesting results.

\section{Results}\label{sec:results}

\subsection{Single Cavity}

In this section only, we assume that the energy flows can be adiabatically eliminated; in which case equation \ref{eqn:mean_field} simplifies to

\begin{align}\label{eqn:single_cavity}
    \hbar \frac{\text{d} N}{\text{d}t} &= -\gamma_\s{R} N + P(t) - G N \left( \rho + 1 \right), \\
    \hbar \frac{\text{d} \rho}{\text{d}t} &= -\gamma_\s{C} \rho + G N \left( \rho + 1 \right) \nonumber,
\end{align}

\noindent where

\begin{equation}\label{eqn:condensation_efficiency}
    G = \frac{g^2 \Gamma}{\delta^2 + (\Gamma/2)^2},
\end{equation}

\noindent is the condensation efficiency. These equations clearly show how the efficiency with which reservoir excitons transfer into the condensate depends on the vibronic transition being in resonance with the energy difference between the reservoir and condensate; i.e. $G$ is a maximum for $\delta = 0$. Figure \ref{fig:single_cavity} shows a close correspondence between the value of $G$ and the number of polaritons, obtained from a steady-state solution of equation \ref{eqn:single_cavity} (with $P(t)=P_0$). We note that the condensation process here is quite distinct to both resonant and non-resonant laser excitation of inorganic systems. We have a condensation mechanism which is incoherent and phase non-preserving (like non-resonant pumping of inorganic materials); but which, due to the vibronic resonance, can be used to selectively condense at a specific energy (like resonant pumping of inorganic materials).

\begin{figure*}[t]
	\centering
	\includegraphics[width=0.8\textwidth]{single_cavity.pdf}
	\caption[]{\textbf{Single Cavity.} The number of condensate polaritons $\rho$, and the parameter $G$ that specifies the efficiency of vibron-mediated condensation. The parameters are given in Methods.}
	\label{fig:single_cavity}
\end{figure*}

\subsection{Dimer}

The simplest example of an interacting array of condensates is a dimer; which can also be thought of as the building block (the unit cell) of more complicated bipartite lattices. The hybridisation between the two condensates leads to a symmetric and anti-symmetric mode with an energy splitting equal to the coupling constant. By solving equation \ref{eqn:mean_field} for the steady-state solutions (where $M=2$ and $\Omega_{12}\equiv \Omega_0$) we obtain solutions for the number of condensate polaritons as a function of detuning; see Figure \ref{fig:double_cavity}. When pumping both sites equally (Figure \ref{fig:double_cavity_sub1}) the population on both sites is, unsurprisingly, the same. Moreover, we see a double resonance corresponding to peaks at the symmetric and anti-symmetric modes. When pumping just the left site (Figure \ref{fig:double_cavity_sub2}) we see a curious effect. For large detunings, far away from resonance, the pumped left (red) site has a larger population than the non-pumped right (black) site. However for detunings in-between the two (symmetric/anti-symmetric) resonances the population flows from the pumped to the non-pumped site and builds up there. We now look at a one-dimensional chain constructed from this building block, the dimerised chain.

\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{double_cavity_1.pdf}
  \caption{Both sites pumped the same}
  \label{fig:double_cavity_sub1}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{double_cavity_2.pdf}
  \caption{Only left ($m=1$) is pumped}
  \label{fig:double_cavity_sub2}
\end{subfigure}
\caption{\textbf{Double cavity}. The number of condensate polaritons in the left cavity ($m=1$, red line) and right cavity ($m=2$, black dashed line) as a function of reservoir-condensate detuning $\delta$. (a) Both sites are pumped the same; $P_1(t)=P_2(t)=P_0$. (b) Only the left cavity is pumped; $P_1(t)=P_0$, $P_2(t)=0$. The parameters are given in Methods.}
\label{fig:double_cavity}
\end{figure}

\subsection{The topological midgap states of a dimerised chain}

We now consider a one-dimensional dimerised chain, with alternating coupling strengths between nearest neighbours. As such the condensate interaction Hamiltonian (equation \ref{eqn:condensate_interaction}) embodies the Su-Schrieffer-Heeger model which features topologically protected midgap states that are exponentially localised at the edges. We encode the geometry by setting the nearest neighbour couplings in equation \ref{eqn:mean_field} to $\Omega_{m,m\pm 1} = \Omega_0 \pm (-1)^m \delta_\Omega $ and all other couplings to zero. We consider a homogeneous pumping profile where all sites are pumped equally.

\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{edge_states_ana_1.pdf}
  \caption{Eigenvalues of $\mathcal{H}$}
  \label{fig:edge_states_ana_sub1}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{edge_states_ana_2.pdf}
  \caption{Eigenvector of the midgap state}
  \label{fig:edge_states_ana_sub2}
\end{subfigure}
\caption{\textbf{Eigenstates of the dimerised chain}. (a) The eigenvalues of the condensate Hamiltonian $\mathcal{H}$, with the fill colour of the circles coded by the corresponding value of $G$ (equation \ref{eqn:condensation_efficiency}) and (b) the eigenvector corresponding to one of the midgap states (Eigenvalue no. 10). The parameters are given in Methods.}
\label{fig:edge_states_ana}
\end{figure}

Before simulating the dynamics, let us gain some insight into the system by finding the eigenstates of the condensate interaction Hamiltonian. We can write equation \ref{eqn:condensate_interaction} in matrix form as $H_\s{I} = \vs{\hat{\psi}}^\dag \mathcal{H} \vs{\hat{\psi}}$, where $\vs{\hat{\psi}} = (\psi_1,\psi_2,\dots,\psi_M)^T$; we show the eigenvalues of $\mathcal{H}$ in Figure \ref{fig:edge_states_ana_sub1}, and the eigenvector of one of the midgap states in Figure \ref{fig:edge_states_ana_sub2}. In addition, the fill colour of the circles that mark the eigenvalues denotes the single cavity condensation efficiency $G$, equation \ref{eqn:condensation_efficiency}. In this way, we can see that we have chosen parameters such that the vibronic transition is in resonance with the midgap states (which also happens to be the energy of a single isolated cavity). As such the condensation process should predominantly excite the midgap states, and we expect to see the formation of edge states similar to the pure eigenvectors; see Figure \ref{fig:edge_states_ana_sub2}.

To simulate the dynamics in the time domain we use the Runga-Kutta method of fourth order (see Methods). We start the simulation with an empty system, $\rho_{mn}=N_m=J_{mn}=0$, and switch on a constant homogeneous pump $P_m(t)=P_0$. As we can see in Figure \ref{fig:edge_states_sim}, the system condenses to achieve a steady state solution that features a strong signature of the midgap edge states. Of course the condensation efficiency $G$ is not a delta-function, but features a broadening due to the total decay rate $\Gamma$; and thus the steady-state solution also features a weaker mix of some of the other bulk states, which acts to wash out the pure edge state. Nonetheless a pronounced edge state is indeed witnessed.

In order to test the robustness of this result, we performed stochastic simulations of the same system with a realistic disorder in both the value of $\omega_\s{C}$ at each site and the coupling strengths $\Omega_{mn}$. In both cases the disorder was sampled from a normal distribution with standard deviations $\sigma_\s{C} = 0.04 \Omega_0$ and $\sigma_\Omega = 0.01 \Omega_0$ respectively; see Methods for details. As can be seen from the error bars in Figure \ref{fig:edge_states_sim_sub2}, the key features of the edge states are not washed out by disorder.

\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{edge_states_sim_1.pdf}
  \caption{Simulation of equation \ref{eqn:mean_field}}
  \label{fig:edge_states_sim_sub1}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{edge_states_sim_2.pdf}
  \caption{Bar chart of final time step}
  \label{fig:edge_states_sim_sub2}
\end{subfigure}
\caption{\textbf{Condensation in the edge states}. (a) A time-domain simulation of equation \ref{eqn:mean_field} for the dimerised chain, and (b) a bar chart of the final time step. The error bars denote the standard deviation of the corresponding stochastic simulations. The parameters are given in Methods.}
\label{fig:edge_states_sim}
\end{figure}

\subsection{Nonreciprocal transport}

In the previous section we showed how the system selectively excites the edge states under a homogeneous pumping of all cavities. Of course we are free to choose the pumping pattern as we like; and so we can explore what interesting effects might arise in other circumstances. As one example, we investigate what happens when pumping just one cavity in the middle of a chain with $M=40$ sites. First, we set $P_m(t)=P_0 \delta_{m,20}$ and simulate the dynamics of equation \ref{eqn:mean_field}. As shown in Figure \ref{fig:nonreciprocal_transport_sub1}, the polaritons propagate predominantly to one side of the excitation spot. If we now move the pump to a neighbouring chain, such that $P_m(t)=P_0 \delta_{m,21}$, the polaritons now propagate in the opposite direction; see Figure \ref{fig:nonreciprocal_transport_sub2}. As such the system is displaying a type of nonreciprocal transport, where the propagation direction is determined by which sublattice is being pumped.

\begin{figure}
\centering
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{nonreciprocal_transport_2.pdf}
  \caption{Pumping site no. 5}
  \label{fig:nonreciprocal_transport_sub1}
\end{subfigure}%
\begin{subfigure}{.5\textwidth}
  \centering
  \includegraphics[width=.95\linewidth]{nonreciprocal_transport_1.pdf}
  \caption{Pumping site no. 6}
  \label{fig:nonreciprocal_transport_sub2}
\end{subfigure}
\caption{\textbf{Nonreciprocal transport}. The parameters are given in Methods.}
\label{fig:nonreciprocal_transport}
\end{figure}

\section{Discussion}\label{sec:discussion}

We have introduced a model for describing the dynamics of organic polaritons trapped in lattices of zero-dimensional microcavities; which interact due to evanescent coupling, and condense by emitting a vibrational phonon. Our treatment assumes that condensation occurs in one mode per microcavity, but this can be easily generalised to multimode cavities. Such a lattice model cannot capture physics related to the details of specific trapping geometries, such as the spatial density patterns that arise when condensing in different order orbital modes. Nonetheless, it is an important step in the theoretical treatment of condensation processes in lattice of organic polaritons.

Our results demonstrate the need for methods in modelling lattices of organic polaritons, other than those typically used in inorganic polaritons, such as generalised complex Gross-Pitaevskii equations which do not capture the new physics that arises due to the vibron-mediated condensation. As such, our model contributes to the theoretical toolbox, and should be readily applied to other organic lattices.

\section{Methods}\label{sec:methods}

The mean-field equations \ref{eqn:mean_field} are obtained from the Hamiltonian in equation \ref{eqn:Hamiltonian} and the master equation

\begin{align}
    \hbar \frac{\text{d}\hat{\vs{\rho}}}{\text{d}t} = \text{i} \left[ \hat{\vs{\rho}},H \right] &+ \sum_n\left[\gamma_0 \mathcal{L}(\psi_n) + \gamma_R \mathcal{L}(a_n)  + P_n(t)\left( \mathcal{L}(a_n) + \mathcal{L}(a_n^\dagger) \right) \right] \\  &+ \sum_n \gamma_b\left[ (1+n_\text{th}) \mathcal{L}(b_n) + n_\text{th} \mathcal{L}(b_n^\dagger) \right],\nonumber
\end{align}

\noindent where $\hat{\vs{\rho}}$ is the density matrix, $n_\s{th}$ is the thermal distribution of vibrons, and $\mathcal{L}(a)=(1/2)(2a\hat{\vs{\rho}} a^\dagger - a^\dagger a \hat{\vs{\rho}} - \hat{\vs{\rho}} a^\dagger a)$. We use the identity $\frac{\s{d}\langle {\hat{A}} \rangle}{\s{d}t} = \s{Tr}(\frac{\s{d}\hat{\vs{\rho}}}{\s{d}t} \hat{A})$ which leads to the mean-field equations. We make two approximations; we assume that at any time the number, and thus flow, of vibrons is negligible $\langle b_m^\dag b_n \rangle \ll 1$, and we discard nonresonant processes. With these assumptions we obtain the mean-field equations \ref{eqn:mean_field}.

\begin{table}[h]
    \centering
    \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|c|}
        \hline
        & $\omega_\s{C}$ & $\omega_\s{R}$ & $\omega_\s{V}$ & $\gamma_\s{C}$ & $\gamma_\s{R}$ & $\gamma_\s{V}$ & $\Omega_0$ & $\delta_\Omega$ & $P_0$ & $g$ & $\sigma_\s{C}$ & $\sigma_\Omega$ \\
        \hline
        Figure \ref{fig:single_cavity} & 2526 & $\omega_\s{R}$ & 194 & 4 & 0.7 & 2.9 & - & - & 30 & 0.5 & - & - \\
        \hline
        Figure \ref{fig:double_cavity} & 2526 & $\omega_\s{R}$ & 194 & 4 & 0.7 & 2.9 & 7 & - & 30 & 0.5 & - & - \\
        \hline
        Figure \ref{fig:edge_states_ana} & 2526 & 2720 & 194 & 4 & 0.7 & 2.9 & 33 & 7 & - & - & - & - \\
        \hline
        Figure \ref{fig:edge_states_sim} & 2526 & 2720 & 194 & 4 & 0.7 & 2.9 & 33 & 7 & 30 & 0.5 & 1.32 & 3.3 \\
        \hline
        Figure \ref{fig:nonreciprocal_transport} & 2526 & 2720 & 194 & 4 & 0.7 & 2.9 & 33 & 7 & 30 & 0.5 & - & - \\
        \hline
    \end{tabular}
    \caption{\textbf{Parameters used}. In this table we list the values of the parameters used in all of the Figures in meV.}
    \label{tab:parameters}
\end{table}

To obtain the steady-state solutions we set the time derivatives to zero in equation \ref{eqn:single_cavity} for the single cavity, and in equation \ref{eqn:mean_field} with $M=2$ for the dimer. We can then solve the resulting equations analytically in the case of a single cavity, and numerically with Mathematica \cite{Mathematica} for the dimer. To simulate the dynamics of equation \ref{eqn:mean_field} we use XMDS2 \cite{XMDS2} which is an open source package for numerically integrating initial value problems, including differential equations. We use the solver which employs the Runga-Kutta method of fourth order. To perform the stochastic simulations of Figure \ref{fig:edge_states_sim_sub2} we modify equation \ref{eqn:mean_field} with the substitutions $\omega_\s{C} \xrightarrow[]{} \omega_\s{C} + \mathcal{R}_m(\sigma_\s{C})$ and $\Omega_{mn} \xrightarrow[]{} \Omega_{mn} + \mathcal{R}_{mn}(\sigma_\Omega)$ where $\mathcal{R}_i(\sigma)$ is a random number sampled from a normal distribution with variance $\sigma^2$. We then perform 100 simulations of different instances of the disorder; the bar chart in Figure \ref{fig:edge_states_sim_sub2} show the average, and the error bars show the standard deviation of simulation results. The source code is available for download; see the Supplementary Material. In Table \ref{tab:parameters} we gather the parameters used in all figures in the manuscript.

\FloatBarrier

\printbibliography

\ack

T.S. and M.S. were supported by the Foundation for Polish Science ``First Team'' project No.\ POIR.04.04.00-00-220E/16-00 (originally: FIRST\ TEAM/2016-2/17). M.S. was also supported by the National Science Centre ‘Sonata Bis’ project No. 2019/34/E/ST2/00273. We are very grateful to Rainer F. Mahrt, Thilo St{\"o}ferle, and Darius Urbonas from IBM Research Zurich for many helpful discussions throughout this project. We would also like to thank Anton Zasedatelev (Skoltech) and Evgeny S. Andrianov (MIPT) for insightful discussions at the beginning of this project.

\clearpage

\setcounter{section}{0}
\setcounter{page}{1}

\title{Supplementary Material: Modelling lattices of organic polaritons}

\author{Thomas J Sturges$^1$ and Magdalena Stobi{\'n}ska$^1$}

\address{$^1$Institute of Theoretical Physics, University of Warsaw, ul. Pasteura 5, 02-093, Warsaw, Poland}
\eads{$^1$\mailto{sturges.tom@gmail.com}}
\vspace{10pt}
\begin{indented}
\item[]June 2020
\end{indented}

\section{More information on the code used for our simulations}

The source code used for the simulations in the main text can be downloaded from GitHub \cite{SturgesCode}. This custom script (\code{org-pol-lat.xmds}) was written using the software package XMDS2 \cite{XMDS2}. As explained in the Methods section, we use the Runga-Kutta method of fourth order to solve equation (4) in the main text. Here we provide some additional guidance on the source code. We encourage anybody who is interested in adapting the code to their own needs to get in touch with us.

Table \ref{tab:translation} provides a translation between the LaTeX symbols used in the main text and the plain text variables used in the XMDS2 language. Let us now run through the main parts of the source code (please also consult the XMDS2 documentation). We use units of $meV$ and $ps$. We define our parameters inside \code{<arguments>}, and inside \code{<globals>} we have the helpful function \code{Kron(m, n)} to represent the Kronecker delta function $\delta_{mn}$. We represent the subscripts $m$ and $n$ as dimensions in our simulation (so-called \code{<transverse\_dimensions>}), in addition to the main temporal dimension (the \code{<propagation\_dimension>}). The number of lattice sites is specified by \code{M}.

\begin{table}[h]
\centering
\begin{tabular}{|c|c|}
\hline
LaTeX & XMDS2 \\
\hline
$\omega_\s{X}$, $\gamma_\s{X}$, $g$, $\Gamma$, $\Omega_0$, $P_0$, $\hbar$ &
\texttt{omegaX, gammaX, g, Gamma, Omega0, P0, hbar} \\
\hline
$\rho_{mn}$, $J_{mn}$, $\Omega_{mn}$ & \texttt{rho(m => m, n => n), J(\dots), Omega(\dots)} \\
\hline
$N_m$, $P_m$, $\delta$ & \texttt{NR(m => m), P(\dots), delta(\dots)} \\
\hline
\end{tabular}
\caption{\label{tab:translation}Translation between the LaTeX symbols of the analytical model and the plain text variables used in the code.}
\end{table}

The variables we are solving for ($\rho_{mn}$, $J_{mn}$ and $N_m$) are specified inside the \code{<vector>} named `\code{variables}', which are indexed as e.g. $J_{xy}\xrightarrow[]{}$\code{J(m => x, n => y)}. We use \code{<noise\_vector>} to provide random numbers that introduce disorder for the the couplings $\Omega$ and the detunings $\delta$, which are specified in the vectors `\code{OmegaVec}' and `\code{deltaVec}'. The strength of the disorder is specified by \code{sigmaOmega} and \code{sigmadelta} respectively.

The differential equations are specified inside \code{<integrate>} where we specify `\code{variables}' as the vector to be integrated, and the other vectors as \code{<dependencies>}. Here, for convenience, we introduce the variables \code{sum1} and \code{sum2} which represent the two summations seen in equations \ref{eqn:ddt}. We also introduce variables such as \code{rhonm} as a shorthand for \code{rho(m => n, n => m)}. Lastly we specify the variables we wish to write to memory in \code{<output>}. We are only interested in the real part of $J_{mm}$ so we introduce a save variable \code{rhoReOut} inside \code{<moments>}, and specify that our \code{<dependencies>} are the vector \code{<variables>}.

\subsection{Recreating the figures in the main text}

In the GitHub repository \cite{SturgesCode} there are also folders for each of the figures, allowing one to easily reproduce them. Figures 1, 2 and 3 consist of Mathematica files with a single cell that needs to be evaluated. The Figure 4 directory consists of sub-directories for Figure4a and Figure 4b. Figure4a contains the XMDS2 code and a Mathematica file for plotting the output. Simply compile the program with \code{xmds2 org-pol-lat.xmds} and then run with \code{./org-pol-lat}. Figure4b contains an additional bash script which allows one to automatically run the simulation several times, saving the output to a new folder for each run. It can be run as \code{./stochastic.sh x} where \code{x} is the number of repeat simulations required. Figure5 contains very similar code to Figure4a, except in this case \code{M = 40} and we specify a single pump site with \code{pump = (m == M / 2) ? P0 : 0;}

\section{Checking the analytical calculations with QuantumAlegra.jl}

To double check our analytical calculations we wrote a custom script with the Julia package QuantumAlgebra.jl \cite{QuantumAlgebra} (the package is maintained by Johannes Feist from the Autonomous University of Madrid). The full script can be downloaded \cite{SturgesCode}, and could be helpful to anybody interested in calculating mean-field equations of motions from Lindblad master equations.

\subsection{Preliminaries}

Consider that we have a known Hamiltonian that is a function of the bosonic operators $a_n$. We have some arbitrary operators $\hat{A}_j = f_j(a_1,a_2,\dots,a_N)$, and we want to calculate the mean-field dynamical equations for $A_j = \langle \hat{A}_j \rangle$. To accomplish this we can use the identity $\langle \hat{X} \rangle = \text{Tr}[\rho \hat{X}]$, which gives us

\begin{equation}
\frac{\text{d}A_j}{\text{d}t} = \text{Tr}\left[ \frac{\text{d}\rho}{\text{d}t} \hat{A}_j \right],
\end{equation}

\noindent as well as the Lindblad master equation

\begin{equation}
\frac{\text{d}\rho}{\text{d}t} = \text{i} \left[ \rho, H \right] + \sum_n \gamma_n \mathcal{L}(L_n,\rho),
\end{equation}

\noindent where $L_n$ are the Lindblad jump operators of the system (in simple cases we encounter just the operators $a_n$ and maybe $a_n^\dagger$), and

\begin{equation}
\mathcal{L}(x,y) = \left( x y x^\dagger - \frac{1}{2} \left\{ x^\dagger x, y \right\} \right).
\end{equation}

\noindent To calculate the dynamical equations we exploit the fact that the trace is invariant under cyclic permutations. After substituting equation (2) into equation (1), this allows us to rewrite the mean-field equations as

\begin{equation}\label{eqn:ddt}
\frac{\text{d}A_j}{\text{d}t} = \left\langle -\text{i} [\hat{A}_j,H] + \sum_n \gamma_n \mathcal{M}(L_n,\hat{A}_j) \right\rangle
\end{equation}

\noindent where

\begin{equation}
\mathcal{M}(x,y) = \left( x^\dagger y x - \frac{1}{2} \left\{ x^\dagger x, y \right\} \right).
\end{equation}

\subsection{The code}

We transcribe the function $\mathcal{M}(x,y)$ to code as

\begin{tcolorbox}[colback=red!5!white,colframe=red!5!white]
\begin{lstlisting}
function M(x,y)
    adjoint(x)*y*x - 0.5*adjoint(x)*x*y - 0.5*y*adjoint(x)*x
end
\end{lstlisting}
\end{tcolorbox}

\noindent and implement the RHS of equation \ref{eqn:ddt} as

\begin{tcolorbox}[colback=red!5!white,colframe=red!5!white]
\begin{lstlisting}
function ddt(H,A,(*@$\gamma$@*),L)
    -im*comm(A,H) + sum( map(((*@$\gamma$@*),L) -> (*@$\gamma$@*) * M(L,A), (*@$\gamma$@*), L))
end
\end{lstlisting}
\end{tcolorbox}

\noindent where \code{H} is the Hamiltonian, \code{A} is the operator for which we want to calculate the mean-field dynamics (e.g. the number of excitations \code{adag(1)a(1)}), and \code{$\gamma$} (resp. \code{L}) is a vector whose components are the $\gamma_n$ (resp. $L_n$) we see in equation \ref{eqn:ddt}.

\subsection{A simple example}

Consider the Hamiltonian

$$H = \omega_1 a_1^\dagger a_1 + \omega_2 a_2^\dagger a_2 + \Omega (a_1^\dagger a_2 + a_2^\dagger a_1).$$

\noindent The system experiences losses with rates $\gamma_1$ and $\gamma_2$, and the 1st site is pumped incoherently with a strength $P$. Thus the master equation can be written as

$$\frac{\text{d}\rho}{\text{d}t} = \text{i} \left[ \rho, H \right] + \gamma_1 \mathcal{L}(a_1,\rho) + \gamma_2 \mathcal{L}(a_2,\rho) + P \left[\mathcal{L}(a_1,\rho) + \mathcal{L}(a_1^\dagger,\rho)\right].$$

\noindent The mean-field equation for the number of excitations on site 1 can be calculated as

\begin{tcolorbox}[colback=red!5!white,colframe=red!5!white]
\begin{lstlisting}
H = param(:(*@$\omega_1$@*))*adag(1)*a(1) + param(:(*@$\omega_2$@*))*adag(2)*a(2)
H = H + param(:(*@$\Omega$@*))*(adag(1)*a(2) + adag(2)*a(1))
(*@$\gamma$@*) = [param(:(*@$\gamma_1$@*)),param(:(*@$\gamma_2$@*)),param(:P)]
L = [a(1),a(2),a(1)+adag(1)];
ddt(H,adag(1)a(1),(*@$\gamma$@*),L)
\end{lstlisting}
\end{tcolorbox}

\section{Derivation from the exciton-photon basis}

In the main text we consider a Hamiltonian that describes an excitonic reservoir, polaritonic modes, and vibrons. This model could also describe upper and lower polaritons, in the regime where the lower polaritons are strongly photonic. In this case, the Hamiltonian in the main text can be derived from the Hamiltonian in the exciton-photon basis $H = H_0 + H_\s{F} + H_\s{I}$, where

\begin{align}
    H_0 &= \sum_m \left( \omega_\s{cav} c_m^\dag c_m + \omega_\s{exc} d_m^\dag d_m + \omega_\s{V} b_m^\dag b_m \right) + \Omega_\s{R} \sum_{m} \left( c_m^\dag d_m + d_m^\dag c_m \right), \\
    H_\s{F} &= g_0 \sum_m d_m^\dag d_m \left( b_m + b_m^\dag \right), \\
    H_\s{I} &= \frac{1}{2} \sum_{m\neq n} \Omega_{mn}^0 \left( c_m^\dag c_n + c_n^\dag c_m \right).
\end{align}

\noindent Here, $c^\dagger$ ($d^\dagger$) creates a photon (exciton) with frequency $\omega_\text{cav}$ ($\omega_\text{exc}$), $\Omega_\text{R}$ is the Rabi-frequency, $g_0$ is the interaction constant between excitons and vibrons, and $\Omega_{mn}^0$ is the coupling between the $m$th and $n$th cavity. The operators that diagonalise the single-cavity light-matter part of the Hamiltonian are $\psi_m = d_m \cos\phi - c_m \sin\phi$ and $a_m = c_m \cos\phi + d_m \sin\phi$, where $\phi=(1/2)\s{atan2}\left(\omega_\s{cav}-\omega_\s{exc}, 2\Omega_\s{R}\right)$. In this basis we have

\begin{equation}
    H_0 = \sum_m \left( \omega_\s{C} \psi_m^\dag \psi_m + \omega_\s{R} a_m^\dag a_m + \omega_\s{V} b_m^\dag b_m \right),
\end{equation}

\noindent where

\begin{equation}
    \omega_{\s{R},\s{C}} = \frac{\omega_\s{exc} + \omega_\s{cav}}{2} \pm \sqrt{\left( \frac{\omega_\s{exc} - \omega_\s{cav}}{2} \right)^2 + \Omega_\s{R}^2}.
\end{equation}

\noindent The Fr{\"o}hlich Hamiltonian becomes

\begin{equation}
    H_\s{F} = g \sum_m \left( a_m b_m^\dag \psi_m^\dag + a_m^\dag b_m \psi_m \right),
\end{equation}

\noindent where we have neglected nonresonant terms, and $g = g_0 \sin\phi \cos\phi$. In this work, we assume that $\cos^2\theta \ll \sin^2\theta$, in which case the inter-microcavity interaction term becomes

\begin{equation}
    H_\s{I} = \frac{1}{2} \sum_{m \neq n} \Omega_{mn} \left( \psi_m^\dag \psi_n + \psi_n^\dag \psi_m \right),
\end{equation}

\noindent where $\Omega_{mn} =  \Omega_{mn}^0 \sin^2\phi$. With these approximations we recover the Hamiltonian in the main text.

\printbibliography

\end{document}