Mathieu Alain Barrat

Mathieu Génois
email: mathieu.genois@gesis.org

Can co-location be used as a proxy for face-to-face contacts?

Keywords: face-to-face contacts, co-presence, digital epidemiology, complex networks

Can co-location be used as a proxy for face-to-face contacts?

Introduction

In the recent years, several methods have been developed to gather quantitative data on human interactions using wearable sensors and complement more traditional methods based on surveys [START_REF] Mossong | Social contacts and mixing patterns relevant to the spread of infectious diseases[END_REF][START_REF] Danon | Social encounter networks: characterizing great britain[END_REF][START_REF] Eames | Six challenges in measuring contact networks for use in modelling[END_REF]. Current data collection methods range from the use of Bluetooth or WiFi signals in mobile phones [START_REF] Eagle | Reality mining: sensing complex social systems[END_REF][START_REF] O'neill | Instrumenting the city: Developing methods for observing and understanding the digital cityscape[END_REF][START_REF] Scherrer | Description and simulation of dynamic mobility networks[END_REF][START_REF] Vu | Joint bluetooth/wifi scanning framework for characterizing and leveraging people movement in university campus[END_REF][START_REF] Zhang | Towards a temporal network analysis of interactive wifi users[END_REF][START_REF] Stopczynski | Measuring large-scale social networks with high resolution[END_REF] to the specific development of dedicated sociometric sensors [START_REF] Olguín | Social sensors for automatic data collection[END_REF][START_REF] Hashemian | Flunet: Automated tracking of contacts during flu season[END_REF][START_REF] Salathé | A high-resolution human contact network for infectious disease transmission[END_REF][START_REF] Cattuto | Dynamics of person-to-person interactions from distributed RFID sensor networks[END_REF][START_REF] Berke | Objective measurement of sociability and activity: Mobile sensing in the community[END_REF][START_REF] Lucet | Electronic sensors for assessing interactions between healthcare workers and patients under airborne precautions[END_REF][START_REF] Hornbeck | Using Sensor Networks to Study the Effect of Peripatetic Healthcare Workers on the Spread of Hospital-Associated Infections[END_REF][START_REF] Lowery-North | Measuring social contacts in the emergency department[END_REF][START_REF] Toth | The role of heterogeneity in contact timing and duration in network models of influenza spread in schools[END_REF][START_REF] Guclu | Social contact networks and mixing among students in K-12 schools in Pittsburgh, PA[END_REF] and enable researchers to record and measure physical proximity events between individuals in various social contexts. Depending on the specific technology considered however, spatial resolution varies and the resulting "contacts" detected can range from co-presence in a room or a part of a building to close face-to-face encounters. The resulting data is often temporally resolved and has been increasingly used in various contexts including the study of human behaviour, the validation of models of human interactions and data-driven models of epidemic spreading [START_REF] Blower | The importance of including dynamic social networks when modeling epidemics of airborne infections: does increasing complexity increase accuracy?[END_REF][START_REF] Barrat | Measuring contact patterns with wearable sensors: methods, data characteristics and applications to data-driven simulations of infectious diseases[END_REF][START_REF] Eames | Six challenges in measuring contact networks for use in modelling[END_REF].

Despite the increasing availability of techniques to measure even high-resolution temporal contact networks however, a number of limitations remain. In particular, measures cannot be carried out for arbitrarily large population sizes. It is thus of crucial interest to infer contacts or build contact proxies from data with lower spatial resolution data or coming from other sources. In this spirit, several studies have considered the issue of inferring social ties from email exchanges [START_REF] De Choudhury | Inferring relevant social networks from interpersonal communication[END_REF], mobile phone data [START_REF] Eagle | Inferring friendship network structure by using mobile phone data[END_REF], or co-location at geographic scale [START_REF] Crandall | Inferring social ties from geographic coincidences[END_REF]. Other works try to infer close proximity in specific settings from individual attributes [START_REF] Scholz | New insights and methods for predicting face-to-face contacts[END_REF] or from a very precise localisation of individuals [START_REF] Lowery-North | Measuring social contacts in the emergency department[END_REF], or, at geographical scale, from the similarity of the WiFi signals received from a large enough number of WiFi routers [START_REF] Sapiezynski | Inferring person-to-person proximity using wifi signals[END_REF].

Here instead, we do not try to infer specific contacts between pairs of individuals but rather investigate if a coarse co-location information on individuals allows us to reach an overall picture of the contact patterns in the population of interest. To this aim, we leverage several data sets collected by the SocioPatterns collaboration [27,[START_REF] Cattuto | Dynamics of person-to-person interactions from distributed RFID sensor networks[END_REF] in various contexts: these data include both detailed information about close, face-to-face encounters between individuals and a location tracking of individuals with low spatial resolution. It is thus possible to build two temporal networks where nodes represent individuals and links correspond respectively to a face-toface contact or to a co-presence event, where co-presence is defined with respect to the localisation of two individuals within the same spatial area. We first compare the structural and statistical properties of these two temporal networks and show that they share some important properties, although the co-presence network is much denser, due to the lower spatial resolution involved in its definition. We thus investigate several methods of down-sampling the co-presence signal in order to create surrogate contact networks, in the spirit of [START_REF] Génois | Compensating for population sampling in simulations of epidemic spread on temporal contact networks[END_REF][START_REF] Sapienza | Estimating the outcome of spreading processes on networks with incomplete information: a mesoscale approach[END_REF], and compare these surrogate data to the actual networks of face-to-face contacts. We focus first on several statistical characteristics of temporal and aggregated networks, and quantify the ability to identify central nodes in the contact network from the surrogate data. We then consider the possibility to use the surrogate data in numerical simulations of data-driven models for epidemic spread. In particular, we compare the outcome of simulations of a standard model of epidemic propagation when surrogate or actual contact data are used, and we explore the possibility to identify the most efficient containment strategies from this limited information [START_REF] Smieszek | A low-cost method to assess the epidemiological importance of individuals in controlling infectious disease outbreaks[END_REF]. Our results turn out to depend strongly on the data collection context, highlighting the limitations of coarse co-presence networks with respect to detailed face-to-face data.

The co-presence network

Data sets

We use data collected by the SocioPatterns collaboration in various contexts. These data were gathered using wearable sensors able to detect face-to-face close range proximity (1.5 m) of participants wearing the sensors on their chests. In addition, the sensors broadcast a signal that can be received by RFID readers located in the environment. In open space, each reader can receive signals from sensors situated within a range of ∼ 30 m, while the actual reception range in a building depends on its specific structure and on the nature of its walls, floors and ceilings. Each reader thus defines a coarse spatial area and the sensors' signals can be followed when the individuals carrying them change area. For each sensor, we define its "spatial location" at each time as the set of readers receiving its broadcasted signal at this time, and we define two individuals to be in co-presence if they share the same spatial location, i.e., the same exact set of readers have received signals from both individuals. We use data sets from various social contexts: a workplace, with data collected in two different years (InVS13, InVS15), a hospital (LH10), a primary school (Lyon-School), a scientific conference (SFHH) and a high school (Thiers13), see Table 1. In each case, we thus consider a temporal network of face-to-face contacts and a temporal network of co-presence between individuals, both at the temporal resolution of 20s. A contact (resp. co-presence) event between two individuals is then defined as a set of successive time-windows of 20s during which the individuals are detected in contact (resp. co-presence), while they are not in the preceding nor in the next 20s time window. While the conference data does not include any other information on the participants and does not exhibit any particular group structure [START_REF] Stehlé | Simulation of an seir infectious disease model on the dynamic contact network of conference attendees[END_REF], the other populations under study can each be divided into groups: departments for the workplace, classes for the school and the high school, and roles (patients, doctors, nurses) in the hospital. In these cases, the overall structure of networks aggregated over a certain time window can be summarised, in addition to usual quantities such as the density, the clustering coefficient or the degree distribution, by contact (resp. co-presence) matrices that give the fraction of pairs of individuals of different groups who have been in contact (resp. in co-presence). Moreover, temporal features of interest include the distributions of durations of contact or co-presence events, of the time elapsed between successive events, of the numbers and aggregated durations of such events between pairs of individuals (the latter quantity yields a natural definition of the weight of a link between individuals in the aggregated network).

We will show in the main text the results corresponding to the InVS15 data set, and we refer to the Supplementary Information for the results obtained with the other data sets. We make also available as Supplementary Files the temporally resolved contact and co-presence networks.

Co-presence and contact networks

We first compare some features of the co-presence and contact networks, both temporal and for networks aggregated either on the whole data gathering period or over daily temporal windows. We show in Fig. 1 the distributions of event and interevent duration, as well as the distributions of number and cumulative duration of events for individual pairs. The co-presence events show broad distributions of these quantities, similarly to the contact events and with similar slopes: using only copresence data yields approximate information on the functional shape of the contact Temporal distributions for co-presence and contact events -InVS15. We show for both the contact and co-presence of the same data set the distributions of event and inter-event duration, link weights (as total contact duration) and number of contacts per link. duration distributions. The distributions of durations and numbers of events are however typically broader for co-presence, with heavier tails, and the distribution of inter-event durations tend to be less broad (see also SI). This is not surprising as the criterion for being in co-presence is less strict than for being in contact. We observe the strongest differences between co-presence and contact distribution functional shapes for the primary school data. This could be related by the fact that the spatial resolution is in that case quite low, with all the schoolyard being covered by one single reader, and some readers covering more than one classroom. Overall, using only co-presence data would lead to over-estimations of the contact durations and aggregate durations. For each data set we compute the cosine similarity between the average daily contact matrix and the co-presence matrix, as well as for the contact matrices obtained for each sampling method of the co-presence data, averaged over 100 realisations for each sampling method.

We compare moreover in Figs. 2 -3 and Tables 2 -3 the overall structures of the contact and co-presence networks, aggregated over daily time windows. The co-presence aggregated networks are much denser than the contact network, with a larger average degree, a larger average clustering coefficient and larger cliques, as expected once again given the lower spatial resolution required for co-presence We compare the average degree ( k) network density (ρ), clique number (ω) and average clustering (c) of daily aggregated networks, for the contact network (c subscript), the co-presence network ( subscript), and the sampled co-presence networks (subscripts 1 to 3 according to the sampling method). Values are averaged over all the days of the study. In the case of SFHH, since on the second day there was activity only during the morning, only the values of the first day are reported. *The network is too large and too dense for the clique number to be determined in reasonable time via the usual algorithm. events. In some cases (school, conference), the aggregated networks are even close to being fully connected (see for illustration Fig. 3). Despite this strong difference in the overall density of links, the contact and co-presence matrices giving the density of links between and within each group, averaged across days, are very similar (Table 2). The similarity is particularly high for the hospital data and, even for the lower value obtained for the high school data, the matrices displayed in the SI show that the overall structure in classes and groups of classes can be inferred from the co-presence data alone. Given the simultaneous discrepancies in density values and similarities in the networks group structures, we investigate if the data exhibits a scaling law between the number of individuals present in an area and their contact activity, as found at geographical scale in phone communication [START_REF]The scaling of human interactions with city size[END_REF] and Twitter data [START_REF] Tizzoni | The scaling of human contacts and epidemic processes in metapopulation networks[END_REF]. Figure 4 and the similar figures shown in SI show the results obtained in the various contexts. Apart from the office cases (InVS13 and InVS15), we observe indeed a correlation between the median of the number of contacts and the number of individuals present. This correlation exhibits a power law shape, with an exponent around 1.5 (see figures in SI). However, huge, context-dependent fluctuations are observed. For instance, in the InVS15 case, the trend is strongly influenced by the numerous instances of an absence of contacts despite potentially large values of the number of individuals present in the area. This is a consequence of the fact that a given reader can receive signals from the sensors of individuals located in different offices. In other areas such as a cafeteria, many more contacts occur with potentially a similar or even smaller number of individuals. Overall, very large fluctuations of the number of contacts, at given number of individuals present, are thus observed, because on the one hand of the low spatial resolution of the co-presence data, and on the other hand of the variety of contexts corresponding to the areas covered by different RFID readers. The stronger correlation is observed for the SFHH conference data, probably because the various areas covered by the readers corresponded to similar contexts, namely different areas of the exhibition and poster rooms.

3 Sampling co-presence data

Sampling methods

As the temporal network of co-presence bears some similarities with the actual contact data, but contains much more events and leads to much denser aggregated networks, we consider the possibility to down-sample the co-presence data: for each pair of individuals, each contact event is indeed included in a co-presence event of the same individuals. Each co-presence event might thus correspond to one or more contact events. As we cannot determine exactly the correct down-sampling to be performed if we have access only to co-presence data, we study here three simple sampling methods. We remind here that we do not try to infer the real contacts but rather to obtain a down-sampled version of the co-presence network that is statistically similar to the real contact data. Moreover, as the total number and duration of actual contacts cannot either be easily guessed from the co-presence data alone, we consider the actual total contact time T c as the (only) parameter of the sampling, and we fix it to its empirical value. The sampling methods we consider are the following:

• Sampling 1: Sampling of co-presence times. We define a co-presence list as a list of individuals present at the same time t in the same area. Each copresence list is thus stamped with its time of occurrence t. We create n copies of each co-presence list , where n is the number of distinct individuals in , and create in this way of a global pool of co-presence lists. We then sample T c lists uniformly at random from the pool without replacement. Each list has thus a probability proportional to the number of individuals it contains to be chosen. From each chosen list, we choose at random a pair i, j of individuals, obtaining a triplet (t, i, j) where t is the time-stamp of the list (we take care to avoid repetitions: if (t, i, j) has already be obtained in a previous random draw, we repeat the random selection). The sampled temporal co-presence network (i.e., the surrogate contact network) is formed by the union of these triplets. • Sampling 2: Sampling of co-presence times with completion. We constitute a pool of lists exactly like in the previous method. We then sample a triplet (t, i, j) as in the previous method, and add all the other triplets (t , i, j) that belong to the same co-presence event to create the surrogate contact event. We iterate this until we reach a cumulative contact time T c , while discarding repetitions.

• Sampling 3: Sampling of co-presence events. We consider directly the list of co-presence events between individuals, (t, i, j, τ ) (co-presence event between individuals i and j, starting at time t and with duration τ ), and sample events from this list, without replacement, adding them to the list of surrogate contact events until we reach a cumulative contact time T c . For each data set, we create 100 instances of surrogate contact networks for each sampling method. We compare in the following the properties of these surrogate contact networks with the real face-to-face contact data. Figure 5 Properties of the sampled co-presence networks -InVS15. We compare several properties of the contact network from the original data set with the surrogate contacts obtained by sampling of the co-presence data: overall timeline of contact activity, distributions of degree, weight w and number of contacts per link n in the network aggregated over the whole data collection period, and distributions of the contact duration τc and inter-contact duration τ i .

Figures 5 -6 and Tables 2 -3 provide elements of comparison between the surrogate contact networks and the empirical data (see also SI). The first observation is that the contact activity timelines are in general broadly recovered, except for the primary school (see SI), while the detailed intra-day activity variations are not always properly reconstructed in the surrogate data (except for the hospital data, see SI). The strongest deviations are observed for the second sampling method for the conference and high school data.

The first sampling method, given it samples separately times of co-presence, yields an exponential distribution of surrogate contact duration, in contrast with actual data and other sampling methods in which broad distributions are observed. Broad distributions of inter-contact durations and of the numbers of contacts between individuals are also obtained, with however slopes that depend on the context. For instance, the second sampling method systematically leads to a distribution of contact durations that is broader than for the real contacts. The third method yields a distribution of contact durations similar to the real one for the InVS13, LH10, and SFHH cases, but gives results similar to the second method in the other cases.

We now turn to the properties of networks aggregated over daily periods or over the whole data collection. At the daily level, we show in Table 2 that the similarity of the contact matrices obtained from the surrogate data with the empirical one is very high, and most often larger than the similarity of the original co-presence matrix. For networks aggregated over the whole data collection, Fig. 5 shows the distributions of degrees and of weights (see also SI). The first sampling method leads to an overestimation of degree values (resulting in a shift of the distribution), the second method tends to shift the distribution to lower degree values (except for the conference case), and the third method yields context-dependent over-or under-estimations of degree values. Note that the distributions of degrees of the copresence networks are not shown in the figure as the degree values are very strongly overestimated. Distributions of weights (aggregated contact durations) recover well the ones of the data for all sampling methods, and are closer than the ones of the co-presence networks. For each data set we compute the cosine similarity between the neighbourhoods of each nodes from each daily network, averaged for all nodes and all pairs of daily networks. The neighbourhood of a node n is defined as the vector of the link weights between n and every other nodes (if the link does not exist the weight is set to zero). We compare the values obtained for the contact data, the co-presence data, and for the networks generated by each sampling method of the co-presence data, averaged over 100 realisations for each sampling method. For reference, we also compute as null model the average similarity when links in the contact data are shuffled randomly within each daily network.

To investigate intermediate timescales of aggregation, Table 4 quantifies the similarity between networks aggregated in different days. The measure is defined as the average cosine similarity between all pairs of instances of a node's neighbourhood, averaged over all nodes. We see that the similarity is higher for the co-presence networks, as expected since the networks are denser. The sampling method 1 generates networks that are more similar than the data, and the other two methods generate networks that are less similar (with the exception of the LH10 case, and the method 2 in the InVS13 case). In the cases of the method 2 for the LyonSchool data, and the methods 2 and 3 for the Thiers13 data, the sampled networks are even almost as different as they would be after a random shuffling of the links.

In addition, Fig. 6 gives the evolution of the average degree and strength for networks aggregated in increasingly long time windows. First, the evolution of the real average aggregated strength is usually better recovered than for the degree by the various sampling better. Second, which sampling method recovers better the evolution of the degree is again context dependent. However, in all cases the sampled data are much closer to the contact data than the co-presence network, which overestimates very strongly these quantities. In a network, more "central" nodes are usually considered as important, as they might play an important role for instance in spreading processes (or other dynamical phenomena) occurring in the network. It is thus of interest to understand whether the most central nodes in the contact network can be identified either in the raw co-presence data or in the surrogate contact data built from the co-presence information. As there are several ways of determining central nodes in a network, we consider here three of the most well-known centrality measures and apply them to the networks aggregated over the whole data collection: degree k, strength s and betweenness b of nodes in the aggregated networks. For each instance of each sampling method, we thus build the resulting surrogate aggregated contact network and rank nodes according to each centrality measure. We then compute the Jaccard similarity index between the top N % nodes in the real contact network and in the surrogate one. We plot in Fig. 7 the median similiarity with the 90 % confidence interval, as a function of N , for the InVS15 case (see SI for the other cases). 28.4 (0.360) 42.6 (0.559) 15.8 (0.117) For each dataset we compute the maximum coreness, and report between parenthesis the Jaccard index between the k-core of the contact network and the k-core in the original and sampled co-presence data (results are averaged over 100 realisations for each sampling method).

Node centralities

In general, no sampling method recovers correctly the most central nodes for low values of N . The best results are obtained for the conference data with similarities around 0.2 -0.4. The similarity values increase as N increases but reach most often only values of ∼ 0.5 when considering the top 50 % nodes, meaning that only 25 % of the most central nodes are identified when using the surrogate data. The best results are obtained for the first sampling method for the LyonSchool case and for the LH10 case, with similarities reaching 0.6 -0.7. Results are typically better than the random baseline but do not outperform the detection of most central nodes based on the whole co-presence network. In terms of the most central nodes as defined by the k-core decomposition (we recall that the k-core of a network is the maximal subgraph such that all nodes in the subgraph have at least degree k, and k is called the coreness), the overestimation of degrees in the co-presence network leads to an overestimation of the maximum coreness, while sampling leads to values closer to the ones of the contact data, but once again in a context-dependent way. The maximum core itself is only partially recovered in the whole and in the sampled co-presence networks (see Table 5).

Using surrogate contact data in epidemic simulations

We have seen in the previous section that none of the three sampling methods yields a perfectly accurate description of all the relevant features of the true contact network: each sampling method yields surrogate data with both interesting similarities and potentially important discrepancies with respect to the original contact data. We now consider the issue of using such surrogate data in simulations of spreading processes: as precise data on face-to-face contacts is not always available, it is important to understand if co-presence information can allow us to obtain on the one hand an accurate prediction of the outcome of an epidemic process, and on the other hand a reliable estimation of the impact of containment measures. In particular, it is important to be able to classify potential containment strategies to determine which one(s) are most adequate.

To this aim, we consider the paradigmatic Susceptible-Infectious-Recovered (SIR) model for epidemic spreading. In this model, susceptible (S) individuals can become infectious (I) at rate β when in contact with an infectious node. Infectious nodes recover spontaneously at rate µ and enter an immune recovered (R) state. Simulations start with a single infectious individual chosen at random and carried out until there are no infectious individuals left in the population, i.e., individuals are either still susceptible or have been infectious and have then recovered. The impact of the epidemics is then quantified by the final fraction n i of individuals in the R state.

We set β = 0.0004 and vary µ by tuning the reproductive number R = β/µ. For each value of R, we measure the fraction P (n i > 20 %) of "large" outbreaks in which the fraction n i of the population that was reached by the outbreak is at least 20 % and the distribution of the sizes n i of these large outbreaks. We average the results over 10 000 simulations performed on the empirical contact network. For each sampling method, we build 100 different instances of the surrogate contact network, and perform 100 simulations on each surrogate network.

We also consider several simple methods to mitigate the spread, namely the vaccination of a number of individuals in the population, under the assumption of a perfect vaccine efficiency: vaccinated individuals cannot become infectious nor transmit the disease and thus slow down and hinder the propagation. We consider the vaccination of (i) 5, 10 or 20 individuals chosen at random (ii) the most central 5, 10 or 20 individuals, where centrality is measured according to either degree, strength or betweenness in either the real or surrogate contact networks (iii) when the population is structured in groups, the vaccination of all individuals in one group.

Figures 8 and9 summarize our results for the InVS15 dataset (see SI for the figures obtained with the other datasets). In terms of the evaluation of the impact of a spreading process, results are context dependent. The simulations performed on the surrogate data obtained with the first method generally lead to an overestimation of the epidemic risk, except for the hospital data. When using the second sampling method, we obtain a good estimation of the risk for the conference, school and highschool data but an underestimation for offices and hospital data. The third method on the other hand leads to a correct estimation for the offices and hospital data but an underestimation for the school and highschool and an overestimation for the conference.

We show in Fig. 9 the impact of the various vaccination strategies, quantified through the ratio of the probabilities of large outbreaks with and without vaccination, as well as the ratio between the median sizes of these large outbreaks. We rank the strategies according to their efficiency in the real contact network, in order to visualize easily whether the surrogate networks lead to the same classification of the We plot the ratio between the vaccination and no vaccination cases of the fraction of the total number of outbreaks that reach at least 20 % of the population (top), and of the median size of these outbreaks (bottom) for different vaccination strategies, for the original data and the reconstructed networks. The vaccination strategies are ordered by decreasing efficiency, based on the effect on the real contact data. The group * strategies consist in vaccinating one or several groups entirely; the group rand strategy vaccinates ng random nodes, where ng is the average group size; the rand n strategies randomly vaccinates a specified fraction n of nodes; the b n, k n, s n strategies vaccinate the top n % nodes according to, respectively, betweenness centrality, degree and stregnth ranking.

strategies: indeed, even when the impact of each specific strategy is not accurately quantified, it would be interesting at least to understand which methods are most efficient. Results are once again uneven and context dependent (see also Table 6). In several cases such as SFHH the ranking of strategies obtained from the sampled co-presence is overall respected (Kendall's tau of 0.818 for the sampling method 1 on the size of outbreaks), while it can be strongly reshuffled in other cases (for instance in the Thiers13 case).

Discussion and conclusion

In this paper, we have investigated whether low resolution co-presence information can be used as a substitute for detailed face-to-face proximity data, both from the point of view of extracting large-scale structural and statistical features of the temporal contact network in a population and in data-driven models of epidemic processes in a population. We have considered several data sets collected in various contexts that contain both high-resolution data on face-to-face contacts between individuals and a coarser location data, both with temporal resolution. The location data can thus be transformed into a co-presence temporal network between individuals. Given its lower spatial resolution, this co-presence data contains much more events than the contact data, leading to much denser aggregated networks: indeed, all individuals in a given area are considered as co-present, while only some of them are typically engaged in a face-to-face contact. Despite this expected issue, a number of properties related to group structure and statistical distributions of temporal properties are similar in contact and co-presence data, with similar matrices of densities of links between groups and broad distributions of (aggregate) contact durations.

We have thus examined several methods to downsample the co-presence networks to create surrogate contact networks with overall the same amount of contact time than the real contact data. The surrogate data statistics are in general closer to the real contact data than the raw co-presence, in particular regarding the distribution of node degrees and link weights (and their evolution in networks aggregated over increasing time windows). These results mean in particular that the distribution of aggregate contact durations, a very important property that has a strong impact on the unfolding of processes on networks such as epidemic processes, could be approximately retrieved from simple sampling processes of the co-presence data and thus fed into data-driven models of populations. Several other properties, such as precise value of the average degree, average clustering or size of largest cliques and cores, turn out however to be strongly context-dependent. Moreover, the most central nodes of the contact network are not better identified than using the bare co-presence information.

We have moreover investigated the use of such surrogate contact data in numerical simulations of spreading processes in a population. Overall, simulations performed on surrogate data obtained with one of the sampling method yield results close to the ones obtained with the real data, while the other methods over-or underestimate these results, but the best method turns out to depend on context (Note however that all these methods give obviously results much closer to the one of the real contact network than if raw co-presence is used, given co-presence overestimates strongly the contacts and thus yields a strongly overestimated epidemic risk). We moreover investigated the possibility to rank containment strategies according to their efficiency, and found that this ranking is once again context dependent: in some cases, simulations on sampled co-presence networks allow us to uncover the most efficient vaccination strategies for containing a spread on the real contact data, while in other cases the rankings differ quite strongly.

In conclusion, we showed that co-presence data, while yielding interesting insights into some of the large scale properties of the contact network, is not easily usable to build in a reliable and systematic fashion surrogate contact data that reproduces detailed features of the real contacts and could be used in numerical simulations to predict the outcome of spreading processes and the impact of containment strategies, at least for processes involving contagion at short distances [START_REF] Stopczynski | Physical proximity and spreading in dynamic social networks[END_REF] (note that, while more sophisticated sampling procedures might be devised, they would most probably involve more parameters and/or more additional information not present in the raw co-presence data, and would also most probably still give context-dependent results). We note however that even coarse location information has been shown to be a useful additional information whenever the precise contact data is incomplete [START_REF] Sapienza | Estimating the outcome of spreading processes on networks with incomplete information: a mesoscale approach[END_REF]. Optimally, data collection with wearable sensors should thus contain both high resolution data about relative positions of individuals, in order to detect face-to-face proximity, and coarser co-presence information to inform for instance on mobility patterns within buildings or complement potential data losses.

Figure 1

 1 Figure1Temporal distributions for co-presence and contact events -InVS15. We show for both the contact and co-presence of the same data set the distributions of event and inter-event duration, link weights (as total contact duration) and number of contacts per link.

Figure 2

 2 Figure 2 Contact and co-presence matrices -InVS15. Comparison between the average matrices of link density for the contacts and the co-presence daily aggregated networks. Values are averaged over all days of the data collection. Both plots have the same colour scale.

Figure 3

 3 Figure 3 Day 2 contact and co-presence networks -InVS15. This figure highlights the difference in terms of link density when comparing contact and co-presence daily networks. Different node colours correspond to the different administrative departments.

Figure 4

 4 Figure 4 Number of contacts as a function of the number of individuals present -InVS15. We plot the number of contacts ca(t) occurring at time t in a certain area a as a function of the number of individuals na(t) present at the same time in a. The red line shows the median of the scatter plot, with error bars defined by the 25 % and 75 % percentiles.

Figure 6

 6 Figure 6 Evolution of the mean aggregated degree and strength. InVS15 data set We compare how the average degree k(t) and the average strength s(t) grow as we aggregate the network on increasing time-windows, for the real contact data and each sampling method.

Figure 7

 7 Figure 7 Node ranking similarity. InVS15 data set We plot for each co-presence sampling method the Jaccard similarity between the top N % nodes in the real and surrogate contact data, when ranked according to their degree k, their strength s or their betweenness centrality b vs. N . The plot shows the median similarity and the shaded areas give the 90 % confidence interval.

Table 1

 1 Characteristics of the data sets. Np is the number of participants, Na the number of RFID readers, T the total duration of the data collection.

	Data set	Location	Year	Np	Na	T	Ref
	InVS13	Fr. Health Obs.	2013	92	27	2 weeks	[31]
	InVS15	Fr. Health Obs.	2015	232	45	2 weeks
	LH10	Hospital	2010	81	8	3 days	[32]
	LyonSchool	Primary school	2009	242	15	2 days	[33]
	SFHH	Conference	2009	403	12	2 days	[34]
	Thiers13	High school	2013	326	18	1 week	[35]

Table 2

 2 Similarity between contact matrices.

		InVS13 InVS15	LH10	LyonSchool Thiers13
	Co-presence	0.790	0.710	0.968	0.706	0.681
	Sampling 1	0.946	0.829	0.960	0.845	0.857
	Sampling 2	0.958	0.901	0.894	0.945	0.937
	Sampling 3	0.888	0.816	0.958	0.738	0.691

Table 3

 3 Characteristics of the contact, co-presence, and sampled co-presence networks.

		InVS13 InVS15	LH10	LyonSchool SFHH	Thiers13
	kc	2.9	6.4	14.0	47.3	28.8	13.5
	k	20.9	35.0	18.2	194.5	234.3	126.8
	k1	5.8	14.2	14.4	101.3	116.7	52.2
	k2	0.9	3.6	7.7	36.9	45.2	5.2
	k3	5.3	5.0	14.0	21.2	40.4	4.1
	ρc	0.030	0.028	0.175	0.196	0.072	0.041
	ρ	0.211	0.152	0.227	0.807	0.807	0.383
	ρ 1	0.058	0.061	0.179	0.420	0.290	0.158
	ρ 2	0.009	0.016	0.097	0.153	0.112	0.016
	ρ 3	0.054	0.022	0.175	0.088	0.101	0.013
	ωc	4.4	7.6	14.3	22.5	11.0	9.4
	ω	18.8	38.7	22.7	141.0	*	74.6
	ω 1	6.6	10.3	17.2	41.5	34.7	33.8
	ω 2	3.0	5.3	8.6	12.8	12.2	3.9
	ω 3	5.5	4.8	17.1	6.3	9.0	3.8
	cc	0.178	0.239	0.428	0.520	0.260	0.379
	c	0.417	0.409	0.491	0.868	0.880	0.581
	c1	0.255	0.266	0.432	0.596	0.442	0.586
	c2	0.045	0.139	0.309	0.370	0.212	0.092
	c3	0.205	0.101	0.426	0.193	0.161	0.047

Table 4

 4 Average similarity between daily networks.

		InVS13 InVS15	LH10	LyonSchool Thiers13
	Contact	0.333	0.305	0.351	0.643	0.431
	Co-presence	0.415	0.348	0.449	0.806	0.683
	Sampling 1	0.361	0.344	0.436	0.749	0.515
	Sampling 2	0.388	0.271	0.403	0.175	0.084
	Sampling 3	0.286	0.205	0.437	0.042	0.071
	Null model	0.022	0.010	0.061	0.046	0.010

Table 5

 5 Comparison of the maximum k-core properties.

		InVS13	InVS15	LH10
	Contact	11	25	23
	Co-presence	78 (0.607)	112 (0.681)	32 (0.682)
	Sampling 1	22.5 (0.660) 57.8 (0.719) 26.0 (0.683)
	Sampling 2	5.23 (0.479) 17.4 (0.692) 16.4 (0.655)
	Sampling 3	28.1 (0.591) 27.4 (0.639) 25.5 (0.693)
		LyonSchool	SFHH	Thiers13
	Contact	47	33	24
	Co-presence	181 (0.615)	320 (0.522)	210 (0.684)
	Sampling 1	99.8 (0.638)	111 (0.575)	76.7 (0.640)
	Sampling 2	39.7 (0.501) 41.4 (0.555) 17.7 (0.375)
	Sampling 3

  Figure 8 Epidemic prevalence. InVS15 data set. We plot the fraction of the total number of outbreaks that reach at least 20 % of the population (crosses) and the distribution of the sizes of these outbreaks (boxplots) for several values of the reproductive number R, for the original and the surrogate contact networks.

			1.0
	i > 0. 2)	0.8
	Epidemic size (n	0.2 0.4 0.6															Data Sampl. 1 Sampl. 2 Sampl. 3
			0.0		1.1	1.5			2.0 R			2.5		3.0 P(ni > 0. 2)
		1.2
	no vacc	0.8 1.0
	vacc /	0.4 0.6			Data
	P	0.2			Sampling 1 Sampling 2
		0.0 1.1	group_SCOM_SSI	group_SFLE	group_DCAR Sampling 3 group_DMCT rand_5 group_DG_SDOC_SRH	group_DMI	group_rand	group_DSE	group_DISQ	group_DST	rand_10	b_5	k_5	s_5	rand_20	b_10	k_10	s_10	b_20	s_20	k_20
		1.0
	vacc	0.9
		0.8
	n vacc /	0.6 0.7			Sampling 1 Data
		0.5			Sampling 2
		0.4	group_SFLE	group_SCOM_SSI	group_DG_SDOC_SRH Sampling 3 group_DMCT rand_5 group_DCAR	group_DSE	group_rand	group_DMI	group_DISQ	group_DST	rand_10	b_5	k_5	s_5	rand_20	k_10	b_10	s_10	b_20	s_20	k_20
	Figure 9 Vaccination strategies. InVS15 data set.

Table 6

 6 Comparison of the vaccination strategy rankings. For each sampling method we compute Kendall's tau between the list of vaccination strategies ranked by increasing efficiency for the contact data and for the sampled co-presence networks, both in terms of the fraction of large outbreaks (P ) and of median sizes of the large outbreaks (n).

		InVS13	InVS15		LH10
		P	n	P	n	P	n
	Sampling 1	0.515	0.235	0.377 0.766	0.397	0.412
	Sampling 2	0.324	0.382	0.377 0.481	0.456	0.412
	Sampling 3	0.279	0.235	0.394 0.706	0.324	0.426
		LyonSchool	SFHH		Thiers13
		P	n	P	n	P	n
	Sampling 1	-0.012 -0.051 0.485 0.818	0.048	0.299
	Sampling 2	0.091	0.083	0.758 0.485 -0.074 -0.108
	Sampling 3	0.020	0.162	0.636 0.545	0.299	-0.108

Availability of data and material

The datasets supporting the conclusions of this article are available on the SocioPatterns website: http://www.sociopatterns.org/datasets/ (see also references in this paper).

Competing interests

The authors declare that they have no competing interests.

Author's contributions MG and AB conceived and designed the study. MG performed the numerical simulations and statistical analysis, created the figures and wrote the first draft of the manuscript. MG and AB wrote the final version of the manuscript.

Author details 1 Aix Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France. 2 GESIS, Leibniz Institute for the Social Sciences, Unter Sachsenhausen 6-8, 50667 Köln, Germany. 3 Data Science Laboratory, ISI Foundation, Torino, Italy.


R Fares
email: rfares@oact.inaf.it

V Bourrier

A A Vidotto

C Moutou

M M Jardine

P Zarka

Ch Helling

A Lecavelier Des Etangs

J Llama

T Louden

P J Wheatley

D Ehrenreich

MOVES I. The evolving magnetic field of the planet-hosting star HD189733

Keywords: stars, magnetic fields -stars, individual, HD189733 -techniques, polarimetry -Planet-Star Interaction

MOVES -I. The evolving magnetic field of the planet-hosting star HD189733

INTRODUCTION

Hot-Jupiters (HJs), i.e., gas giant exoplanets orbiting close ( 0.1 au) to their host stars, are useful laboratories to study the complex interactions (e.g. magnetospheric, tidal, ionisa-lier des Etangs 2013; Koskinen et al. 2013), which can be probed through transmission spectroscopy technique (Vidal-Madjar et al. 2003;Fossati et al. 2010;Lecavelier Des Etangs et al. 2010;Ehrenreich et al. 2012;Lecavelier des Etangs et al. 2012;Bourrier et al. 2013). In addition to the extreme radiation environment, HJs are also immersed in stellar winds, whose physical characteristics are unparalleled to those felt by the planets in our own solar system (Preusse et al. 2005;Grießmeier et al. 2007;Vidotto et al. 2009Vidotto et al. , 2015)). The extreme particle and magnetic environments of the stellar wind could lead to powerful reconnection events between the stellar and planetary magnetospheres (e.g. Ip, Kopp & Hu 2004), which may have the potential to enhance the stellar activity (Cuntz, Saar & Musielak 2000;Shkolnik, Walker & Bohlender 2003;Shkolnik et al. 2008;Scandariato et al. 2013). In addition, as a result of the interaction between the planet's magnetosphere and the coronal material of the star, bow shocks might also form around HJs (e.g. Vidotto, Jardine & Helling 2010;[START_REF] Cohen | [END_REF]Llama et al. 2013;Matsakos, Uribe & Königl 2015), explaining asymmetric UV transit features (Fossati et al. 2010;Vidotto, Jardine & Helling 2010). Radio emission resultant from the interaction between the stellar magnetised wind and a magnetised exoplanet is expected to be several orders of magnitude larger than that of the largest emitter in our own solar system, Jupiter (Zarka 2007;Grießmeier, Zarka & Spreeuw 2007;Jardine & Cameron 2008;Fares et al. 2010;Vidotto et al. 2012Vidotto et al. , 2015;;See et al. 2015). A varying stellar magnetic field has also implication for the atmosphere of the orbiting planet leading to variability of the Cosmic Ray flux reaching the planetary atmosphere (Helling et al. 2013;Rimmer & Helling 2013). Cosmic rays have been studied as one example for external ionisation sources that also effect the chemical compositions by possibly opening reaction paths to carbo-hydrate molecules (Rimmer, Stark & Helling 2014;Rimmer, Helling & Bilger 2014).

HD189733 is an ideal system to study these different types of interactions. At a distance of 19.3 pc, the bright (V=7.7) and active K2V star hosts a transiting HJ orbiting at 0.031 ± 0.001 au, i.e., at an orbital distance of 8.84 ± 0.27 R (see Table 1). The interactions within this system have been the subject of many studies in the literature (e.g. Smith et al. 2009;Pillitteri et al. 2010;Berdyugina et al. 2011;Lecavelier Des Etangs et al. 2011;Ben-Jaffel & Ballester 2013;Bourrier & Lecavelier des Etangs 2013;Bourrier et al. 2013;Poppenhaeger, Schmitt & Wolk 2013;Louden & Wheatley 2015;Brogi et al. 2016;Bott et al. 2016; Barnes et al. 2016). Because of its active host, HD189733b is in fast-changing radiation, particle and magnetic environments. Recent UV observations of the planetary transit showed that the properties of the hydrogen exosphere surrounding the planet are varying over time (Lecavelier des Etangs et al. 2012). An X-ray flare was detected 8h prior to the UV detection of atmospheric escape, suggesting that variations in the XUV emission of the host-star and/or its magnetised outflowing plasma might have been the cause of the observed variability.

Variations in the quiescent stellar wind are also expected for this system. Through spectropolarimetric observations performed by Moutou et al. (2007) and Fares et al. (2010), it was shown that the large-scale magnetic field of HD189733 varied for the observed epochs. Recently, using et al. (2009) the large-scale magnetic maps reconstructed by Fares et al. (2010), Llama et al. (2013) showed that the stellar wind of HD189733 presents inhomogeneities at the position of the planet, which could cause short timescale variations in the UV transit lightcurve (on the order of an orbital period), as well as longer timescale variations (on the order of a year, due to the change in the stellar magnetic field). These authors concluded that multi-wavelength data acquired simultaneously would provide the best tool for a comprehensive characterisation of the system. In this context, we started a multi-wavelength observational campaign of this system, in the frame of the MOVES collaboration (Multiwavelength Observations of an eVaporating Exoplanet and its Star, PI V. Bourrier). Observations of the star and the planet were obtained at similar epochs with ground-based and space-borne instruments in X-rays with Swift and XMM-Newton, UV with HST and XMM-Newton, optical spectropolarimetry with NARVAL and ES-PaDOnS, and radio with LOFAR. In the present (first) paper of this collaboration, we use optical spectropolarimetric observations to reconstruct the surface magnetic field of HD189733 at five different epochs, contemporaneous to other sets of observations. These magnetic field maps will provide a crucial input to the analysis, modelling and interpretation of the multi-wavelength data sets that will follow in forthcoming papers.

This paper is organised as follows. Section 2 presents our observations. In Section 3, we describe the magnetic imaging method we use. The results are shown in Section 4, where we present the reconstruction of the magnetic field of HD189733 at five different epochs. We discuss the magnetic field evolution of HD189733 in Section 5, based on this paper's results and results of Moutou et al. (2007) and Fares et al. (2010). Section 6 presents the summary and conclusions of this work.

Our spectropolarimetric data were obtained using both NARVAL spectropolarimeter at the TBL (2m) and ES-PaDOnS at CFHT (3.6m). A spectropolarimeter is a spectrograph with a polarisation section, allowing measurement of the polarisation of the spectral lines. A circular polarisation spectrum is extracted from 4 subexposures taken each at a different angle of the polarisation waveplates. ES-PaDOns and NARVAL are twin instruments, both operating in the optical domain (370 to 1000 nm) at a resolution of 65,000 in the polarisation mode. Data reduction is done automatically using Libre-Esprit, a fully automated reduction tool installed at TBL and CFHT (Donati et al. 1997). The spectra are normalised to a unit continuum, their wavelengths referring to the Heliocentric rest frame. Telluric lines are used to correct from spectral shifts due to instrumental effects. This correction secures a radial velocity (RV) precision of about 15 m s -1 (Moutou et al. 2007;Morin et al. 2008b).

Observations with both NARVAL and ESPaDOnS were obtained in service mode. Our NARVAL 2013 data (PI Bourrier) were collected as follow: 5 spectra in May (11-13), 2 spectra in June (12-15), 5 spectra in July (01-04-05-08-14), 14 spectra in August (02 to 23), 13 spectra in September (02 to 24) and finally 5 spectra in October (08-09-12-13-18). ES-PaDOnS' data were obtained via a filling program targeting planet-hosting stars (PI Moutou), 7 spectra were collected in September 2013. Table 2 presents the log of observations of all our data in 2013. In 2014, 15 spectra were collected between 01 September and 20 October. Another 18 spectra were collected between 20 May 2015 and 20 July 2015. The log of the 2014 and the 2015 (PI Bourrier) observations are presented in Table 3.

The equatorial rotation period of HD189733 is ∼ 12 days (Fares et al. 2010). For the sake of homogeneity with our previous studies of this star, we used the same ephemeris as in Moutou et al. (2007) and Fares et al. (2010). Rotational (ERot) and orbital (E Orb ) cycles are calculated by: T0 = HJD 2, 453, 629.389 + 12 ERot T0 = HJD 2, 453, 629.389 + 2.218575 E Orb

(1)

For cool stars, the polarisation signature in single lines is typically within the noise level. Using the polarisation information of many spectral line simultaneously, we can extract the polarisation signal of the spectrum, this is known as multi-line technique Least-Square Deconvolution (LSD, Donati et al. 1997). LSD assumes that all lines have the same polarisation information, and extracts the polarisation signature by deconvolving the observed spectrum with a line mask. LSD calculates intensity (Stokes I) profiles, circular polarisation (Stokes V) profiles, as well as a null polarisation profiles (labelled N). These profiles are calculated using a combination of the sub-exposures taken at different angles of the waveplates. The Null profile is a check profile, it is calculated such to cancel out the stellar polarisation signature, and thus should contain no polarisation. It helps check for spurious or instrumental signatures. For more details on how these profiles are calculated, see Donati et al. (1997) and Mengel et al. (2017).

We compute a line mask for HD 189733 using Kurucz's lists of atomic line parameters (Kurucz CD-Rom 18) and a Kurucz model atmosphere with solar abundances, temperature set to 5000 K and logarithmic gravity (in cm s -2 ) set to 4.0. Only moderate-to-strong lines, featuring central depths larger than 40% of the local continuum, are included in the mask (before any macroturbulent or rotational broadening); the strongest and broadest features (such as Balmer lines) are excluded. In the optical domain, this mask contains about 4000 lines. The LSD profiles calculated by deconvolving the stellar spectra with the mask have a S/N ∼ 30 times higher than the S/N in single lines with average magnetic sensitivity (see Tables 2 and3), allowing for the detection of the polarisation signature.

We calculated the RV of each Stokes I profile by fitting a Gaussian profile to it. Since the star hosts a hot-Jupiter, the RV signatures varies over the planetary orbit. Our RV data agree (within the error bars) with the orbital solution of Boisse et al. (2009). There is an offset between our values and theirs, +0.06 km s -1 for June-July 2013, +0.13 km s -1 for August 2013, +0.10 km s -1 for September 2013, +0.13 km s -1 for September 2014, and +0.11 km s -1 for June 2015. Such offsets are due to the use of different reduction pipelines (ESPaDOnS and Narval vs Sophie). For the tomographic Imaging, we shift each profile by its RV to centre all profiles around 0 km s -1 .

MODEL DESCRIPTION AND IMAGING METHOD

Magnetic fields, when present in the photosphere, cause splitting of the magnetically sensitive spectral lines via the Zeeman effect. The lines that form in such environments are polarised. Measuring the wavelength shift between spectral components of the line (when possible, which depends, among others, on the amplitude of this shift, the rotational broadening and the spectrograph resolution) allows us to calculate the longitudinal magnetic field of the star (see, e.g., Shulyak et al. (2014) for M dwarfs). The field topology (i.e. distribution, polarity, configuration), on the other hand, can not be determined by the wavelength shift alone, but is determined using a tomographic imaging technique, Zeeman-Doppler Imaging ZDI (Donati et al. 1997). Polarisation of the spectral lines depends on the orientation of the magnetic field relative to the line of sight (see Fig 2 of Reiners 2012). To map the field, spectra are collected during one or more stellar rotations. ZDI inverts these spectra/profiles into a magnetic topology that can produce the observed polarisation signatures. Because this problem is ill-posed, ZDI uses Maximum-Entropy regularisation to get the simplest magnetic map compatible with the data.

In the newest version of ZDI (Donati et al. 2006), the field is described by its poloidal and toroidal components [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF], with all the components expressed in terms of spherical harmonics expansion. The highest degree of spherical harmonic expansion used to map the field represents the ZDI map resolution around the equator. For slow rotators, such as HD189733, we limit the spherical harmonics to the lowest degrees (l ≤ 5 , see Fares et al. 2010 for more details). ZDI follows an iterative procedure, it compares synthetic Stokes V profiles to the observed profiles col-Table 2. List of observations in 2013. The columns list, respectively, the dates of observation, Heliocentric Julian Date and UT time of observations (at mid-exposure), the peak S/N (per 2.6 km s -1 velocity bin) of each observation (around 700 nm), the rotational cycle of the star and orbital cycle of the planet calculated using the ephemeris of Eq. 1, and the radial velocity (RV) of the star at each exposure, the rms noise level (relative to the unpolarized continuum level Ic and per 1.8 km s -1 velocity bin) in the circular polarisation profile produced by LSD. The data were taken using NARVAL spectropolarimeter, except for 6 spectra collected using ESPaDOnS (marked by * next to the date of observation). The exposure time of each observation is 4×900 s. Dates marked with † were not used for the mapping of the stellar field (see text for more details). lected during the stellar rotational cycles. In practice, the stellar surface is divided into 9000 grid cells of similar area. The contribution of each grid cell is then calculated in the weak field regime, and a synthetic Stokes profile for each observed rotation phase is produced.

Date (UT

In addition to the field intensity and distribution, ZDI also gives an indication on the stellar inclination (up to ∼ 10 • accuracy), on the stellar rotation period (see, e.g. Alvarado-Gómez et al. 2015), and on the differential rotation (DR). We describe DR using a solar-like DR law, where the rotation rate at a latitude θ is defined by Ω(θ) = Ωeq -dΩ sin 2 (θ), where Ωeq and dΩ are respectively the rotation rate at the equator, and the difference in rotation rate between the pole and the equator. In practice, to measure DR, we reconstruct the magnetic image at a given magnetic energy for a pair of fixed (Ωeq, dΩ), and obtain the reducedchi squared (χ 2 r ) of the fit. We investigate the χ 2 r values of the 2D parameter space of (Ωeq, dΩ). The optimum DR parameters are the ones minimizing χ 2 r . They are obtained by fitting the surface of the χ 2 r map with a paraboloid around the minimum value of χ 2 r .

The Null profiles, for each epoch, are used as a test to check for spurious polarisation. In practice, we fit these profiles with a zero magnetic field configuration. Since the Null profiles should contain no polarisation, the χ 2 r of the fit should be one or less. A χ 2 r greater than one indicates either a spurious polarisation signature, or an underestimation of the error bars. If a systematic signature is found in the profiles, it indicates a spurious origin of the signature. We calculate a mean signature of the Null profiles, and subtract it from the Stokes V profiles to eliminate this spurious feature.

RESULTS

ZDI requires the use of data covering one or more stellar rotations to derive the photospheric magnetic map. In some field configurations (see, e.g. Morin et al. 2008a) the largescale magnetic field of the star is stable for many stellar rotations, and also over many years. In these cases, one can combine data collected over many rotations together as one dataset. To estimate whether the field is stable over many stellar rotations, we compare the quality of the fit and the shape of the polarisation profiles at the same rotational phases. Our observations cover many stellar rotations at each observing epoch, from June 2013 to July 2015 (see Tables 2 and3). We performed a series of tests on the data, reconstructing the maps using a combination of datasets for each epoch. We find that using Stokes V profiles spread over more than two consecutive stellar rotations to reconstruct the map worsens the quality of the fit. We therefore adopted datasets of up to two stellar rotations for each reconstructed map. Some spectra were not used in the reconstruction because they were collected with a rotational phase gap of more than a stellar rotation in respect to the main dataset (see Tables 2 and3).

Differential rotation

Differential rotation distorts the magnetic regions on the stellar surface. Spectropolarimetric data can therefore be used to estimate the level of DR. We applied the same technique as in Fares et al. (2010) to search for differential rotation in our data (explained in Section 3). In this study, we were able to detect DR for August 2013 dataset. The χ 2 r map in the (Ωeq, dΩ) space is shown in Fig 1 . A well defined χ 2 r minimum is found for dΩ = 0.11 ± 0.05 rad d -1 and Ωeq of 0.535 ± 0.004 rad d -1 , corresponding to a rotation period at the equator 11.7±0.1d. The shear value dΩ is consistent with the one measured by Fares et al. (2010) at dΩ = 0.146 ± 0.049 rad d -1 . These measurements, within their uncertainties, do not reveal a variation of DR. HD 189733 is a slow rotator for which the Fourier Transform of the Intensity profile technique presented by Reiners & Schmitt (2002) can not be applied. There are no measurements of the DR of this star in the literature, apart from the recent work of Cegla et al. (2016). Cegla et al. (2016) modelled the Rossiter-McLaughlin effect to probe planetary parameters and stellar differential rotation, and found a dΩ > 0.12 rad d -1 for HD 189733, in agreement with our findings. Stellar differential rotation is not constant across spectral types and rotation rates. Many observational and theoretical studies have addressed the question of its variation (see, e.g., Barnes et al. 2005;Reiners 2006 2016) propose an empirical relation between the shear value, the effective temperature of the star and its rotation. For HD 189733, the predicted value of the dΩ is 0.076 +0.004 -0.001 . Our value is slightly higher, their value is, however, within the error bars of our measurement.

Magnetic maps

We reconstruct the magnetic maps for five observing epochs. The observed Stokes V (in red) and a 1-sigma error-bar are shown in Fig 2, along with the ZDI reconstructed circular polarisation (in black). The observed profiles show a detection of the polarisation signature (false-alarm probability less than 10 -5 , see Donati et al. (1997) for the classification of the detections). The reconstructed maps are shown in Fig 3 . As mentioned previously, each dataset used to reconstruct a magnetic map consisted of no more than two stellar rotations.

• June/July 2013

We reconstructed the map of June/July 2013, although our phase coverage is not optimal. The χ 2 r of the reconstruction is 1.15. The mean magnetic field is of 36 G. 61% of the magnetic energy is in the toroidal component. For the poloidal component, we define axisymetric modes as having m = 0 and m<l/2 (m and l being the order and the degree of the associated Legendre polynomial, used to describe the field as spherical harmonics expansion). About 40% of the poloidal field is axisymetric.

• August 2013

For August 2013, the data samples the stellar rotation well. However, the Null profile shows a systematic signature in the core of the profile. Fitting these Null profiles with no magnetic field configuration leads to a χ 2 r of 1.3. We calculated a mean profile for all the Null profiles. Subtracting the mean Null profile from the August 2013 Null spectra leads to a χ 2 r of 0.7 when fitting the corrected Null profiles with no magnetic configuration. Given that the signatures in the Null profiles show a systematic trend, it is more likely to be a spurious signature rather than an underestimation of the error bars. We subtract the mean Null profile signature of the Stokes V profiles. The map is reconstructed with a χ 2 r c 0000 RAS, MNRAS 000, 000-000

The evolving magnetic field of HD189733 7 of 1.85. The average surface magnetic field is 40 G. 50% of the total energy is in the toroidal component of the field. 2% of the poloidal energy is in the axisymetric modes, and mainly in modes with m = 0. The poloidal component is thus mainly non-axisymetric. Spherical harmonics modes with l>3 (i.e. modes higher than the dipole, quadrupole and octupole) contribute to ∼ 38% of the poloidal energy.

• September 2013

In September 2013, we have data from ESPaDOnS and NARVAL, covering two stellar rotations. We reconstructed a map using those data with a χ 2 r of 2.15. We notice a systematic difference in the Intensity profile depth between Narval's spectra and ESPaDOnS' spectra. The difference in depth is not significant enough to require different modelling of the intensity profiles between instruments. This difference is due to the normalisation of the spectra at the telescope.

The average surface magnetic field is 42 G. 59% of the energy in the toroidal component. The poloidal field is still mainly non-axisymetric (88%, see Table 4). About 50% of the poloidal component is in the octupolar mode, higher orders contribute to ∼ 45%.

• September 2014

In September 2014, the average magnetic field drops to 31 G (see Fig 3). The χ 2 r of the fit is 1.25. The field remains mainly toroidal, with 78% of the total energy stored in the toroidal components. The poloidal field is strongly non-axisymetric, and only 12% of the energy lies within axisymmetric modes.

• July 2015

In July 2015, the magnetic field, and in particular the radial component, has changed significantly relative to September 2014. The radial component changes polarity around the equator between September 2014 and July 2015. The poloidal component contributes to 15% of the total energy. The average magnetic field is of 37 G.

The characteristics of each reconstructed map are listed in Table 4. ZDI does not provide error bars for the reconstructed map. However, statistical errors can be calculated by varying the input parameters (e.g., v sin i, stellar inclination, Ωeq, dΩ) within their error bars. Error bars on field characteristics are calculated by comparing the characteristics of the reconstructed field for the best set of input parameters to the characteristics of maps reconstructed by varying those best input parameters within their error bars. We follow Mengel et al. (2016) for error bar calculations: we fix (Ωeq,dΩ) and vary v sin i within its error bars (2.75 km s -1 to 3.2 km s -1 , Table 1). We then do a new set of maps fixing v sin i and varying (Ωeq, dΩ) within their error bars. Error bars in table 4 represent the highest values of our procedure.

Extrapolation of the magnetic field to investigate the corona and the planetary orbit

The interactions between the stellar wind and the planetary magnetosphere might trigger planetary emission, such as radio emission or bow shock formation (see, e.g. Fares et al. 2010;Llama et al. 2013), and it might influence the ionisa-tion state of the planetary atmosphere (Rimmer & Helling 2013;Rimmer, Helling & Bilger 2014). In this Section, we examine the stellar magnetic field in the corona up to the planetary orbit using a potential field extrapolation of the surface magnetic field. For this purpose, we use the Potential Field Source surface (PFSS) code of Jardine, Collier Cameron & Donati 2002, originally developed for the Sun (van Ballegooijen et al. 1998) based on Altschuler & Newkirk (1969). The potential field extrapolation assumes that there is no electric current in the corona. The components of magnetic field in the corona are described using a spherical harmonics decomposition. This technique delivers satisfactory results when compared to wind modelling of the solar corona (Riley et al. 2006).

The PFSS model extrapolates the magnetic field considering two boundary conditions: the first one being that the field is purely radial at a surface called the Source Surface, the second boundary condition is the observed field geometry at the surface of the star. We assume that the Source Surface is at 3.4 R . Our radial magnetic maps are used as a boundary condition at the surface of the star, thus giving a realistic model of the radial field.

The extrapolated field in the stellar corona for the five maps presented in this paper are shown in Fig 4 . This figure shows how different surface field configurations produce different field line configuration in the corona. Since the magnetic field lines in the corona are not those of a very simple configuration, and since the planet and the star are not synchronised, the planet crosses different field configuration on its orbit, as well as from one orbit to the other (the planet crosses in front of the same stellar field configuration after one beat period (of rotation and orbital periods), rather than after one orbital period). We calculate, for each observing epoch, the footprints of field lines connecting the stellar surface to the position of the planet on its orbit . They are shown in Fig 5.

The PFSS model allows the calculation of the field value and orientation at each position in the corona, out to the Source Surface. Since the planet is outside the Source Surface, we proceed as follows to calculate the field at its position. First, we calculate the position of the sub-planetary point on the Source Surface, we then calculate the energy budget at this point. We remind the reader here that at the Source Surface, the stellar field is purely radial, the meridional and azimuthal components are negligible. Assuming that the magnetic flux is conserved over spherical shells from the source surface to the planetary orbit, we calculate the decay of the magnetic field between the Source Surface and the orbit. In Fig 6, we plot the field value at the planetary orbit for each epoch of observations (including June 2007 andJuly 2008). The maximum value the field reaches at the planetary orbit can change by 100% between different epochs, which supports the importance of simultaneous observations when studying Star-Planet Interactions. We investigate the effect the error bars on the maps could have on the calculated magnetic field at the planetary orbit. To do so, for each epoch, we use as boundary condition each of the maps calculated for a range of v sin i, Ωeq, and dΩ (see Section 4.2), extrapolate the magnetic field for each map, and calculate the magnetic field at the position of the planet. We find that the mean difference between the field values presented in this paper and those calculated for the the maps reconstructed for sets of v sin i, Ωeq, and dΩ is of the order of 3mG. The difference can reach up to 70% of the mean value over a small fraction of the orbit (for one epoch, all other epochs showed a smaller maximum difference value). This shows the robustness of our results: at the position of the planet, the variation of the stellar magnetic field is real and not affected by the error-bars on the maps. The data presented in this Paper (5 observing epochs spanning two years) shows that the field of HD 189733 can evolve over a few stellar rotations. Combining datasets spanning more than two stellar rotations systematically reduces the quality of the fit. Fares et al. (2010) have merged data obtained in 2006, spread over 5 stellar rotations. We revisited the summer 2006 data and found that the map was overfitted. We adopt the results of Moutou et al. (2007) in this paper.

Due to the 5-year gap in the observations, we can not investigate the presence of cyclic variations in the stellar magnetic field, in a similar way as those reported for the planet-hosting star τ Boo (Donati et al. 2008;Fares et al. 2009Fares et al. , 2013)). Nevertheless, an evolution in the stellar magnetic field intensity and topology can be seen. Variations in both the axisymmetric contribution to the poloidal field For all epochs (apart from June 2006), the toroidal component dominates over the poloidal one. Petit et al. (2008) suggest, studying a sample of solar-like stars, that the toroidal energy dominates over the poloidal one for stars with rotation periods less than ∼ 12 days. HD 189733, having an equatorial rotation period of 12 days, does not contradict their findings. On the other hand, Donati & Landstreet (2009) suggest that stars with Rossby number (Ro) < 1 develop toroidal fields. HD 189733 has a Ro = 0.403 (Vidotto et al. 2014). HD 189733 field's geometry is therefore compatible with the geometries observed for stars with similar masses and Ro numbers. The fraction of axisymmetric field is almost always less than 50%.

In addition, it is also interesting to compare the evolution of the stellar magnetic field within the 2013 data sets, that are separated by just 9 rotation periods. We note that little variation was seen for the toroidal component. The percentage of the contributors (dipole, quadrupole, octupole) to the poloidal field, on the other hand, has changed. The main contributor to the field (i.e. azimuthal component) does not change polarity. The radial component evolved, with negative and positive magnetic features appearing at the surface.

Epoch

Bmean Figure 6. The stellar magnetic field value at the position of the planetary orbit (at 0.031 au). Different colours represent different epochs of observation. This plot shows that the planet is in a non-homogeneous environment, and that this environment varies from one epoch to the other. The error bar shown here is a mean error bar value per orbit, valid for all epochs (see text for more details).

CONCLUSIONS AND FUTURE PLANS

This paper is part of the MOVES collaboration (Multiwavelength Observations of an eVaporating Exoplanet and its Star), which aims to characterise comprehensively the complex environment of the exoplanet HD189733b. Orbiting a bright and active K dwarf at short distance, this transiting hot-Jupiter has been subjected to many stellar and planetary atmosphere studies. The main objectives of MOVES are to probe the different regions of the extended planetary atmosphere, its interactions with the host-star, and their temporal variability. The wider set of multi-wavelength observations (X-ray with Swift and XMM-Newton, UV spectroscopy with HST and XMM-Newton, and radio observations with LOFAR) were taken contemporaneously with the magnetic field mapping presented here.

In this first paper, we presented a detailed spectropolarimetric study of HD189733 and studied the evolution of its magnetism. Stellar magnetism is an important ingredient in stellar evolution, and also has important effects on planets surrounding these stars. The star was observed at five epochs (Jul 2013, Aug 2013, Sept 2013, Sept 2014and Jun 2015), during which we also collected X-ray and UV observations (Wheatley et al, in prep). Using Zeeman Doppler Imaging, we reconstructed the magnetic maps of the star. With a strength up to 45 G, the magnetic field is dominated by the toroidal component at the five epochs. The toroidal component is mainly axisymmetric during all observing epochs. In contrast, the poloidal component is mainly non-axisymmetric. We will continue monitoring this system to study the magnetic evolution on time-scales longer than 2 years and look for a potential magnetic cycle. These reconstructed magnetic maps are crucial for analysing multiwavelength observations. They allow us, modelling the stellar wind in the corona and at the planetary orbit, to reconstruct the X-ray emission and irradiation of the planet, as well as the spatial distribution of X-ray that will be absorbed by the extended atmosphere of the exoplanet.

  ; Collier Cameron 2007; Augustson et al. 2012; Reinhold, Reiners & Basri 2013; Distefano et al. 2016; Balona & Abedigamba 2016). Using Kepler photometric data, Reinhold, Reiners & Basri (2013) and Balona & Abedigamba (2016) studied the DR for a large sample of stars. HD 189733 DR is within the range of the shear values of Kepler stars both as a function of temperature and as a function of rotation. Balona & Abedigamba (

Figure 1 .

 1 Figure 1. Variations of χ 2 r as a function of Ωeq and dΩ , derived from the modelling of the Stokes V data set for August 2013. The outer colour contour corresponds to a 3.5% increase in the χ 2 r , and traces a 3 σ interval for both parameters taken as a pair.

Figure 2 .

 2 Figure 2. Observed Stokes V of HD189733 in red lines, and the fitted Stokes V profiles using ZDI in black, for June-July 2013 (top left), August 2013 (top middle), September 2013 (top right), September 2014 (bottom left) and July 2015 (bottom right). The rotational cycle of each observation (as listed in Table 2) and 1σ error bars are also shown next to each profile.

5

  DISCUSSION: MAGNETIC FIELD EVOLUTIONIn this Section, we take advantage of previous spectropolarimetric observations of HD189733 to investigate the evolution of its magnetic field (intensity and topology) over a 9 year-timespan. The first reconstructed magnetic maps of HD189733, based on ESPaDOnS/CFHT data of June and August 2006, were presented byMoutou et al. (2007). Later

Figure 3 .

 3 Figure 3. The reconstructed maps of HD 189733 for June-July 2013 (top row) August 2013 (second row), September 2013 (third row), September 2014 (fourth raw) and July 2015 (bottom row). The maps are in a polar flattened projection down to latitudes of -30 • , the equator is represented by the bold circle. The radial, azimuthal and meridional field components are shown. The magnetic flux values are labelled in G. Radial ticks around each map indicate the rotational phase of our observations.

Figure 4 .

 4 Figure 4. The extrapolated magnetic field of HD 189733 for June-July 2013 (top left), August 2013 (top middle), September 2013 (top right), September 2014 (bottom left) and July 2015 (bottom right). White lines correspond to the closed magnetic lines, blue ones to the open field lines (reaching the source surface). The star is shown at the same rotational phase (0.5) to better visualise the differences in magnetic field topology at each observing epoch. The star is viewed almost equator on (∼ 5 • ), the inclination of the system on the sky (as seen from Earth).

Figure 5 .

 5 Figure 5. The radial magnetic field of HD 189733 for June-July 2013 (top left), August 2013 (top right), September 2013 (middle left), September 2014 (middle right) and July 2015 (bottom). White dots represent the footprints of the field lines connecting the stellar surface to the position of the planet on its orbit.

Table 4 .

 4 Magnetic field characteristics of HD189733 for different epochs. The columns are: the epoch of the observations, the mean magnetic field at the surface of the star, the percentage of the toroidal energy relative to the total one, the percentage of the energy contained in the axisymmetric modes of the poloidal component relative to the poloidal energy, the percentage contribution of the dipolar, quadrupolar and octupolar components to the poloidal energy, and the mean stellar field at the position of the planetary orbit (see Section 4.3). Results of 2006 are from Moutou et al. (2007) and the results of 2007 and 2008 are from Fares et al. (2010). Error bars are calculated as in Mengel et al. (2016), i.e. by varying the input parameters within their error bars. The error bar on B orbit is of the order of 3 mG (see text for more details).

Figure 7 .

 7 Figure 7. Magnetic field evolution of HD189733. From left to right: the square-root of the total magnetic energy, the toroidal energy relative to the total one, and the energy in the axisymmetric modes of the poloidal component relative to the energy of the poloidal component. Error bars are calculated as stated in the text.

Table 1 .

 1 Summary of the physical properties of HD189733 and its hot Jupiter, HD189733b.

	Physical property	value	reference
	Star:
	Sp. type	K2V
	V (mag)	7.7
	T eff (K)	5050 ± 50	Bouchy et al. (2005)
	M (M )	0.92 ± 0.03	Bouchy et al. (2005)
	R (R )	0.76 ± 0.01	Winn et al. (2007)
	v sin i (km s -1 )	2.97 ± 0.22	Winn et al. (2006)
	Prot (d)	11.94 ± 0.16	Fares et al. (2010)
	i ( • )	∼ 85	Fares et al. (2010)
	dΩ (rad s -1 )	0.146 ± 0.049	Fares et al. (2010)
	Planet:
	i orb ( • )	85.76 ± 0.29	Boisse et al. (2009)
	Mp(M )	1.13 ± 0.03	Boisse et al. (2009)
	Rp(R )	1.154 ± 0.032	Boisse et al. (2009)
	P orb (d)	2.2185733 ± 0.0000019 Boisse et al. (2009)
	a (au)	0.031 ± 0.001	Boisse et al. (2009)
	i spin-orbit	0.85 +0.32 -0.28

Table 3 .

 3 Same as Table2for the observations in 2014 and 2015. The exposure time is 4 × 900 s, expect for 28 May 2015 with an exposure time of 4 × 800 s, and 20 July 2015 where only two sub-exposures of 900 s each were taken. All spectra were obtained using NARVAL. Dates marked with † were not used for the mapping of the stellar field.

	)	HJD	UT	S/N	Rot. Cycle Orb. Cycle	v rad	σ LSD
	(2013)	(2,456,000+)	(h:m:s)		(239+)	(1297+)	(km s -1 ) (10 -4 Ic)
	12 May †	424.58158	01:55:55	630	-6.0673	-37.0955	-2.134	0.50
	12 May †	424.62595	02:59:48	670	-6.0636	-37.0755	-2.154	0.48
	13 May †	425.59553	02:15:52	680	-5.9828	-36.6385	-2.340	0.46
	13 May †	425.63991	03:19:46	680	-5.9791	-36.6185	-2.325	0.46
	14 May †	426.62177	02:53:31	660	-5.8973	-36.1759	-2.029	0.47
	13 Jun	456.53057	00:38:37	440	-3.4049	-22.6948	-2.385	0.81
	16 Jun	459.52985	00:37:17	540	-3.1549	-21.3429	-2.045	0.62
	02 Jul	475.60489	02:24:02	650	-1.8153	-14.0973	-2.075	0.51
	05 Jul	478.60456	02:23:22	650	-1.5654	-12.7452	-2.477	0.47
	06 Jul	479.59405	02:08:10	580	-1.4829	-12.2992	-2.045	0.56
	08 Jul	482.52265	24:25:11	690	-1.2389	-10.9792	-2.218	0.47
	15 Jul	488.54927	01:03:14	380	-0.7366	-8.2627	-1.906	0.92
	02 Aug	507.50215	23:55:01	400	0.8428	0.2801	-2.298	0.87
	04 Aug	509.37413	20:50:4	660	0.9988	1.1239	-2.277	0.50
	05 Aug	510.51950	24:20:01	440	1.0942	1.6401	-2.00	0.71
	08 Aug	513.50839	24:04:03	530	1.3433	2.9873	-2.124	0.63
	09 Aug	514.50144	23:54:04	680	1.426	3.4349	-2.241	0.46
	10 Aug	515.50391	23:57:39	650	1.5096	3.8868	-2.016	0.49
	11 Aug	516.52369	24:26:09	660	1.5946	4.3464	-2.294	0.49
	14 Aug	518.53989	0:49:32	680	1.7626	5.2552	-2.326	0.46
	15 Aug	520.38524	21:06:53	530	1.9164	6.0870	-2.268	0.59
	18 Aug	523.37830	20:56:59	630	2.1658	7.4361	-2.225	0.51
	19 Aug	524.46529	23: 2:18	670	2.2564	7.926	-2.058	0.47
	21 Aug	526.49842	23:50:05	610	2.4258	8.8425	-2.027	0.53
	22 Aug	527.40235	21:31:48	650	2.5011	9.2499	-2.37	0.50
	23 Aug	528.50256	23:56:09	640	2.5928	9.7458	-1.938	0.51
	02 Sep	538.34180	20:05:15	680	3.4127	14.1807	-2.352	0.47
	03 Sep	539.40804	21:40:42	660	3.5016	14.6613	-1.985	0.47
	04 Sep	540.37393	20:51:39	650	3.5821	15.0967	-2.287	0.51
	08 Sep	544.47392	23:15:58	610	3.9237	16.9447	-2.125	0.53
	10 Sep	546.42285	22: 2:36	670	4.0862	17.8232	-2.014	0.48
	11 Sep	547.39768	21:26:26	660	4.1674	18.2626	-2.387	0.48
	12 Sep	548.39552	21:23:24	660	4.2505	18.7123	-2.082	0.30
	13 Sep*	548.86339	08:37:12 1020	4.2895	18.9232	-2.082	0.30
	13 Sep	549.40705	21:40:06	560	4.3348	19.1683	-2.346	0.60
	15 Sep	551.36003	20:32:34	420	4.4976	20.0486	-2.211	0.76
	17 Sep*	552.89635	09:25:02	900	4.6256	20.7410	-1.986	0.37
	19 Sep	555.38839	21:13:48	670	4.8333	21.8643	-2.242	0.34
	20 Sep*	555.73218	05:28:55	960	4.8619	22.0193	-2.242	0.34
	20 Sep	556.39231	21:19:33	640	4.9169	22.3168	-2.39	0.49
	21 Sep	557.35420	20:24:46	560	4.9971	22.7504	-2.183	0.34
	22 Sep*	557.88116	09:03:4	990	5.0410	22.9879	-2.182	0.34
	24 Sep* †	559.89429	09:22:47	870	5.2088	23.8953	-2.017	0.39
	24 Sep	560.32644	19:45:07	670	5.2448	24.0901	-2.36	0.33
	25 Sep*	560.80689	07:17:02	970	5.2848	24.3066	-2.36	0.33
	27 Sep*	562.89300	09:21:16	880	5.4587	25.2469	-2.379	0.37
	08 Oct †	574.28920	18:53:07	690	6.4084	30.3837	-2.332	0.46
	09 Oct †	575.28808	18:51:38	670	6.4916	30.8339	-2.007	0.47
	12 Oct †	578.28726	18:50:50	610	6.7415	32.1857	-2.329	0.53
	13 Oct †	579.29042	18:55:31	530	6.8251	32.6379	-2.008	0.63
	18 Oct †	584.28323	18:45:50	560	7.2412	34.8884	-2.037	0.60
							c 0000 RAS, MNRAS 000, 000-000

c 0000 RAS, MNRAS 000, 000-000

ACKNOWLEDGMENTS

The authors thank an anonymous referee for their useful comments. This work is based on observations obtained with ESPaDOnS at the Canada-France-Hawaii Telescope (CFHT) and with NARVAL at the Télescope Bernard Lyot (TBL). CFHT/ESPaDOnS are operated by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (INSU/CNRS) of France, and the University of Hawaii, while TBL/NARVAL are operated by INSU/CNRS. We thank the CFHT and TBL staff for their help during the observations, and in particular R. Cabanac and P. Petit. We also thank J.-F. Donati and E. Hébrard for useful comments on the data analysis. RF acknowledges financial support by WOW from INAF through the Progetti Premiali funding scheme of the Italian Ministry of Education, University, and Research. VB and AL acknowledge the support of the French Agence Nationale de la Recherche (ANR), under program ANR-12-BS05-0012 'Exo-Atmos'. Part of VB work has been carried out in the frame of the National Centre for Competence in Research "PlanetS" supported by the Swiss National Science Foundation (SNSF). V.B. also acknowledges the financial support of the SNSF. AAV acknowledges partial support from an Ambizione Fellowship of the Swiss National Science Foundation. ChH highlights financial support of the European Community under the FP7 by an ERC starting grant number 257431. P.W. is supported by a STFC consolidated grant (ST/L000733/7)