Michal Valko
email: michal.valko@inria.fr

Branislav Kveton
email: branislav.kveton@technicolor.com

Hamed Valizadegan

Gregory F Cooper

Milos Hauskrecht

Conditional Anomaly Detection with Soft Harmonic Functions

Keywords: conditional anomaly detection, outlier and anomaly detection, graph methods, harmonic solution, backbone graph, random walks, health care informatics

In this paper, we consider the problem of conditional anomaly detection that aims to identify data instances with an unusual response or a class label. We develop a new non-parametric approach for conditional anomaly detection based on the soft harmonic solution, with which we estimate the confidence of the label to detect anomalous mislabeling. We further regularize the solution to avoid the detection of isolated examples and examples on the boundary of the distribution support. We demonstrate the efficacy of the proposed method on several synthetic and UCI ML datasets in detecting unusual labels when compared to several baseline approaches. We also evaluate the performance of our method on a real-world electronic health record dataset where we seek to identify unusual patient-management decisions.

I. INTRODUCTION

Anomaly detection is the task of finding unusual elements in a set of observations. Most existing anomaly detection methods in data analysis are unconditional and look for outliers with respect to all data attributes [START_REF] Markou | Novelty detection: a review, part 1: statistical approaches[END_REF], [START_REF]Novelty detection: a review, part 2: neural network based approaches[END_REF]. Conditional anomaly detection (CAD) [START_REF] Chandola | Anomaly detection: A survey[END_REF], [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF], [START_REF] Valko | Conditional anomaly detection methods for patient-management alert systems[END_REF] is the problem of detecting unusual values for a subset of variables given the values of the remaining variables. In other words, one set of variables defines the context in which the other set is examined for anomalous values.

CAD can be extremely useful for detecting unusual behaviors, outcomes, or unusual attribute pairings in many domains [START_REF] Das | Anomaly pattern detection in categorical datasets[END_REF]. Examples of such problems are the detection of unusual actions or outcomes in medicine [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF], [START_REF] Valko | Conditional anomaly detection methods for patient-management alert systems[END_REF], [START_REF] Hauskrecht | Conditional outlier detection for clinical alerting[END_REF], investments [START_REF] Rubin | An auctioning reputation system based on anomaly detection[END_REF], law [START_REF] Aktolga | Detecting outlier sections in us congressional legislation[END_REF], social networks [START_REF] Heard | Bayesian anomaly detection methods for social networks[END_REF], politics [START_REF] Kolar | Estimating time-varying networks[END_REF] and other fields [START_REF] Das | Anomaly pattern detection in categorical datasets[END_REF]. In all these domains, the outcome strongly depends on the context (patient conditions, economy and market, case circumstances, etc.), hence the outcome is unusual only if it is compared to the examples with the same context.

In this work, we study a special case of CAD that tries to identify the unusual values for just one target variable given the values of the remaining variables (attributes). The target variable is assumed to take on a finite set of values which we also refer to as labels, because of its similarity to the classification problems.

In general, the concept of anomaly in data in the existing literature is somewhat ambiguous and several definitions have been proposed in the past [START_REF] Markou | Novelty detection: a review, part 1: statistical approaches[END_REF], [START_REF]Novelty detection: a review, part 2: neural network based approaches[END_REF]. Typically, an example is considered anomalous when it is not expected from some underlying model. For the practical purposes of this paper, we define the conditional anomaly detection as follows:

Problem statement ( ): Given a set of n past examples (x i , y i ) n i=1 (with possible label noise), identify and rank instances i in recent m examples (x i , y i ) n+m i=n+1 that are unusual. In this statement, we do not assume that the labels {y i } n i=1 are perfect; they may also be subject to the label noise.

In order to assess the anomalousness of an example, we typically output an anomaly score. One way to define the score is to use a probabilistic model M and calculate the anomaly score as the probability of a different label: P (y = y i |x i , M ). However, the probabilistic model M is not known in advance and must be estimated from available data. This may lead to two major complications that are illustrated in Figure 1: First, an instance may be far from the past observed data points. Because of the lack of the support for alternative responses, it is difficult to assess the anomalousness of these instances. We refer to these instances as isolated points. Second, the examples on the boundary of the class distribution support may look anomalous due to their low likelihood. These boundary examples are also known as fringe points [START_REF] Papadimitriou | Cross-outlier detection[END_REF].

One approach to the CAD task is to construct a classification model on the past data (x i , y i ) n i=1 and apply it to (x i , y i ) n+m i=n+1 to check if the assigned labels (y i ) n+m i=n+1 are correct. However, in that way we would disregard the labels (y i ) n+m i=n+1 which are available and can be utilized to improve the performance of the CAD task by leveraging the interaction between the labels of past examples and the new observed examples.

Because the underlying conditional distribution of the data is unknown, a non-parametric approach that looks for the label consistency of the instances on their neighborhood (e.g., k-nearest neighbor or k-NN) can be very useful [START_REF] Papadimitriou | Cross-outlier detection[END_REF]. One problem with relying on models such as k-NN is that they fail to detect clusters of anomalous instances.

In this paper, we develop and present a new nonparametric method to tackle CAD and its challenges. Our method relies on the similarity graph of instances and attempts to assess whether the response to an input variable agrees with responses of data points in its neighborhood using propagation of labeling information through the graph. Our method differs from typical local neighborhood methods in two important aspects. First, it respects the structure of the manifold and accounts for more complex interactions in the data. Second, it solves the problem of isolated and fringe points by decreasing the confidence in predicting an opposite label for such points through regularization.

Similar to other graph-based approaches for propagation of information (e.g., semi-supervised learning), the solution to our approach demands the computation of the inverse of the similarity matrix that can be challenging for large number of instances. We address this problem by proposing a method for building a smaller backbone graph that approximates the original graph.

In summary, the main contributions of this paper are:

• We introduce a label propagation approach on the data similarity graph for conditional anomaly detection to estimate the confidence of the labels. • We propose a specific regularization that avoids unconditional outliers and fringe points (Figure 1). • We present a compact computation of unconstrained regularization to account for the approximate backbone graph with different node weights. • We propose a scaling approach that adjusts for multitask predictions and makes anomaly scores comparable. • We verify the efficacy of the proposed algorithm on some synthetic datasets, several UCI datasets and a challenging real-world dataset of patients' health records. In the following, we first review the related work in Section II-A, and label propagation on the data similarity graph in Section II-C. In Section III-A, we adopt the label propagation on the data similarity graph for CAD problem and propose a regularization that addresses the isolated and fringe points problems. Next, in Section III-B, we show how to create a smaller backbone graph to deal with more than a few thousand examples. We report the results of our approach on the synthetic datasets, UCI ML datasets, and a real-world medical dataset in Section IV.

II. BACKGROUND

A. Related work

While traditional anomaly detection has been studied for a long time [START_REF] Markou | Novelty detection: a review, part 1: statistical approaches[END_REF], [START_REF]Novelty detection: a review, part 2: neural network based approaches[END_REF], the methods for CAD are relatively new and the research in this area is only emerging [START_REF] Chandola | Anomaly detection: A survey[END_REF], [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF], [START_REF] Song | Conditional anomaly detection[END_REF]. Typically, the existing CAD methods adopt and use classification methods, such as generative models [START_REF] Hauskrecht | Evidence-based anomaly detection[END_REF], [START_REF] Song | Conditional anomaly detection[END_REF], [START_REF] Valko | Distance metric learning for conditional anomaly detection[END_REF], or max-margin classifiers [START_REF] Valko | Conditional anomaly detection methods for patient-management alert systems[END_REF], [START_REF] Hauskrecht | Conditional outlier detection for clinical alerting[END_REF]. In Section IV, we compare our method to these approaches.

The work on CAD, when the target variable is restricted to discrete values only, is closely related to mislabeling detection [START_REF] Brodley | Identifying mislabeled training data[END_REF] and cross-outlier detection [START_REF] Papadimitriou | Cross-outlier detection[END_REF]. The objective of mislabeling detection is to 1) make a yes/no decision on whether the examples are mislabeled, and 2) improve the classification accuracy by removing the mislabeled examples.

Brodley et al. [START_REF] Brodley | Identifying mislabeled training data[END_REF] used different classification approaches to remove mislabeled samples including single and ensemble classifiers. Bagging and boosting are applied in [START_REF] Verbaeten | Ensemble methods for noise elimination in classification problems[END_REF]. Jiang et al. [START_REF] Jiang | Editing training data for knn classifiers with neural network ensemble[END_REF] used an ensemble of neural nets to enhance the performance of a k-NN classifier. Sanchez et al. [START_REF] Sanchez | Analysis of new techniques to obtain quality training sets[END_REF] introduced several other k-NN based approaches including depuration, nearest centroid neighborhood (NCN) and iterative k-NCN. Finally, Valizadegan and Tan [START_REF] Valizadegan | Kernel based detection of mislabeled training examples[END_REF] introduced an objective function based on the weighted nearest neighbor approach and solved it with Newton method.

The main difference between mislabeling detection and CAD is that mislabeling detection identifies and removes the mislabeled examples in order to learn better classifiers, while CAD is interested in ranking examples according to the severity of conditional anomalies in data. This is the main reason our evaluations in Section IV measure the rankings of the cases being anomalous and not the improved classification accuracy when we remove them. Nevertheless, we do compare (Sections IV-A and IV-B) to the methods typically used in mislabeling detection.

Papadimitriou and Faloutsos [START_REF] Papadimitriou | Cross-outlier detection[END_REF] define cross-outliers as examples that seem normal when considering the distribution of examples from the assigned class, but are abnormal when considering the examples from the other class. For each sample (x, y), they compute two statistics based on the the similarity of x to its neighborhood from the samples belonging to class y and samples not belonging to class y. An example is considered anomalous if the first statistic is significantly smaller than the second statistic. However, their method is not very robust to fringe points (Figure 1) [START_REF] Papadimitriou | Cross-outlier detection[END_REF]. The conditional anomaly approach developed in this paper addresses this problem.

B. Notation

We use the following notation throughout the paper. Let (x i , y i ) n+m i=1 be the collection of n past and m recent observed examples. Without loss of generality, we limit y to binary class labels, i.e., y ∈ {±1}. Let G be the similarity graph constructed on the nodes {x i } n+m i=1 with weighted edges W . The entries w ij of W encode the pairwise similarities between x i and x j . We denote by L(W ) the unnormalized graph Laplacian defined as L(W ) = D -W where D is a diagonal matrix whose entries are given by d ii = j w ij .

C. Label Propagation

Label propagation on the graph is widely used for semisupervised learning (SSL). The general idea is to assume the consistency of labels among the data which are 1) close to each other and 2) lie on some structure (a manifold or a cluster). The two examples are the consistency method of Zhou et al. [START_REF] Zhou | Learning with local and global consistency[END_REF] and the harmonic solution of Zhu et al. [START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF]. The inference of the labels by the approach of Zhu et al. [START_REF] Zhu | Semi-supervised learning using gaussian fields and harmonic functions[END_REF] can be interpreted as a random walk on G with the transition matrix P = D -1 W . The harmonic solution satisfies the harmonic property i = 1 dii j∼i w ij j . 1 Harmonic solution and consistency method are the instances of a bigger class of the optimization problems called the unconstrained regularization [START_REF] Cortes | Stability of transductive regression algorithms[END_REF]. In the transductive setting, the unconstrained regularization searches for soft (continuous) label assignment such that it maximizes fit to the labeled data and penalizes for not following the manifold structure:

= min ∈R n ( -y) T C( -y) + T K , (1)
where K is a symmetric regularization matrix and C is a symmetric matrix of empirical weights. C is usually diagonal and the diagonal entries often equal to some fixed constant c l for the labeled data and c u for the unlabeled. In a SSL setting, y is a vector of pseudo-targets such that y i is the label of the i-th example when the example is labeled, and y i = 0 otherwise. Many methods can be derived from (1). For example, for the (hard) harmonic solution K = L(W ), c l = ∞, and c u = 0. Consistency method has K equal to the normalized graph Laplacian K = I -D -1/2 W D -1/2 and c u = c l is set to a non-zero constant. The appealing property 1 j ∼ i denotes that j and i are neighbors in G of ( 1) is that its solution can be computed in closed form as follows [START_REF] Cortes | Stability of transductive regression algorithms[END_REF]:

= (C -1 K + I) -1 y (2)
III. METHODOLOGY

In this section we show how to solve the CAD problem ( ) using label propagation on a data similarity graph and how to compute the anomaly score. In particular, we will build on the harmonic solution approach (Section II-C) and adopt it for CAD in the following ways: 1) show how to compute the confidence of mislabeling, 2) add a regularizer to address the problem of isolated and fringe points, 3) use soft constraints to account for a fully labeled setting, and 4) describe a compact computation of the solution from a quantized backbone graph.

A. Conditional Anomaly Detection

The label propagation method described in Section II-C can be applied to CAD by considering all observed data as labeled examples with no unlabeled examples. The setting for matrix C is dependent on the quality of the past observed data. If the labels of the past observed data (or any example from the recent sample) are guaranteed to be correct, we set the corresponding diagonal elements of C to a large value to make their labels fixed. Notice that specific domain techniques can be used to make sure that the collected examples from the past observed data have correct labels. In this paper, we assume that we do not have access to such prior knowledge and therefore, the observed data are also subject to label noise.

We now propose a way to compute the anomaly score from [START_REF]Novelty detection: a review, part 2: neural network based approaches[END_REF]. The output of (1) for the example i can be rewritten as:

i = | i | × sgn( i ) (3)
SSL methods use sgn( i ) in (3) as the predicted label for i. For an unlabeled example, when the value of i is close to ±1, then the labeling information that was propagated to it is more consistent. Typically, that means that the example is close to the labeled examples of the respective class. The key observation, which we exploit in this paper, is that we can interpret | i | as a confidence in the label. Our situation differs from SSL, as all our examples are labeled and we aim to assess the confidence of already labeled example. Therefore, we define the anomaly score as the absolute difference between the actual label y i and the inferred soft label i :

s i = | i -y i |. (4)
We will now address the problems illustrated in Figure 1. Recall that the isolated points are the examples that are (with respect to some metric) far from the majority of the data. Consequently, they are surrounded by few or no nearby points. Therefore, no matter what their label is, we do not want to report them as conditional anomalies. In other words, we want CAD methods to assign them a low anomaly score. Even when the isolated points are far from the majority data, they still can be orders of magnitudes closer to the data points with the opposite label. This can make a label propagation approach falsely confident about that example being a conditional anomaly. In the same way we do not want to assign a high anomaly score to fringe points just because they lie on a distribution boundary. To tackle these problems we set K = L(W )+γ g I, where we diagonally regularize the graph Laplacian. Intuitively, such a regularization lowers the confidence value | | of all examples; however it reduces the confidence score of far outlier points relatively more. To see this, notice (Section IV-C) that the similarity weight metric is an exponentially decreasing function of the Euclidean distance. In other words, such a regularization can be interpreted as a label propagation on the graph with an extra sink. The sink is an extra node in G with label 0 and every other node connected to it with the same small weight γ g . The edge weight of γ g affects the isolated points more than other points because their connections to other nodes are small.

In the fully labeled setting, the hard harmonic solution degenerates to the weighted k-NN. In particular, the hard constraints of the harmonic solution do not allow the labels spread beyond other labeled examples. However, despite the fully labeled case, we still want to take the advantage of the manifold structure. To alleviate this problem we allow labels to spread on the graph by using soft constraints in the unconstrained regularization problem (1). In particular, instead of c l = ∞ we set c l to a finite constant and we also set C = c l I. With such a setting of K and C, we can solve (1) using (2) to get:

= (c l I) -1 (L(W ) + γ g ) + I -1 y (5) = c -1 l L(W ) + 1 + γ g c l I -1 y. (6)
To avoid computation of the inverse, 2 we calculate (6) using the following system of linear equations:

c -1 l L(W ) + 1 + γ g c l I = y (7)
We then plug the output of ( 7) into (4) to get the anomaly score. We will refer to this score as the SoftHAD score. Intuitively, when the confidence is high but sign( i ) = y i , we will consider the label y i of the case (x i , y i ) conditionally anomalous.

2 due to numerical instability

B. Backbone Graph

The computation of the system of linear equations (7) scales with complexity 3 O(n 3 ). This is not feasible for a graph with more than several thousand nodes. To address the problem, we use data quantization [START_REF] Gray | Quantization[END_REF] and sample a set of nodes from the training data to create G. We then substitute the nodes in the graph with a smaller set of k n distinct centroids which results in O(k 3 ) computation.

We improve the approximation of the original graph with the backbone graph, by assigning different weights to the centroids. We do it by computing the multiplicities (i.e., how many nodes each centroid represents). In the following we will describe how to modify [START_REF] Hauskrecht | Conditional outlier detection for clinical alerting[END_REF] to allow for the computation with multiplicities.

Let V be the diagonal matrix of multiplicities with v ii being the number of nodes that centroid x i represents. We will set the multiplicities according to the empirical prior. Let W V be the compact representation of the matrix W on G, where each node x i is replicated v ii times. Let L V and K V be the graph Laplacian and regularized graph Laplacian of W V . Finally, let C V be the C in (1) with the adjustment for the multiplicities. C V accounts for the fact that we care about "fitting" to train data according to the train data multiplicities. Then:

W V = V W V L V = L(W V ) K V = L V + γ g V C V = V 1/2 CV 1/2
The unconstrained regularization (1) now becomes:

V = min ∈R n ( -y) T C V ( -y) + T K V (8)
and subsequently (6) becomes:

V = C V -1 K V + I -1 y = V -1/2 C -1 V -1/2 (L V + γ g V ) + I -1 y = (c l V ) -1 (L V + γ g V ) + I -1 y = 1/c l V -1 L V + c l γ g + I -1 y
With these adjustments the anomaly score that accounts for the multiplicities is equal to | V -y|. 3 The complexity can be further improved to O(n 2.376 u ) with the Coppersmith-Winograd algorithm.

IV. EXPERIMENTS

To evaluate the performance of our SoftHAD method, we compare it to the following baselines:

• 1-class SVM approach in which we cover each class by a separate 1-class SVM [START_REF] Scholkopf | Estimating the support of a highdimensional distribution[END_REF] with RBF kernel and the anomaly score equals to the distance of the example from the learned boundary of its own class. This method is an example of the traditional anomaly detection method adopted for CAD. • Quadratic discriminant analysis (QDA) model [START_REF] Hastie | The Elements of Statistical Learning[END_REF],

where we model each class by a multivariate Gaussian, and the anomaly score is the class posterior of the opposite class. • SVM classification model [START_REF] Vapnik | The nature of statistical learning theory[END_REF] with RBF kernel where we consider an example anomalous if it falls far on the opposite side of the decision boundary. This method is an example of the classification method adopted for CAD and was used by Valko et al. [START_REF] Valko | Conditional anomaly detection methods for patient-management alert systems[END_REF]. • Weighted k-NN approach [START_REF] Hastie | The Elements of Statistical Learning[END_REF] that uses the same weight metric W as SoftHAD, but relies only on the labels in the local neighborhood and does not account for the manifold structure.

A. Synthetic data

The evaluation of a CAD is a very challenging task when the true model is not known. Therefore, we first evaluate and compare the results of different CAD methods on three synthetic datasets (D1, D2, and D3) with known underlying models that let us compute the true anomaly scores.

We show the three datasets we used in our experiments in We then calculate the true anomaly score as

P (y = y i |x i ) = 1 -P (y = y i |x i ),
reflecting how anomalous the label of the example is with respect to the true model. Each of the methods outputs a score which orders the examples according to the belief of the anomalous labeling. For each of the CAD methods, we assess how much this ordering is consistent with the ordering of the true anomaly score. We did it by counting the number of swapped pairs between the true and predicted ordering which is equivalent to the Wilcoxon score or AUC score. 4The less number of swapped pairs, the higher the agreement score.

Table I compares the agreement scores (or equivalently, the AUCs) of the experiment for all methods. The results demonstrate that our method outperforms all other baselines and comes closest to the order induced by the true model. We also evaluated the linear versions of SVM and 1-class SVM, but the results were inferior to the ones with the RBF kernel. Figure 3 shows the top 5 anomalies identified by each method on D3. We see that only the soft harmonic method was able to identify the top conditional anomalies, which correspond to the examples with switched labels in the middle region that carries a lot of counter-support and hence leads to the highest anomaly score.

B. UCI ML Datasets

We also evaluated our method on the three UCI ML datasets [START_REF] Asuncion | UCI machine learning repository[END_REF], for which an ordinal response variable was available to calculate the true anomaly score. In particular, we selected 1) Wine Quality dataset with the response variable quality 2) Housing dataset with the response variable median value of owner-occupied homes and 3) Auto MPG dataset with the response variable miles per gallon. In each of the dataset we scaled the response variable y r to the [-1, +1] interval and set the class label as y := y r ≥ 0. As with the synthetic datasets, we randomly switched the class labels for three percent of examples. The true anomaly score was computed as the absolute difference between the original response variable y r and the (possibly switched) label. Table II compares the agreement scores to the true score for all methods on (2/3, 1/3) train-test split. Again, we see that SoftHAD either performed the best or was close to the best method.

C. Medical data

In this experiment, we evaluated our method on the problem of detecting unusual patient-management actions [START_REF] Hauskrecht | Conditional outlier detection for clinical alerting[END_REF]. We asked a panel of clinical experts to judge the outputs of the CAD methods for the clinical relevance. 1) Data: We used the data extracted from electronic health records (EHRs) of 4,486 patients as described in [START_REF] Hauskrecht | Conditional outlier detection for clinical alerting[END_REF]. The patients were divided into a train set (2646 patients) and a test set (1840 patients). Patient records were segmented in time (every day of a patient's visit at 8:00am) to obtain 51,492 patient-state instances, such that 30,828 were train and 20,664 test instances. The data in EHRs for these instances were then converted into 9,282 features -a vector representation of the patient state. For every patient-state instance we had 749 decision labels (or tasks) which were possible lab-order and medication decisions with true/false values, reflecting whether a particular lab was ordered or a particular medication was given within a 24-hour period.

2) Evaluation: We evaluated our CAD method on 222 patient-instance/action pairs. We selected these 222 cases such that they represented a wide range of low, medium and high anomaly scores according to the baseline SVM method [START_REF] Valko | Conditional anomaly detection methods for patient-management alert systems[END_REF], [START_REF] Hauskrecht | Conditional outlier detection for clinical alerting[END_REF]. Each instance/action pair was evaluated by three different clinical experts determining whether the action is anomalous and whether this anomaly is clinically relevant. To assess the example, we used the majority rule (two out of three experts). We then evaluated the quality of CAD methods using the area under the ROC (AUC) metric. We compared SoftHAD method to the three baselines 1) weighted k-NN on the same graph 2) SVM with RBF kernel 3) 1-class SVM with RBF kernel described at the beginning of this section.

3) Parameters for the graph-based algorithms: To construct G, we computed the similarity weights as:

w ij = exp -||x i -x j || 2 2,ψ /σ 2 ,
where ψ is a weighing of the features and σ is a length scale parameter. The reason for the different feature weights is the high dimensionality of the data. Without any feature scaling, a distance based on 9K features would make any two points almost equidistant and thus meaningless. Therefore, we weighted the features based on their discriminative power according to the univariate Wilcoxon score [START_REF] Hanley | The meaning and use of the area under a receiver operating characteristic (roc) curve[END_REF]. Next, σ is chosen so that the graph is reasonably sparse [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]. We followed [START_REF] Valizadegan | Kernel based detection of mislabeled training examples[END_REF] and chose σ as 10% of the empirical variance of the Euclidean distances. Based on the experiments, our algorithm is not sensitive to the small perturbations of σ; what is important is that the graph does not become disconnected by having all edges of several nodes with weights close to zero. For each label, we sampled an equal number of positive and negative instances to construct a k-NN graph. We set k = 75, c l = 1 and varied γ g and the graph size. 4) Scaling for multi-task anomaly detection: So far, we have described CAD only for a single task (anomaly in a single label). In this dataset, we have 749 binary tasks (or labels) that correspond to 749 different possible orders of lab tests or medications. In our experiments, we compute the CAD score for each task independently. Figure 4 shows the CAD scores for two of them. CAD scores close to 1 indicate that the order should be done, while the scores close to 0 indicate the opposite. The ranges for the anomaly scores can vary among the different labels/tasks, as one can notice in Figure 4. However, we want to output an anomaly score which is comparable among the different tasks/labels so we can set a unified threshold when the system is deployed in practice. To achieve this score comparability, we propose a simple approach, where we take the minimum and the maximum score obtained for the training set and scale all scores for the same task linearly so that the score after the scaling ranges from 0 to 1.

5) Results: In Figure 5, we fixed γ g = 1 and vary the number of examples we sample from the training set to construct the similarity graph, and also compare it to the weighted k-NN. The error bars show the variances over 10 runs. Notice that the both of the methods are not too sensitive to the graph size. This is due to the multiplicity Figure 4. Histogram of anomaly scores for 2 different tasks. The scores for the top and bottom task range from 0.1 to 0.9 and from 0.25 and 0.61, respectively. The arrow in both cases points to the scores of the evaluated examples, both with negative labels. Despite the score is lower for the bottom task, we may believe that it is more anomalous because it is more extreme within the scores for the same task. adjustment for the backbone graph (Section III-B). Since we use the same graph both for SoftHAD and weighted k-NN, we anticipate that we are able to outperform weighted k-NN due to the label propagation over the data manifold and not only within the immediate neighborhood. In Figure 6, we compare SoftHAD to the CAD using SVM with an RBF kernel for different regularization settings. We sample 200 examples to construct a graph (or train an SVM) and vary the γ g regularizer (or cost c for SVM). We outperform the SVM approach over the range of regularizers. The AUC for the 1-class SVM with an RBF was consistently below 55%, so we do not show it in the figure. We also compared the two methods with scaling adjustment for this multi-task problem (Figure 6). The scaling of anomaly scores improved the performance of both methods and makes the methods less sensitive to the regularization settings.

V. CONCLUSION

We presented a non-parametric graph-based algorithm for conditional anomaly detection. Our algorithm goes beyond exploring just the local neighborhood (nearest neighbor approach) and uses label propagation on the data manifold structure to estimate the confidence of labeling. We evaluated our method on synthetic data, where the true model was known to confirm that the anomaly score from our method outperforms the others in ordering examples according to the true anomaly score. We also presented the evaluation of our method on the real-world data of patient health records, where the true model is not know, but when we used the experts in clinical care to evaluate the severity of our alerts.

In future, we plan to work on the structure anomalies where instead of computing the anomaly score independently for each label, we compute it jointly. With such a structured approach we can avoid the necessity of the score scaling.

Figure 1 .

 1 Figure 1. Challenges for CAD -the squares and the plus signs represent the examples from the two different classes: 1) Fringe points are the examples on the boundary of the class distribution support 2) Isolated points are the examples far from the majority but not within a different class.

Figure 2 .

 2 Figure 2. Synthetic Data. Top: A sample of datasets D1, D2, and D3. Bottom: Synthetic datasets after changing the labels of 3% of the examples.

Figure 2 .

 2 All datasets are modeled as the mixtures of multivariate Gaussians and the class densities we used to generate these datasets vary in locations, shapes and mutual overlaps. Dataset D1 is similar to XOR type of data with one of the classes modeled by a single elongated Gaussian. In D2, the classes overlap and the form of D3 is close to the concentric circles but the clusters are non-overlapping. For each dataset, we sampled 500 examples from the class +1 and 500 examples from the class -1 for the training set and the same number of the examples for the testing set. For each experiment we sample the datasets 100 times. After the sampling, we randomly switch the class labels for three percent of examples for both training and testing set.

Figure 3 .

 3 Figure 3. Black dots depict the top five conditional anomalies based on the score for each of the methods on D3. The top five conditional anomalies according to the true model are in the middle (top left).

Figure 5 .

 5 Figure 5. Medical Dataset: Varying graph size. Comparison of 1) SoftHAD and 2) weighted k-NN on the same graph.

Figure 6 .

 6 Figure 6. Medical Dataset: Varying regularizer 1) γg for SoftHAD 2) cost c for SVM with RBF kernel.

The research work presented in this paper was done when Michal Valko was a PhD student at the University of Pittsburgh and was supported by funds from NIH awarded to University of Pittsburgh.

AUC commonly used for classification is a special case with the true score being ±1.

VI. ACKNOWLEDGMENTS This research work was supported by the grants R21LM009102, R01LM010019, and R01GM088224 from the NIH and the grant IIS-0911032 from the NSF. Its content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH or NSF. Moreover, it was also supported by Ministry of Higher Education and Research, Nord-Pas de Calais Regional Council and FEDER through the "contrat de projets état region 2007-2013", and by PASCAL2 European Network of Excellence. The research leading to these results has also received funding from the European Community's Seventh Framework Programme (FP7/2007-2013) under grant agreement n o 231495 (project CompLACS).


M A Ramos

B Kabtoul

M Hassaine

Calorimetric and thermodynamic study of glass-forming monohydroxy alcohols

Keywords: glass transition, specific heat, polymorphism, boson peak, crystallization monohydroxy alcohols, molecular glasses Glass transition, Monohydroxy alcohols, Specific heat, Polymorphism, Molecular glasses, Crystallization kinetics

Calorimetric and thermodynamic study of glass-forming monohydroxy alcohols

Miguel Angel Ramos, Bisher Kabtoul, Merzak

Introduction

The very nature of glasses (and disordered condensed matter, in general) and of glasstransition related phenomena still represents an open challenge in condensed-matter physics [1] and chemical physics sciences, after more than one hundred years of research on these subjects [START_REF]A non-exhaustive list of good review textbooks is[END_REF].

In particular, much interest has been recently paid to molecular glass-forming liquids. One practical reason for such an interest is that these liquids usually solidify either into crystal or glassy states at temperatures between liquid-nitrogen and room temperatures, therefore providing easy access to the different states of the substance for many experimental techniques. They often also include manifestations of polymorphism and even polyamorphism [START_REF] Senker | [END_REF] (i.e. the apparent existence of first-order transitions between two liquid states of a singlecomponent substance or, more generally, between two distinct amorphous states of that substance), that makes them even more interesting to be studied.

In this paper, we briefly review the main calorimetric and thermodynamic results found in the family of monohydroxy alcohols -or monoalcohols-H(CH 2 ) n OH for n = 1, 2, 3, 4, both from our own experiments during the last decade and from the literature. Firstly, we will discuss how the gradual elongation of the molecule with increasing number of ethyl groups affects the vitrification vs crystallization kinetics for primary monoalcohols (that is, when the hydrogen bond is fixed at the end of the molecule). We will also address briefly the reported presence of crystalline polymorphism in ethanol and polyamorphism in butanol. Then, we will show the noticeable changes that occur in the lattice dynamics and hence in the low-temperature properties of both glass and crystal states when chemical isomers are introduced, namely by comparing 1-propanol vs 2-propanol and 1-butanol vs 2-butanol. Finally, available thermodynamic data of all these monoalcohols, at the glass and crystal melting transitions, will be compiled and globally discussed, checking whether proposed correlations in the literature with liquid fragility could shed light on the observed behaviour.

Primary monohydroxy alcohols: elongating the molecule

Methanol

Pure methanol, CH 3 OH, has a small, 'quasi-spherical' molecule and readily crystallizes, being very difficult to obtain methanol as bulk glass for standard conditions and cooling rates in typical calorimetric experiments [4]. Some researchers have studied methanol with a few molar percent of water as to avoid crystallization, but that is certainly another substance. The structure of both the liquid and the two low-temperature crystal phases found at ambient pressure consists of linear chains of alternating hydrogen bonded molecules [5].

Ethanol

For n = 2, we have ethanol, CH 3 CH 2 OH, which is not an excellent glass-former but at least the glass state can be readily obtained by supercooling the liquid at 30 K/min [6]. However, somewhat slower cooling rates above 6 K/min could suffice [7][8][9][10][11], since the critical cooling rate for ethanol seems to depend on several factors such as minor water impurity, mechanical or geometrical details of the experimental cell driving heterogeneous crystallization and even on the previous thermal history of the liquid [12].

Pure ethanol (see figure 1) exhibits a very rich polymorphism presenting four different solid states [6][7][8]: the ordinary amorphous solid or glass state when the liquid is cooled sufficiently fast, a body-centred cubic plastic crystal (PC) with rotational disorder which by quenching becomes an orientationally-disordered crystal (ODC), and a fully-ordered (monoclinic) crystal state with 4 atoms per unit cell [13], which in turn exists in four distinct phases [9], depending on temperature and thermal history. Therefore, a relatively simple and easily available chemical substance as ethanol offers an extremely useful benchmark to study the glass transition phenomenon and the very nature of glasses and amorphous solids, since one can choose and discriminate the roles played by translational, rotational and orientational disorder.

[Insert figure 1 about here] Very interestingly, both the standard glass transition (amorphous solid ↔ supercooled liquid) and that dynamical freezing of the PC into the ODC were found to occur at the same temperature (T g ≈ T g ´ ≈ 97 K) and with comparable discontinuities in the specific heat [6,8]. Furthermore, the ODC state exhibited, qualitatively and almost quantitatively, the same universal glassy properties at low temperatures such as tunnelling states, the plateau in the thermal conductivity or the boson peak as the truly amorphous, structural glass state [8,10,11,14]. For those reasons, such ODCs are sometimes also named 'glassy crystals' or 'orientational glasses'.

Propanol

By putting a third carbon atom (i.e. adding another ethyl group) to the monoalcohol chain, one gets for n = 3 the best glass-forming substance of the series: 1-propanol or n-propanol. This fact is evidenced by a narrower temperature range for the supercooled liquid than in ethanol. The glass transition temperatures are essentially the same, T g ≈ 98 K for 1-propanol, but the melting point is now 11 K lower: T m = 148 K. Indeed complete crystallization of 1propanol requires annealing over several hours at temperatures within a narrow range about 135 K [15,16]. The crystalline state of 1-propanol was also found to present a monoclinic structure but with 6 atoms per unit cell [16].

Butanol

By further elongating the molecule (n = 4), the glass-forming tendency begins to reverse: 1-butanol glass can be easily obtained when cooling the liquid at a moderately slow rate, but it crystallizes when heating the glass or if the supercooled liquid is kept at certain temperatures long enough. Interestingly, a metastable phase (between the glass and crystal ones) was found in 1-butanol [17], very similar to that previously found in triphenyl phosphite (TPP) [18,19]. In both cases, a new solid phase denoted as glacial phase or glacial state [18], obtained by a first-order, exothermic transformation from the supercooled liquid state was observed. Among other interpretations, this controversial 'glacial phase' was firstly attributed in TPP to a defect-ordered phase within the theory of frustration-limited domains [18], and later was considered (both for TPP and 1-butanol) as a manifestation of liquid-liquid transitions between two distinct liquid (with their corresponding amorphous solid) states [19].

A brief review about polyamorphism with special emphasis on the case of TPP can be found in Ref. [START_REF] Senker | [END_REF]. However, Hédoux et al. [20] have conducted several experiments on TPP, and have shown that the atypical transformation kinetics found is a consequence of a time lag of nucleation, inducing a high nucleation rate which frustrates the crystallization and leads to the so-called "glacial state". By using Raman-scattering, the origin of the frustration of the crystallization in 1-butanol has also been explained [21] from the consideration that the rearrangement of the hydrogen-bonded molecules is unavoidably prevented by the very low molecular mobility of the system only a few degrees above T g . The "glaciation" process would be an abortive crystallization driven by the formation of a layered, two-dimensional hydrogen-bond network. This particular transformation kinetics can nonetheless be described [20,21] by a classical nucleation-growth law. In fact, the temporal evolution of the transformation degree in 1-butanol at two different temperatures above T g were found [21] to follow a clear Avrami sigmoidal shape, indicative of a first-order character for the isothermal transformation. The authors estimated degrees of crystallinity at the end of the transformation processes typically around 70-80% [21]. Our recent experiments [22,23] support this view that the "glacial phase" of 1-butanol or n-butanol is just a mixture of many nanocrystallites embedded in a disordered matrix, essentially due to an aborted crystallization originated by a high nucleation rate in a temperature range where the crystal growth is low. The stable crystalline state was found to present triclinic symmetry with 6 atoms per unit cell [23]. From the case of propanol on, i.e. for n ≥ 3, a new possibility arises to play with a different variable to study its influence on the physical and chemical behaviour of their corresponding liquid and solid phases, namely the existence of some chemical isomers. The chemical isomer of 1-propanol is 2-propanol or iso-propanol, that is also a good glass-former, with glass transition and melting points temperatures, T g ≈ 115 K and T m = 185 K, surprisingly very much higher than those of 1-propanol. Moreover, 2-propanol remarkably exhibits a much higher boson peak (seen in figure 2 as a broad maximum in the Debye-reduced specific heat C p /T 3 ) than in 1-propanol. Interestingly, the same softening of the vibrational lattice was observed in the corresponding crystal states [11,15] with obtained Debye coefficients C Deb = 1.10 mJ•mol -1 K -4 for 1-propanol and C Deb = 1.55 mJ•mol -1 K -4 for 2-propanol. Therefore, both crystalline and glassy states of 2-propanol exhibit a significantly higher density of Debye-like acoustic excitations than those of 1-propanol.

[Insert figure 2 about here] On the other hand, 2-butanol is a better glass-former than 1-butanol and we have not been able to obtain its crystalline state. Nevertheless, the specific heat of the glass states of both isomers has been measured [24,25] at low temperatures, and we have found again that the secondary monoalcohol presents a much higher boson peak than the corresponding primary one, as shown in figure 2 (right panel). To change the hydrogen-bonding location from the end to the centre of the molecule clearly soften the low-frequency vibrational spectrum and enhance its density of states, for both Debye-like acoustic phonons and glassy 'soft modes'.

Finally, we want to mention another possible strategy to investigate the microscopic nature of low-frequency excitations in hydrogen-bonded glasses. It is to assess isotopic effects on their physical properties by employing fully-or partially-deuterated samples. Indeed fullydeuterated ethanol has been measured and compared to normal hydrogenated ethanol [8,10,11], but this method has not been fully exploited yet. In those experiments with fullydeuterated ethanol we showed that the lattice-vibrational contributions for all solid states (crystal, glass and ODC) are well accounted for by Debye's theory, since the cubic coefficients of the low-temperature specific-heat C p data nicely scaled with M 3/2 (M being the molecular mass), as expected [10]. On the other hand, the height of the broad maximum in C p /T 3 (i.e. the boson peak) and the linear term of C p (i.e. the density of tunneling states) for both the glass and the ODC were found to roughly scale with the inverse square root of the moment of inertia of the molecule [11], as previously suggested by some models. Furthermore, these similar correlations in ethanol for glass and ODC phases suggested a common "orientational" origin for the low-temperature excitations, at least in ethanol. Therefore, it would be especially interesting to compare the dynamical and thermal behaviour of monoalcohol glasses with partial deuteration, by selectively replacing hydrogen atoms by deuterium for bonding (OH) groups and non-bonding (methyl or ethyl) groups. By observing different isotopic effects, depending on the location of the deuterated atoms, on relevant experimental data, one could shed light on the unknown microscopic nature of those anomalous low-energy excitations, responsible for the universal low-temperature properties of glasses. As we have shown above, by gradually increasing the 'length' of the molecule for primary monoalcohols, one modifies the balance between crystallization kinetics and glass-forming ability. Whereas pure methanol (a small, 'quasi-spherical' molecule) readily crystallizes and we are not able to obtain it in the glass state for standard conditions and cooling rates, ethanol seems to possess a 'critical' aspect ratio for its molecule, with competing interactions and driving forces which are likely at the origin of its interesting and rich polymorphism, depicted in figure 1. By adding another ethyl group to the hydrocarbon backbone and hence further elongating the 'ellipsoidal' molecule, the best glass-forming primary monoalcohol is observed: 1-propanol. Propanol molecules seem to be long enough as to seriously hinder the crystallization arrangement. With the next in the series, the glass-forming tendency begins to reverse: 1-butanol glass can be easily obtained when cooling the liquid at a moderately slow rate, but it readily crystallizes when heating the glass, though often in a two-step process that could be wrongly taken as a polyamorphic transition, as discussed above in section 2.4. We have recently performed calorimetric, thermal conductivity, Brillouin-scattering and x-ray diffraction experiments [22,23] which strongly support the alternative interpretation that the "glacial phase" of 1-butanol is just a mixture of many nanocrystallites embedded in a disordered matrix. Specifically, Brillouin spectra [22] of the "glacial phase" were composed of two different phonon-like peaks, one of them clearly being the same acoustic peak of the stable crystal. Furthermore, x-ray diffraction experiments [23] showed that the "glaciation" transformation process corresponded to an emergence of the crystalline Bragg peaks onto the amorphous diffraction pattern. When the "glacial phase" was heated to obtain the stable triclinic crystal state, the main Bragg peaks were observed to lie in exactly the same positions, only that much sharper and better defined than in the previous metastable phase, whereas the amorphous background disappeared and allowed the observation of the secondary peaks of the crystalline diffraction pattern. The size of these nanocrystallites was estimated to be of several tens of nanometre. To which extent our finding about the nanocrystalline nature of the so-called glacial state of 1-butanol can be applied to other cases such as TPP remains to be investigated, though earlier works [20] gave evidence in the affirmative.

Interestingly, this glass-forming tendency found in primary monohydroxy alcohols is exactly the same previously observed for the first four members of the alkyl-cyclohexanes family [26]. It suggests therefore that when the hydrocarbon molecules are long enough, they begin to find easier to align themselves and crystallize, and this crossover typically occurs when passing from 3 to 4 carbon atoms in the linear chain.

On the other hand, we have compiled in table 1 the most relevant calorimetric and thermodynamic data concerning the glass transition and the crystal melting for all the considered primary and secondary monoalcohols up to n = 4, as well as reported values for the so-called fragility index [27] ascribed to average relaxation times τ obtained from either viscosity or dielectric measurements. The fragility or steepness index m is defined as the slope of a fragility plot log 10 <τ> vs T g /T evaluated at T g and is usually considered as an essential characteristic of the kinetic and relaxational behaviour of undercooled liquids. Many attempts have been made to correlate these kinetic properties of undercooled liquids with their thermodynamic behavior around the glass transition. Among others, an empirical relation recently suggested [27,28] tried to correlate the fragility m with a dimensionless combination of the specific-heat discontinuity measured at the glass transition temperature ∆C p (T g ) and the melting enthalpy ∆H m , through the equation: Alternatively, by making use of the typical rule T g /T m ≈ 2/3, one could approximately write [27] ( )

m g p S T C m ∆ ∆ = 40 (2)
where ∆S m is the entropy of melting, thereby establishing a direct correlation between a kinetic fragility and a thermodynamic one. Based upon the random first-order transition (RFOT) theory, Lubchenko and Wolynes [34] found a relation very similar to equations ( 1) and ( 2), with only slightly smaller numerical prefactors, hence giving theoretical support to that empirical finding. In Ref. [27], Wang et al. found a good correlation between kinetic and thermodynamic fragilities for more than 50 nonpolymeric glass-forming materials, with a very few exceptions. As can be seen in table 1, primary and secondary monohydroxy alcohols exhibit intermediate values of fragility and, more importantly, such a correlation is not observed at all. Furthermore, one would expect a lower fragility index (i.e. a stronger glassforming liquid) for the case of propanol, the best glass former, and the opposite trend for methanol. This is clearly not the case, neither for the (kinetic) measured fragility nor for the (thermodynamic) calculated one. This failure of the general correlation found in other glassforming systems may be traced back to the earlier idea by Angell [START_REF] Angell | Hydrogen-Bonded Liquids[END_REF], that the intermediate fragility for the viscosity found in hydrogen-bonded liquids could be the compensating result of a very strong kinetic behaviour (with high energy barriers between adjacent minima in the energy landscape) and a very fragile thermodynamic behaviour (with a high density of configurational states accounting for the large values of ∆C p (T g ) typically found in alcohols).

The kinetic-thermodynamic correlation fails even more dramatically for the glass-like transition ODC-PC in ethanol, which in other respects is very similar to the canonical glass transition. The enthalpy of melting for plastic crystals is always very low and consequently they are predicted by eq. ( 1) to be very fragile, what is in clear disagreement with the results found from dielectric measurements in ethanol of the relaxation time, giving m = 60 for the structural glass and only m = 40 for the plastic crystal [29]. In those cases, as well as for the fragility indices m meas of 1-butanol [30] and 2-butanol [33] in table 1, the relaxation-time curves have been determined from the dielectric-loss peak positions of the non-Debye, primary (α) relaxation process, that is the one normally attributed to the glass transition and hence equivalent to the relaxation times obtained from viscosity measurements.

Another interesting concept traditionally used to rationalize the thermodynamic magnitudes for the glass transition and crystal melting processes is that of the number of "beads", i.e. more or less spherical and compact units that constitute the molecule and account for their configurational degrees of freedom. Wolynes and co-workers [34,[START_REF] Stevenson | [END_REF] have proposed to count beads from the entropy of melting. In particular, they found an entropy of fusion per particle for Lennard-Jones spheres, s LJ = 1.68k B , and therefore they proposed a general relation

The corresponding estimations of the number of configurational beads per molecule for the monohydroxy alcohols are also shown in table 1. Using eq. ( 3), Stevenson and Wolynes [START_REF] Stevenson | [END_REF] reported calculated values of 1.3 for methanol, 2.3 for ethanol and 2.6 for 1propanol, whereas we have determined values of 1.6, 2.2 and 2.6, respectively, the small differences arising from the different sources of experimental data employed for the entropy of melting. In addition, we obtain 3.6 beads per molecule of 1-butanol. There is therefore a reasonable increase of N beads with increasing chain length (basically, the number of carbon atoms) for primary alcohols, but this trend is broken for secondary alcohols, where the intermolecular connections change significantly. Eq. ( 3) obviously includes several approximations and oversimplifications, and one should not expect to obtain exact integer numbers. In the RFOT theory, the entropy of fusion comes in only as a way of relating the entropy costs of localizing real molecules that are not spherical to their locations in the crystal and the entropy loss for localizing "beads" which are assumed to be spheres.

In contrast, the effective number of beads for the particular case of the plastic crystal (usually ascribed to the rotation of the methyl group) would be only of around 1 bead per 3 molecules of ethanol. This small number is due to the very small enthalpy of fusion typical of plastic crystals.

Conclusions

Putting together our more recent specific-heat measurements and calorimetric studies of butanol isomers with our earlier experiments on ethanol and propanol, as well as with other data from the literature, we have been able to discuss how either the gradual elongation of the monoalcohol molecule with increasing number of ethyl groups or the hydrogen-bonding location have a strong influence on the glass-forming ability, as well as on low-frequency dynamics and low-temperature properties of molecular solids.

The competition between crystallization rate and glass-forming ability seems to depend mainly on the aspect ratio of the molecules, ranging from the easiest crystal former methanol to the best glass former propanol. Intermediate glass formers as ethanol or butanol has been argued to present the added value of a rich polymorphism -the former-and the possibility of comparing different chemical isomersthe latter-. In particular, it was found that the boson peaks of the structural (amorphous) and orientational (ODC) glasses of ethanol were very similar, whereas those of secondary alcohols were unexpectedly larger than those of their corresponding primary alcohols.

By comparing kinetic and thermodynamic data for this model system of monohydroxy alcohols, we have found that proposed phenomenological correlations for the fragility of undercooled liquids is not fulfilled at all for these -and probably others-hydrogen-bonded liquids, which seem to behave kinetically as strong liquids but thermodynamically as fragile ones. For them, the entropy of crystal melting and the specific-heat jump at the glass transition do not directly correlate with the kinetic fragility index of the supercooled liquid. Furthermore, the glassy transition for the plastic crystal (i.e. the dynamical freezing of rotational disorder into static orientational disorder), that behaves so similarly to the canonical glass transition in many other respects, exhibits a completely disparate behaviour. Table 1. Calorimetric and thermodynamic data of simplest monohydroxy alcohols. T g is glass-transition temperature and ∆C p (T g ) its specific-heat discontinuity; T m and ∆H m are temperature and enthalpy of melting, respectively; m meas and m calc are the measured and calculated (using eq. ( 1)) values for the fragility index; N beads is the estimated number of configurational units or "beads" (see text) as defined by eq. ( 3).

substance T g (K)

∆ ∆ ∆ ∆C p (T g ) (J mol

3

  Chemical isomerism in monohydroxy alcohols: the role of hydrogen-bond location

N

  beads = ∆S m /1.68R

Figure 1 .Figure 2 .

 12 Figure1. Schematic Gibbs-energy phase diagram for ethanol. When the supercooled liquid (SCL) below the melting point (T m = 159 K) is quenched fast enough, the glass state is obtained at T g ≈ 97 K. By heating the glass above T g or by cooling the SCL more slowly, a plastic crystal (PC) is obtained, which becomes an orientationally-disordered crystal (ODC) below a glass-like transition at T g ´ ≈ 97 K. The PC melts at around T m ´ ≈ 125 K and immediately transforms into a fully-ordered monoclinic crystal.

-1 K -1 ) T m (K)

	∆ ∆ ∆ ∆H m
	(

kJ mol -1 ) m meas m calc N beads Refs.

	Page 13 of 13
	methanol	103	30	176	3.85	39	45	1.6	4,27
	ethanol (glass)	97	35.3	159	4.93	55-60 39	2.2 6,12,27,29
	1-propanol	98	45	148	5.40	40	46	2.6	15,27
	1-butanol	111	48	184	9.28	59	32	3.6	22,30
	ethanol (ODC)	97	22.8	127.5	0.66	40	188 0.37	6,12,29
	2-propanol	115	46	185	5.37	?	55	2.1	15,31
	2-butanol	118	39	185	5.97	63	43	2.3	25,32,33

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgments

We acknowledge financial support by the Spanish Ministry of Science within both project FIS2009-09811 and program CONSOLIDER Nanociencia Molecular CSD2007-00010.


Colin R Mcinnes
email: colin.mcinnes@strath.ac.ukwww.strath.ac.uk

Displaced Non-Keplerian Orbits Using Impulsive Thrust

Keywords: Non-Keplerian orbits, circular restricted three-body problem, artificial equilibria, artificial satellites

Displaced non-Keplerian orbits using impulsive thrust

Colin R.

Introduction

The use of continuous low thrust propulsion to generate artificial equilibria was apparently first proposed by [START_REF] Dusek | Motion in the vicinity of libration points of a generalized restricted three-body model[END_REF] who pointed out that a spacecraft could be held at an artificial equilibrium point some distance from a natural three-body equilibrium point. Later, others such as [START_REF] Austin | The ubiquitous solar electric propulsion stage[END_REF], [START_REF] Nock | Rendezvous with Saturn's rings. Anneaux des Planetes[END_REF] and [START_REF] Yashko | Analysis of thruster requirements and capabilities for local satellite clusters[END_REF] noted that propulsive thrust can in principle be used to displace a two-body orbit using continuous thrust normal to the orbit plane. [START_REF] Forward | Light levitated geostationary cylindrical orbits using perforated light sails[END_REF], [START_REF] Baig | Artificial three-body equilibria for hybrid low-thrust propulsion[END_REF] and [START_REF] Heiligers | Displaced geostationary orbits using hybrid low-thrust propulsion[END_REF] have considered displaced geostationary orbits using continuous thrust due to solar radiation pressure. Such displaced two-body orbits will be termed non-Keplerian orbits since the orbit plane does not contain the central mass.

Large families of displaced non-Keplerian two-body orbits have been identified for spacecraft utilising both solar sail (McInnes, 1992a(McInnes, , 1992b) ) and generic low thrust propulsion [START_REF] Dankowicz | Some special orbits in the two-body problem with radiation pressure[END_REF][START_REF] Mcinnes | The existence and stability of families of displaced two-body orbits[END_REF]. More recently displaced orbits have been investigated for a range of applications and propulsion technologies [START_REF] Mengali | Non-Keplerian orbits for electric sails[END_REF]Quarta, 2009, Gong et al, 2009). Studies have determined the conditions for such orbits to exist, their stability properties and their controllability [START_REF] Mcinnes | The existence and stability of families of displaced two-body orbits[END_REF][START_REF] Mcinnes | Dynamics, stability, and control of displaced non-Keplerian orbits[END_REF][START_REF] Xu | Nonlinear dynamical analysis for displaced orbits above a planet[END_REF]. Recently, a large catalogue of such orbits was provided by [START_REF] Mckay | Non-Keplerian orbits using low thrust, high I SP propulsion systems[END_REF] for motion around planetary bodies.

Non-Keplerian two-body orbits are generated using low thrust propulsion by finding conditions for artificial equilibria in a rotating frame of reference.

These equilibria then appear as displaced circular orbits when viewed from an inertial frame [START_REF] Mcinnes | The existence and stability of families of displaced two-body orbits[END_REF]. [START_REF] Scheeres | Stability of hovering orbits around small bodies[END_REF] considered such hover orbits in a frame of reference rotating with an asteroid. Later, work by [START_REF] Sawai | Control of hovering spacecraft using altimetry[END_REF] investigated conditions for stability and controllability of such orbits at irregular shaped asteroids.

Conditions for artificial equilibria in the circular restricted three-body problem can be found using solar sails [START_REF] Mcinnes | Solar sail parking in restricted three-body systems[END_REF], generic low thrust propulsion [START_REF] Morimoto | Artificial equilibrium points in the low-thrust restricted three-body problem[END_REF] and hybrid solar sail propulsion [START_REF] Baig | Artificial three-body equilibria for hybrid low-thrust propulsion[END_REF], while periodic orbits can be found about such artificial equilibria were investigated by [START_REF] Baoyin | Solar sail halo orbits at the Sun-Earth artificial L 1 point[END_REF], [START_REF] Waters | Periodic orbits above the ecliptic plane in the solar sail restricted 3-body problem[END_REF], [START_REF] Baig | Artificial halo orbits for low-thrust propulsion spacecraft[END_REF] and [START_REF] Farres | Periodic and quasi-periodic motions of a solar sail close to SL 1 in the Earth-Sun system[END_REF]. Of particular interest is the case of low thrust propulsion where regions of stable artificial equilibria were found by [START_REF] Morimoto | Artificial equilibrium points in the low-thrust restricted three-body problem[END_REF] which do not exist for the solar sail problem.

Further investigation of the stability of such artificial equilibria has recently been provided by [START_REF] Bombardelli | On the stability of artificial equilibrium points in the circular restricted three-body problem[END_REF].

In this paper, new families of displaced non-Keplerian orbits are generated using impulsive, rather than continuous thrust. The use of impulse thrust was considered by [START_REF] Nock | Rendezvous with Saturn's rings. Anneaux des Planetes[END_REF], [START_REF] Yashko | Analysis of thruster requirements and capabilities for local satellite clusters[END_REF], [START_REF] Mcinnes | Dynamics, stability, and control of displaced non-Keplerian orbits[END_REF], [START_REF] Hope | Pulsed thrust method for hover formation flying[END_REF] and [START_REF] Spilker | Saturn ring observer[END_REF] to generate small displacements away from a circular Keplerian orbit at linear order in the two-body problem. In this paper, the more general problem is considered with arbitrarily large displacements. The displaced non-Keplerian orbits are approximated by a sequence of individual Keplerian arcs connected by impulses. The sequence of arcs is therefore continuous in position, with discontinuities in velocity which are removed using the impulses. For frequent impulses the families of orbits approximate the continuous thrust non-Keplerian orbits found in previous studies. The conditions for displaced non-Keplerian orbits using impulse thrust can be found by linearising the dynamics of the two and three-body problems in the vicinity of an arbitrary reference point in a rotating frame of reference. While the motion of the spacecraft relative to the reference point is considered at linear order, the reference point can be selected arbitrarily so that, when viewed from an inertial frame, the resulting orbit has arbitrarily large displacements away from a Keplerian orbit. Finally, some applications are considered for displaced geostationary orbits.

Two-Body Non-Keplerian Orbits

Continuous thrust

Consider the dynamics of a spacecraft at position r in a frame of reference R(x,y,z) rotating with constant angular velocity ω ω ω ω relative to an inertial frame I(X,Y,Z), as shown in Fig. 1 (McInnes, 1997[START_REF] Mcinnes | Dynamics, stability, and control of displaced non-Keplerian orbits[END_REF]. The spacecraft experiences thrust induced acceleration a and a two-body gravitational potential V. The dynamics of the spacecraft S at position ( )

z y x , , = r
can be defined by an equation of motion in the rotating frame R such that

( ) ( ) r a r ω ω r ω r V dt d dt d ∇ - = × × + × + 2 2 2 (1)
The gravitational potential V of the central body and the thrust-induced acceleration a are defined as

n a r       = - = m T V , µ (2)
where m is the mass of the spacecraft, n is the direction of the thrust vector and µ is the two-body gravitational parameter. The condition for equilibrium in the rotating frame of reference R is therefore given by ( ) ( )

r ω ω r a × × + ∇ = V (3)
The required thrust magnitude and direction for can then be obtained as

( ) ( ) ( ) r ω ω r ω r × × + ∇ = V m T , (4a) ( ) ( ) ( ) ( ) ( ) r ω ω r r ω ω r ω r n × × + ∇ × × + ∇ = V V , (4b)
Using the thrust magnitude and direction defined by Eqs. ( 4 Eq. (4a) are shown in Fig. 2. The units are non-dimensionalised with µ=1 so that the unit of acceleration is the gravitational acceleration at unit distance from the origin. The angular velocity of the rotating frame of reference is selected as 1 = o ω so that a circular Keplerian orbit is found at unit distance in the X-Y plane, as can be seen in the X-Z section shown in Fig. 2. Each point on a contour in Fig. 2 corresponds to a unique displaced circular orbit of period 2π. In addition, radially displaced orbits are generated at the intersection of the contours with the plane. These orbits exist for orbit radii both greater or less than the circular Keplerian orbit a unit distance from the origin, again with period 2π. If the spacecraft S is located along the z-axis then it is in static equilibrium with the thrust induced acceleration balancing the local inverse square gravitational acceleration.

Previous studies have investigated families of orbits with different functional forms of the free parameter o ω , their stability properties and their controllability [START_REF] Mcinnes | The existence and stability of families of displaced two-body orbits[END_REF][START_REF] Mcinnes | Dynamics, stability, and control of displaced non-Keplerian orbits[END_REF][START_REF] Xu | Nonlinear dynamical analysis for displaced orbits above a planet[END_REF]. It will now be shown that the families of continuous thrust non-Keplerian orbits can be approximated using impulsive thrust by creating a sequence of connected Keplerian arcs, with impulses used to remove discontinuities in velocity at the junction between neighbouring arcs.

Impulse thrust

The analysis of Section 2.1 provides the required thrust magnitude and direction for generic families of continuous thrust non-Keplerian orbits. In this section, the analysis is extended to investigate the use of periodic impulses to generate a sequence of connected Keplerian arcs, which approximate continuous thrust non-Keplerian orbits. The S spacecraft will be considered at some arbitrary position r in the rotating frame of reference R, which can be decomposed into a position

( ) ζ η ξ , , = ρ
, relative to an arbitrary point P at position

( ) o o o z y x , , ~= r
, as shown in Fig. 1. Viewed from an inertial frame of reference I, the point P will trace a displaced circular path. The goal is now to find conditions for a sequence of Keplerian arcs ρ ρ ρ ρ, each of which begins and terminates at P with impulses connecting the arcs, as shown in Fig. 3. The point P will therefore trace a displaced non-Keplerian orbit, while the trajectory of the spacecraft S will approximate this displaced orbit by following a sequence of Keplerian arcs.

From Eq. ( 1), the dynamics of the spacecraft S at position

ρ r r + = ~ is described by ( ) ( ) ( ) ( ) ( ) ρ r ρ r ω ω ρ r ω ρ r + -∇ = + × × + + × + + 2 ~2 2 V dt d dt d (5)
which can be written as

( ) ( ) ( ) r ω ω ρ r ρ ω ω ρ ω ρ 2 2 2 × × - + -∇ = × × + × + V dt d dt d (6)
If the displacement ρ ρ ρ ρ of the spacecraft relative to the arbitrary point P at r ~is small, then the gradient of the gravitational potential can be expanded to first order as

( ) ( ) . . . . ~~+       ∂ ∂∇ + ∇ = + ∇ = ρ r r ρ r r r V V V (7)
where the Hessian matrix is evaluated at P. The linearised dynamics of the spacecraft S relative to P are now described by

( ) ( ) ( ) r ω ω r ρ r ρ ω ω ρ ω ρ r r 2 2 2 × × - ∇ -       ∂ ∂∇ - = × × + × + = V V dt d dt d (8)
It can be noted that the last two terms of Eq. ( 8) do not vanish since the point P does not lie on a Keplerian orbit. For the classical linear Clohessey-Wiltshire equations (Clohessey and Wiltshire, 1966), the term

( ) ( ) 0 ~= × × + ∇ r ω ω r V on
a circular Keplerian reference orbit so that Eq. ( 8) becomes homogeneous. In addition, the Hessian matrix of the gravitational potential simplifies to diagonal form.

For the displaced non-Keplerian orbit problem considered in this paper, the linear dynamics are non-homogeneous and can be written in compact form as

Q Mρ ρ N ρ + + = dt d dt d 2 2 (9)
where the skew symmetric gyroscopic matrix N is given by

          - = 0 0 0 0 0 2 0 2 0 o o ω ω N (10)
and the matrix M and non-homogeneous part Q are given by

          + + = zz zy zx yz o yy yx xz xy o xx V V V V V V V V V 2 2 ω ω M (11)           = 3 2 1 Q Q Q Q (12)
Due to the symmetry of the 2-body problem the condition y o =0 can be set to simplify the subsequent analysis. However, this constraint is not used in the analysis of the 3-body problem in Section 3 where the cylindrical symmetry of the problem is broken. The components of the Hessian matrix of the gravitational potential are then given by

( ) ( ) 2 3 2 2 2 5 2 - 3 o o o o o xx z x z x x V + + = µ µ (13a) ( ) 2 5 2 o o o o zx xz z x z x V V + = = µ (13b) ( ) 2 3 2 2 - o o yy z x V + = µ (13c) ( ) ( ) 2 3 2 2 2 5 2 2 2 - 3 o o o o o zz z x z x z V + + = µ µ (13d) where 0 = = yx xy V V and 0 = = zy yz V V
. The non-homogeneous terms are then given by ( )

2 2 3 2 2 1 o o o o o x z x x Q ω µ + + - = (14a) 0 2 = Q (14b) ( ) 2 3 2 2 3 o o o z x z Q + - = µ (14c)
The resulting system of equations is similar to that considered by [START_REF] Sawai | Control of hovering spacecraft using altimetry[END_REF] developed for stability and control analysis of a spacecraft hovering at an asteroid.

Defining the system state vector

( ) T ρ ρ X ɺ , =
for the spacecraft S it can be seen that Eq. ( 9) can then be written in first order form as

B AX X + = dt d (15)
where the system matrix A and non-homogeneous part B are given by

                    - + + = 0 0 0 0 0 0 2 0 0 0 2 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 2 2 zz zx o o yy o xz o xx V V V V V ω ω ω ω A (16)                     = 3 1 0 0 0 0 Q Q B (17)
In order to obtain a general solution to Eq. ( 15) it will be necessary to diagonalise the system matrix A. This can be achieved using the eigenvectors v i (i=1-6) associated with the eigenvalues λ i (i=1-6) of the matrix. The eigenvalues can be determined from the characteristic polynomial of the system matrix A from

= + + + K K K η η η (18) where 2 λ η = and 2 1 2 o K ω = (19a) ( ) ( ) ( ) ( ) 2 5 2 2 2 2 2 3 2 2 2 2 2 2 2 4 2 2 3 3 6 o o o o o o o o o o z x z x z x z x K + - + + + + - = µω µ µ ω (19b) ( ) ( ) ( ) ( ) ( ) ( ) 2 3 2 2 4 2 5 2 2 2 4 2 9 2 2 3 3 2 2 2 2 4 2 2 2 2 2 2 3 3 2 2 2 3 o o o o o o o o o o o o o o o o o z x z x z z x z x z x z x K + + + - + - + - + + = µω µω µ ω µ ω µ (19c)
In general, the eigenvalue spectrum obtained from Eq. ( 18) consists of two pairs of conjugate imaginary eigenvalues and a pair of real eigenvalues. The eigenvectors v i corresponding to each eigenvalue λ i can be used to assemble a

matrix [ ] . . . , . . . , 6 1 v v v V i =
from the eigenvectors which can then be used to diagonalise the system matrix A [START_REF] Betounes | Differential equations: theory and applications[END_REF]. The diagonal form of the system matrix A is then defined by

AV V D 1 - = where             = 6 1 . . λ λ D (20)
Defining a new state vector

X V U 1 - =
it can be seen that Eq. ( 15) can be transformed to ( )

B V U AV V U 1 1 - - + = dt d (21)
Finally, defining a matrix E as ( )

( )             = t exp t exp 6 1 . . λ λ E (22)
it can be shown that Eq. ( 21) can be integrated directly and so the general solution to the inhomogeneous linear system defined by Eq. ( 15) is given by ( ) ( ) ( )

0 VX VE B V D I E V X 1 1 1 - - - + - = t (23)
where ( )

0 X
is a set of arbitrary initial conditions.

This general solution can then be more usefully written in terms of an effective state transition matrix ( )

V VE Φ 1 - = t as ( ) ( ) ( ) Γ X Φ X + = 0 t t (24)
where ( )

B V D I E V Γ 1 1 - - - =
. Decomposing the state transition matrix and nonhomogeneous part as

( ) ( ) ( ) ( ) ( )       = t t t t t 22 21 12 11 Φ Φ Φ Φ Φ (25a)       = 2 1 Γ Γ Γ (25b)
the dynamics of the spacecraft S relative to the point P are defined by

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )       +             =       2 1 22 21 12 11 0 0 Γ Γ ρ ρ Φ Φ Φ Φ ρ ρ ɺ ɺ t t t t t t (26)
which provides the general solution to Eq. ( 8). The trajectory of S relative to P is therefore given by

( ) ( ) ( ) ( ) ( ) 1 12 11 0 0 Γ ρ Φ ρ Φ ρ + + = ɺ t t t (27)
In order to generate a displaced non-Keplerian orbit the condition is set that ( ) ( )

0 ρ ρ = τ
, where τ is the duration of the Keplerian arc described by Eq. ( 27).

This condition ensures that a sequence of connected Keplerian arcs are generated which are continuous in position, but discontinuous in velocity. As will be seen later τ can be selected as a simple fraction of the orbit period of point P. Since the point P is arbitrary, the initial conditions of S can be selected such that ( ) 0 0 = ρ .

From Eq. ( 26) the required initial velocity to ensure that ( ) ( )

0 ρ ρ = τ is then given by ( ) ( ) 1 1 12 0 Γ Φ ρ τ - - = ɺ (28)
The velocity of S on return to P after time τ is then

( ) ( ) ( ) 2 22 0 Γ ρ Φ ρ + = ɺ ɺ τ τ (29)
and so the impulse required to connect neighbouring Keplerian arcs is given by

( ) ( ) ( ) [ ] 2 Γ Γ Φ Φ Φ v - - = ∆ - - τ τ τ (30)
as shown schematically in Fig. 3.

In order to illustrate the use of impulses to generate families of displaced non-Keplerian orbits, the family of continuous thrust orbits shown in Fig. 2 will be considered. Again, the angular velocity of the rotating frame of reference is selected as

1 = o ω
, so that a circular Keplerian orbit will found at unit distance in the plane. The reference point P in the rotating frame will be selected as

=

. Here N=10 is selected for illustration. Using the condition ( ) ( )

0 ρ ρ = τ
a sequence of 10 arcs are generated, where each arc returns to the point P in the rotating frame, as can be seen in Fig. 4. The arcs are generated using Eq. ( 27), the linear solution to Eq. ( 8), with the impulses connecting the Keplerian arcs defined by Eq. ( 30).

In order to evaluate the accuracy of the linear approximation, the twopoint boundary value problem defined by ( ) ( )

0 ρ ρ =
τ can be solved numerically using a shooting method for the full non-linear dynamics defined by Eq. ( 6). Fig. 5 shows the components of the displacement

( ) ζ η ξ , , = ρ
of the spacecraft S from the point P. It can be seen that the linear solution provides a good approximation to the full non-linear dynamics of the problem in this case. Cleary for large displacements about a planetary body, the required impulse magnitude will be extremely large. However, Section 4 provides an example of a displaced geostationary orbit where the impulse magnitude is modest.

Reduction to Clohessy-Wiltshire equations

As the number of Keplerian arcs increases it is expected that non-Keplerian orbits using continuous and impulse thrust will match asymptotically. To demonstrate this, the linearised dynamics at an arbitrary point P can be reduced to the

                          - - + - + = 0 0 0 0 0 0 0 2 0 0 0 2 0 0 0 2 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 3 2 3 2 3 o o o o o o o x x x µ ω ω µ ω ω µ A (31) and                     + - = 0 0 0 0 0 2 2 o o x ω µ B (32)
However, if P is then defined on a circular Keplerian orbit of radius

x o then 3 2 o o x µ ω =
and so Eq. ( 8) reduces to

ξ ξ ω η ω ξ a o o = - - 2 3 2 ɺ ɺ ɺ (33a) η ξ ω η a o = + ɺ ɺ ɺ 2 (33b) ζ ζ ω ζ a o = + 2 ɺ ɺ (33c)
where ( )

ζ η ξ a a a , ,
correspond to the components of continuous thrust induced acceleration a in the local frame attached to P. The continuous thrust acceleration is re-introduced to allow a comparison of continuous and impulse thrust non-Keplerian orbits. For continuous thrust, the required acceleration components for an equilibrium solution at some point

( ) o o ζ ξ , 0 , = ρ
in the rotating frame of reference can be obtained from Eqs. (33) as

a ξ = -3ω o 2 ξ o (34a) 0 = η a (34b) o o a ζ ω ζ 2 = (34c)
which can also be obtained from Eq. ( 3) for small displacements from a circular Keplerian orbit in the plane.

If there is no thrust induced acceleration, then the Clohessy-Wiltshire equations Eq. ( 33) have a general solution (Clohessy and Wiltshire, 1966)

of the form ( ) ( ) ( )       + +       + -       = o o o o o o o o o o t t t ξ ω η ω ω η ξ ω ω ξ ξ 4 2 cos 2 3 sin ɺ ɺ ɺ (35a) ( ) ( ) ( ) [ ]       - + + -       +       + = o o o o o o o o o o o o o t t t t ω ξ ξ ξ ω η ω ω ξ ω ξ ω η η ɺ ɺ ɺ ɺ 2 6 3 cos 2 sin 6 4 (35b) ( ) ( ) ( ) t t t o o o o o ω ω ζ ω ζ ζ sin cos ɺ + = (35c)
A displaced non-Keplerian orbit can then be obtained using periodic impulses, as discussed in Section 2.2. To maintain an out-of-plane displacement, repeated vertical impulses are required such that ( )

o ζ τ ζ =
where τ is again the period between impulses [START_REF] Mcinnes | Dynamics, stability, and control of displaced non-Keplerian orbits[END_REF]. Then, using Eq. ( 35c), the required out-ofplane velocity at the start of the Keplerian arc is given by

      = 2 tan τ ω ζ ω ζ o o o o ɺ (36)
Due to the symmetry of the problem it can be seen that ( ) ( )

0 ζ ζ ɺ ɺ - = T
and so the effective out-of-plane acceleration ζ a provided by the repeated impulses is then

( ) τ ζ / 0 2 ɺ such that       = 2 tan 2 τ ω τ ζ ω ζ o o o a (37)
If the time τ between impulses is small, Eq. ( 37) can then be expanded to yield . . . 12 1 2 4 2

+ + = τ ζ ω ζ ω ζ o o o o a (38)
Comparing Eq. ( 38) and Eq. (34c) it can be seen that to first order the required accelerations are equal, with the higher order terms representing the difference in required mean acceleration between the continuous thrust and impulse cases.

( ) ( ) ( ) ( ) 2 sin 8 2 cos 3 2 sin 3 0 2 τ ω τ ω τ ω τ ω τξ ω ξ o o o o o o - = ɺ (39a) ( ) ( ) ( ) ( ) ( ) 2 sin 8 2 cos 3 2 sin 2 2 cos 6 0 τ ω τ ω τ ω τ ω τ ω τ ω ξ ω η o o o o o o o o - - - = ɺ (39b)
However, it is found that ( ) ( )

0 η τ η ɺ ɺ = and ( ) ( ) 0 ξ τ ξ ɺ ɺ - =
so that only repeated radial impulses are required [START_REF] Mcinnes | Dynamics, stability, and control of displaced non-Keplerian orbits[END_REF]. The effective radial acceleration x a provided by the repeated impulses is then

( ) τ ξ / 0 2 ɺ
. Expanding Eq. (39a) for small periods τ it is found that .

. . 4 3 3 2 4 2 + + - = τ ξ ω ξ ω ξ o o o o a (40)
Again, comparing Eq. ( 40) and Eq. (34a) it can be seen that to first order the required accelerations are equal, again with the higher order terms representing the difference in mean acceleration between the continuous thrust and impulse cases. These approximations provide a simple means of assessing the required acceleration and accumulated effective ∆v for small displacements in the vicinity of a Keplerian orbit, along with the difference in cost of the continuous and impulse thrust modes of operation.

Three-Body Non-Keplerian Orbits

Continuous thrust

The analysis provided in Section 2 can be extended to consider artificial equilibria in the circular restricted three-body problem. Again, artificial equilibria can be generated using continuous low thrust or periodic impulses to connect a sequence of 3-body ballistic arcs. The analysis for the 3-body problem is similar to that provided in Section 2. Here, the 3-body gravitational potential is defined as Using this new potential, Eq. ( 4) can then provide the required thrust magnitude and direction for equilibria in the rotating frame of reference, as used by [START_REF] Morimoto | Artificial equilibrium points in the low-thrust restricted three-body problem[END_REF]. Contours of equal thrust induced acceleration are shown in Fig. 7 for a mass ratio of

6 10 04 . 3 ~-× = µ (Earth-Sun system). The condition ( ) ( ) 0 = × × + ∇ r ω ω r V
can be seen at the L 1 and L 2 natural equilibrium points where 0 = a . For a modest thrust induced acceleration, nested surfaces of artificial equilibria are generated in the vicinity of the natural equilibrium points.

Impulse thrust

Again, the procedure detailed in Section 2.2 can be followed to define the requirements for periodic motion near artificial equilibria using impulsive thrust.

Periodic impulses are now used to generate a sequence of connected ballistic three-body arcs, which approximate continuous thrust artificial equilibria. The spacecraft S will again be considered at some position r in the rotating frame of reference, which can be decomposed into a position

( ) ζ η ξ , , = ρ
, relative to an arbitrary point P at position

( ) o o o z y x , , ~= r
, as shown in Fig. 6. If conditions can be found for the spacecraft S to return to P then an approximation can be found to the families of artificial equilibria using continuous thrust. The goal is therefore to find conditions for a sequence of ballistic three-body arcs ρ ρ ρ ρ, each of which begins and terminates at P with impulses connecting the arcs. The spacecraft S will therefore remain in the vicinity of P.

Using Eq. ( 41), the components of the Hessian matrix of the 3-body gravitational potential defined by Eq. ( 11) are found to be

( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) 2 3 2 2 2 2 5 2 2 2 2 2 3 2 2 2 2 5 2 2 2 2 ~1 ~1 3 1 - 1 1 3 o o o o o o o o o o o o o o xx z y x z y x x z y x z y x x V + + + + - + + + + + + - - + + + + - + + + + - + + - = µ µ µ µ µ µ µ µ µ µ (42a) ( ) ( ) ( ) ( )( ) ( ) ( ) 2 5 2 2 2 2 5 2 2 2 ~1 3 1 1 3 o o o o o o o o o o xy z y x y x z y x y x V + + + + + - - + + + + - + + - = µ µ µ µ µ µ (42b) ( ) ( ) ( ) ( )( ) ( ) ( ) 2 5 2 2 2 2 5 2 2 2 ~1 3 1 1 3 o o o o o o o o o o xz z y x z x z y x z x V + + + + + - - + + + + - + + - = µ µ µ µ µ µ (42c) ( ) ( ) ( ) ( )( ) ( ) ( ) 2 5 2 2 2 2 5 2 2 2 ~1 3 1 1 3 o o o o o o o o o o yx z y x y x z y x y x V + + + + + - - + + + + - + + - = µ µ µ µ µ µ (42d) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3 2 2 2 2 5 2 2 2 2 2 3 2 2 2 2 5 2 2 2 2 ~1 ~1 3 1 - 1 3 o o o o o o o o o o o o o o yy z y x z y x y z y x z y x y V + + + + - + + + + + - - + + + + - + + + + - = µ µ µ µ µ µ µ µ (42e) ( ) ( ) ( ) ( ) ( ) 2 5 2 2 2 2 5 2 2 2 ~1 3 1 3 o o o o o o o o o o yz z y x z y z y x z y V + + + + - - + + + + - = µ µ µ µ (42f) ( ) ( ) ( ) ( )( ) ( ) ( ) 2 5 2 2 2 2 5 2 2 2 ~1 3 1 1 3 o o o o o o o o o o zx z y x z x z y x z x V + + + + + - - + + + + - + + - = µ µ µ µ µ (42g) ( ) ( ) ( ) ( ) ( ) 2 5 2 2 2 2 5 2 2 2 ~1 3 1 3 o o o o o o o o o o zy z y x z y z y x z y V + + + + - - + + + + - = µ µ µ µ (42h) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 3 2 2 2 2 5 2 2 2 2 2 3 2 2 2 2 5 2 2 2 2 ~1 ~1 3 1 - 1 3 o o o o o o o o o o o o o o zz z y x z y x z z y x z y x z V + + + + - + + + + + - - ( ) ( ) ( ) ( )( ) ( ) ( ) o o o o o o o o o x z y x x z y x x Q + + + + + + - + + + + + - + + - = 2 5 2 2 2 2 5 2 2 2 1 ~1 1 1 ~µ µ µ µ µ µ (43a) ( ) ( ) ( ) ( ) ( ) o o o o o o o o o y z y x y z y x y Q + + + + + - + + + + + - = ~1 1 ~o o o o o o o o z y x z z y x z Q + + + + - + + + + + - = µ µ µ µ (43c)
Then, the conditions for the spacecraft S to return to the point P are again defined by Eq. ( 30) using the matrix entries defined by Eq. ( 42) and ( 43).

In order to illustrate the use of impulses to generate families of artificial 3body equilibria, the family of continuous thrust equilibria shown in Fig. 7 will be considered. Eight reference points P 1-8 in the rotating frame will be selected at 26), the linear solution to Eq. ( 8) with the potential defined in Eq. ( 41) and the impulses connecting the arcs defined by Eq. (30). For the example provided in Fig. 8 the impulse magnitude is 1082.2 ms -1 for the Earth-Sun system. Clearly this is large, but is more modest for points closer to the natural equilibria.

Lastly, in order to evaluate the accuracy of the linear approximation, the two-point boundary value problem defined by ( ) ( ) 0 ρ ρ = τ can be solved numerically using a shooting method for the full non-linear dynamics. Fig. 9 shows the components of the displacement

( ) ζ η ξ , , = ρ
of the spacecraft S from the point P. It can be seen that the linear solution provides a reasonable approximation to the full non-linear dynamics of the problem. However, the match is poorer than the 2-body problem due to the large trajectory loops generated by using a smaller number of arcs with N=5, rather than N=10. A larger number of arcs will generate a smaller displacement of S from P and so better accuracy for the linear solution.

Applications

Highly non-Keplerian orbits have a range of potential application, as discussed in Section 1. Here, an example is provided for a displaced geostationary orbit.

Geostationary spacecraft typically use a station-keeping box sized at 0.05 deg of orbit arc. This corresponds to a physical size of order 35x35 km. Figure 10 shows a displaced orbit with a displacement distance of 35 km using both continuous thrust and periodic impulses (N=10). Applications of such displaced orbits include stacking geostationary satellite at crowded longitudes in a similar manner to the solar sail displaced orbits considered by [START_REF] Forward | Light levitated geostationary cylindrical orbits using perforated light sails[END_REF] and later [START_REF] Baig | Light-levitated geostationary orbits are feasible[END_REF] and [START_REF] Heiligers | Displaced geostationary orbits using hybrid low-thrust propulsion[END_REF]. The required impulse magnitude for the displaced orbit shown in Fig. 10 can be determined from Eq. ( 30) as 1.654 ms -1 , so that the ∆v per orbit is 16.54 ms -1 . The required continuous thrust induced acceleration can be determined from Eq. (4a) as 1.851x10 -4 ms -2 , so that the effective ∆v per orbit for the continuous thrust case is 15.99 ms -1 . There is therefore a modest penalty of order 0.5 ms -1 for using the impulse mode rather than the continuous thrust mode. In addition, continuous low thrust allows the use of high specific impulse propulsion, such as an electric thruster with a specific impulse of 3000 s, compared to a bi-propellant thruster with a specific impulse of 320 s used for the impulse mode [START_REF] Heiligers | Displaced geostationary orbits using hybrid low-thrust propulsion[END_REF]. For a large 4000 kg geostationary platform, the propellant mass per orbit for the continuous thrust case is 2.15 kg, which rises to 21.02 kg for the impulse mode.

Conclusions

Families of displaced 2-body non-Keplerian orbits and 3-body artificial equilibria have been extended by considering impulse rather than continuous thrust. The displaced non-Keplerian orbits comprise a sequence of individual Keplerian arcs connected by impulses. The sequence of arcs are therefore continuous in position, with discontinuities in velocity which are removed using the impulses. For frequent impulses it has been shown that the new families of orbits approximate continuous thrust non-Keplerian orbits. The conditions for displaced non-Keplerian orbits using impulse thrust were found by linearising the dynamics of the two and three-body problems in the vicinity of an arbitrary reference point in a rotating frame. While the motion of the spacecraft relative to the reference point was considered at linear order, the reference point was be selected arbitrarily so that when viewed from an inertial frame of reference the non-Keplerian orbit has arbitrarily large displacements.

  ), surfaces of artificial equilibria are generated in the rotating frame of reference, parameterised by the angular velocity of the rotating frame ω = constant thrust induced acceleration generated by

  manner. The characteristic polynomial then reduces to a bi-cubic equation of the form

  -Keplerian orbit will therefore have a radius and displacement distance of 0.5 in non-dimensional units, as shown in Fig.4. In order to generate a sequence of Keplerian arcs which connect to form the displaced orbit, the duration of each arc is selected as a simple fraction of the orbit period

Clohessy

  A and non-homogeneous part B from Eqs. (16-17) reduce to

For

  Again using Eqs. (35) it can be seen that the required initial velocity components are given by

  L 2 points. In order to generate a sequence of ballistic 3-body arcs which return to the point P, the duration of each arc is again selected as a simple fraction N of the orbit period of the two primary masses. Here N=5 is selected for illustration to generate large arcs (where N=10 was used in Section 2.2arcs are generated, where each arc returns to each point P in the rotating frame, as can be seen in Fig.8. The arcs are again generated using Eq. (

  Figure 1.

.

  Figure 6.

1 12 1 12

Acknowledgments

The work presented in this paper was funded through European Research Grant 227571 VISIONSPACE.

( )

where the position of the spacecraft S relative to the two primary masses m 1 and m 2 is defined by ( )

( )

, as shown in Fig. 6. The mass ratio of the problem is defined by

where the gravitational constant and the sum of the primary masses are taken to be unity. Unlike the families of 2-body orbits, the angular velocity of the rotating frame of reference is defined uniquely by the orbital angular velocity of the primary masses and so is taken to be unity.

while the non-homogeneous terms are given by