Elias Egho

Dino Ienco

Nicolas Jay

Amedeo Napoli

Pascal Poncelet

Catherine Quantin

Chedy Raïssi

Maguelonne Teisseire

Oi E N C O

Healtcare Trajectory Mining by Combining Multi-dimensional Component and Itemsets

Keywords: Sequential Patterns, Multi-dimensional Sequential Patterns, Data Mining

Introduction

Data warehouses are constituting a large source of informations that are used and exploited to extract useful knowledge for expert analysis and decision makers [START_REF] Cohen | Online expansion of largescale data warehouses[END_REF]. In temporal data warehouses, every bit of information is associated with a timeline describing a total order over events. This particular ordering introduces more complexity to the extraction process and more precisely to mining processes that enumerate patterns that encompass interesting transient events. Many efficient approaches were developed to mine these patterns (i.e., sequential patterns) like PrefixSpan [START_REF] Pei | Prefixspan: Mining sequential patterns by prefixprojected growth[END_REF], SPADE [START_REF] Zaki | Spade: An efficient algorithm for mining frequent sequences[END_REF], SPAM [START_REF] Ayres | Sequential pattern mining using a bitmap representation[END_REF], PSP [START_REF] Masseglia | The psp approach for mining sequential patterns[END_REF], DISC [START_REF] Chiu | An efficient algorithm for mining frequent sequences by a new strategy without support counting[END_REF], PAID [START_REF] Yang | Paid: Mining sequential patterns by passed item deduction in large databases[END_REF], FAST [START_REF] Salvemini | Fast sequence mining based on sparse id-lists[END_REF]. However, all these techniques and algorithms, without any exception, focus solely on sequences of set valued data (i.e., itemsets) and contrast with real-world data that have multiple dimensions. To overcome this problem, Pinto et al. [START_REF] Pinto | Multi-dimensional sequential pattern mining[END_REF] introduced the notion of multi-dimensionality in sequences and proposed an efficient algorithm. Later works, like Zhang et al. [START_REF] Zhang | Approxmgmsp: A scalable method of mining approximate multidimensional sequential patterns on distributed system[END_REF] or Yu et al. [START_REF] Yu | Mining sequential patterns from multidimensional sequence data[END_REF] extended the initial Pinto's approach for different scenarios and use-cases. While in set valued approaches the events are represented by itemsets, in multidimensional temporal databases the events are defined over a fixed schema in which all the attributes are mandatory in the extracted patterns. Furthermore, and this is particularly true in the data warehouse environment, background knowledge is usually available and can be represented as a hierarchy over the values of the attributes. Following this logic, Plantevit et al. introduced M3SP [START_REF] Plantevit | Mining multidimensional and multilevel sequential patterns[END_REF], an efficient algorithm that is able to incorporate different dimensions and their ordering (organization) in the sequential pattern mining process. The benefit of this approach is to extract patterns with the most appropriate level of granularity. Still, this idyllic representation of uniform data is very uncommon in real-world applications where heterogeneity is usually elevated to a foundational concept. In this study, we focus on extracting knowledge from medical data warehouse representing information about patients in different hospitals. The successive hospitalizations of a patient can be expressed as a sequence of multidimensional attributes associated with a set of medical procedures. Our goal is to be able to extract patterns that express patients stays along with combinations of procedures over time. This type of pattern is very useful to healthcare professionals to better understand the global behavior of patients over time. Unfortunately this full richness and complexity of the data cannot be exploited by any of the traditional sequential pattern mining techniques. In this paper, we propose a new approach to extract patterns from sequences which include multi-dimensional and set valued data at the same time. In addition, the proposed approach incorporates background knowledge in the form of hierarchies over attributes.

The remainder of this paper is organized as follows, Section 2 introduces the problem statement as well as a running example. The method for extracting multi-dimensional-itemset frequent patterns is described in Section 3. Section 4 presents experimental results from both quantitative and qualitative point of views and Section 5 concludes the paper.

Problem Statement

In this section we list some preliminary definitions needed to formalize and present our problem. First of all, we introduce a motivating example from a real data set related to the PMSI (Program of medical information systems). This French nationwide information system describes hospital activities from both economical and medical points of view. In this system, each hospitalization is related to the recording of administrative, demographical and medical data. A multi-dimensional-itemset data sequence is composed of events. Definition 4. (Multi-dimensional-itemset Sequence) A multi-dimensional-itemset sequence s =<e 1 ,e 2 ,...,e l > is an ordered list of events e i . Given two Multidimensional-itemset Sequences s =<e 1 ,e 2 ,...,e l > and s =<e 1 ,e 2 ,...,e l >, s is included in s , denoted by s ⊆ s s , if there exist indices 1 ≤ i 1 <i 2 <. . .< i l ≤ l such that e j ⊆ e e ij for all j =1...l and l l . Example 4. The sequence s = ((UH Paris ,Cancer), {p 1 ,p 2 }), ((GH Lyon ,R 1 ), {p 2 }) is a sequence of two events. It expresses the fact that a patient was admitted to the University Hospital of Paris UH Paris for a cancer disease Cancer and underwent procedures p 1 and p 2 , then he went to the General Hospital of Lyon GH Lyon for pneumonitis R 1 and underwent procedure p 2 .As e q u e n c es =

((UH Paris ,C 1 ), {p 1 }) is included in s, s ⊆ s s, because ((UH Paris ,C 1 ), {p 1 }) ⊆ e ((UH Paris ,Cancer), {p 1 ,p 2 }).

Definition 5. (Patient Trajectory) A patient trajectory is defined as am u l t idimensional-itemset sequence.

Example 5. In Table 1, the sequence s = ((UH Paris ,C 1 ), {p 1 ,p 2 }), ((UH Paris ,C 1 ), {p 1 }), ((GH Lyon ,R 1 ), {p 2 }) represents the trajectory for the patient P 1 . Let supp(s) be the number of sequences that includes s in S DB . Furthermore σ be a minimum support threshold specified by the end-user. Definition 6. (Most Specific Frequent Sequence) Let s be multi-dimensionalitemset sequences, we can say that, s is the most specific frequent sequences in S DB ,i fa n do n li f :supp(s) ≥ σ and s ∈ S DB , where supp(s)=supp(s ) and s ⊆ s s.

Problem 1. The problem of mining multi-dimensional-itemset sequences is to extract the set of all most specific frequent sequences in S DB such as supp(s) ≥ σ. By using the taxonomies we can extract more or less general or specific patterns and overcome problems of excessive granularity and low support. Example 6. Let σ =0.75 (i.e. a sequence is frequent if it appears at least three times in S DB ). The sequence

s 1 = (UH Paris ,C 1 ), {p 1 ,p 2 }), ((GH Lyon ,R 1 ), {p 2 } is frequent. s 2 = (UH,Cancer), {p 1 ,p 2 }), ((GH, Respiratory), {p 2 }} is also frequent. Nevertheless, s 2 is not kept since it is too general compared to s 1 .
In this section, we present the MMISP (Mining Multi-dimensional-Itemset Sequential Patterns) algorithm for extracting multi-dimensional-Itemset sequential patterns with different levels of granularity over each dimension. MMISP follows a bottom-up approach by first focusing on extracting mdc that can exist at different level of granularity, then it considers the itemset part of the events and compute the support of every item is S DB . After these two steps, frequent multidimensional components and frequent items are combined to generate events. In the final step, the frequent events are mapped to a new representation and a standard sequential mining algorithm is applied to enumerate multi-dimensional itemset sequential patterns.

In the next subections, we provide the details of each step of our work and discuss the different challenges.

Generating Multi-dimensional Components

MMISP starts by processing the multi-dimensional components of the sequences. Basically it considers three types of dimensions: a temporal dimension D t ,as e t of analysis dimensions D A and a set of reference dimensions D R .M M I S Ps p l i t s S DB into blocks according to reference dimension D R . Then, MMISP sorts each block according to the temporal dimension D t . The tuples of multi-dimensional component appearing in an event are defined w.r.t. analysis dimensions D A . The support of the multi-dimensional component is computed according to dimensions of D R . It is the ratio of the number of blocks supporting the multidimensional component over the total number of blocks. This a classic way of partitioning the database and was introduced in [START_REF] Plantevit | Mining multidimensional and multilevel sequential patterns[END_REF]. Following this partitioning step, MMISP generates all the frequent multidimensional components. Firstly, we generate the most general multi-dimensional component, that is (T 1 ,...,T m ). In our running example, we have two dimensions (hospital and disease), so the most general multi-dimensional component is (T hospital ,T disease ). Then, our approach generates all multi-dimensional components of the form (T 1 ,. . . T i-1 , d i , T i+1 ,. . . ,T m )w h e r ed i ∈ down(T i ). We take only the frequent multi-dimensional component which has support greater than σ. In our running example and for σ = 75% (3 blocks from 4), we have four new frequent multidimensional components: (UH,T disease ), (GH, T disease ), (T hospital ,Respiratory) and (T hospital ,Cancer).

We continue the recursive generation of the new multidimensional components by using each previously generated frequent multidimensional component (a). This is done with a pivot method that identifies an integer z which is the position of the last dimension in a and is not top T . For example if a=(UH,T Disease ), z is the first dimension (hospital) because the value for the hospital dimension (UH) and the second dimension (disease) has the value T disease .

For each dimension d k in a,w h e r ek ∈ [z, m], we replace d k with one of its specialization from the set down(d k ). For example, if a=(UH,T Disease ), we have z =1 and we can generate four new mdc s : {(UH Paris ,T Disease ), (UH Nancy ,T Disease ), (UH, Respiratory), (UH,Cancer)}. The first and the second multidimensional components are generated by replacing UH by down(UH)={UH Paris ,UH Nancy }, the third and the forth multidimensional components are generated by replacing T Disease by down(T Disease )={Respiratory, Cancer}.

We select only the frequent multidimensional components. For our previously example with σ = 75%, {(UH Paris ,T Disease ), (UH,Cancer)} are the new frequent multidimensional components generated by (UH,T Disease ).

Finally, from all frequent multi-dimensional components generated, we select only the most specific multi-dimensional component.

Definition 7. (Most Specific Multi-dimensional Component) Let a be multidimensional component, we can say that, a is the most specific multi-dimensional component, if and only if a multi-dimensional component

Generating Frequent Itemset

In this step, MMISP focuses on the itemset part of the sequences. Basically, this step aims at extracting the set of all items that are frequent in a sequence of length 1. Let us remind that usually, in level-wise approaches, either itemsetextension or sequence-extension can be considered. For example, if we have a sequence s 1 = {1, 2, 3},t h e ns 2 = {1, 2, 3}{4} is a sequence-extended sequence of s 1 and s 3 = {1, 2, 3, 4} is an itemset-extended sequence of s 1 . In our context focusing on sequence of length 1 we will only consider itemset-extension. Such an operation can be easily done by using any standard sequential pattern algorithm.

Patients Sequences of procedures

P1 {p1,p2}{p1}{p2} P2 {p1}{p1,p2}{p2} P3 {p1,p2}{p2} P4 {p2}{p3}{p2}

Sequences of procedures

Frequent Itemset Candidates {p1} {p2} {p1,p2}

The frequent itemset Fig. 4. The frequent itemset generated Example 9. Figure 4 shows the sequence of medical procedures for patients and the frequent itemset candidates for σ =0.75.

Generating Frequent Events

Generating frequent events is achieved by combining frequent multi-dimensional components with frequent itemsets. This task is achieved by building a prefix tree such that the first level in this tree is composed of the frequent multidimensional components and the second level is composed the frequent itemsets. More precisely, each branch in the tree represents an event. Then a scan is performed over the database to prune irrelevant events from the tree. For example, Figure 5 illustrates the tree before and after pruning infrequent events for σ =0.75. Combination tree before pruning Combination tree after pruning Fig. 5. An example of the tree for generating frequent events before and after the pruning

Extracting Frequent Multi-dimensional Itemset Patterns

Frequent sequences can then be mined by using any standard sequential pattern mining algorithm. As these algorithms require that the dataset to be mined is composed of pairs in the form (id, seq), where id is a sequence identifier and seq is a sequence of itemsets, we transform the initial dataset as follows:

-Each branch in the combination tree after pruning is assigned a unique id which will be used during the mining operation. This is illustrated in Table 3.

-Each block (patient) is assigned a unique id of the form P i .

-Every block b is transformed into a pair (P i , S(p i )), where S(P i ) is built according to the date and the content of the blocks. The final result is reported in Table 4.

A standard sequence mining algorithm can be applied on the transformed database.

event-id Frequent Event e1 (UH Paris ,C1), {p1} e2 (UH Paris ,C1), {p 2} e3 (UH Paris ,C1), {p1,p2} e4 (GH Lyon ,R1), {p 2} Table 3. Identification each branch (Event) in T id Sequence data P1 {e1,e2,e3}{e1}{e4} P2 {e1}{e1,e2,e3}{e4} P3 {e1,e2,e3}{e4} P4 {e2} Table 4. Transformed database
Then, the extraction of frequent sequences can be carried out. With σ =0.75, the pattern {e 3 }{e 2 } is frequent. This sequence transforms to ((UH Paris ,C 1 ), {p 1 ,p 2 }), ((GH Lyon ,R 1 ), {p 2 }) by using the identification in Table 3.

4E x p e r i m e n t s

We conduct experiments on both real and synthetic datasets. The algorithms is implemented in Java and the experiments are carried out on a MacBook Pro with a 2.5GHz Intel Core i5, 4GB of RAM Memory running OS X 10.6.8. The extraction of sequential patterns is based on the public implementation of CloSpan algorithm [START_REF] Yan | Clospan: Mining closed sequential patterns in large datasets[END_REF]. We use the implementation supplied by the IlliMine5 toolkit.

In order to assess the effectiveness of our approach, we run several experiments on the PMSI dataset. This database includes the following informations: Patients (id, gender . . . ), Stays (id, hospital, principal diagnosis, . . . ) and Medical Procedures (id, date,. . . ). Our dataset contains 486 patients suffering from lung cancer and living in the eastern region of France. The average length of data sequences is 27. The data is encoded using controlled vocabularies. In particular, diagnoses are encoded with the International Classification of Diseases (ICD10) 6 . This classification is used as an input taxonomy for MMISP. The ICD10 can be seen as a tree with two levels. As illustrated in Figure 6, 3-characters codes such as C34 (Lung cancer) have specializations: C340 is cancer of the main bronchus, C341 is cancer of upper lobe etc. Table 5. Care trajectories of 4 patients Table 5 is an example of care trajectories described over two dimensions (diagnosis, hospital ID) coupled with a set of medical procedures. For example (C341, 750712184), {ZBQK002} represents the stay of a patient in the Univer-sity Hospital of Dijon (coded as 210780581) treated for a lung cancer (C341), where the patient underwent chest radiography (coded as ZBQK002).

The experiments are designed to extract some multi-dimensional sequential patterns for helping the medical experts to describe some healthcare patients trajectories. For this experiment the support value is set to 15 (i.e. σ =0 .03). MMISP generates 121 different patients trajectories. Table 6 shows some patients trajectories obtained by our approach. Pattern 2 can be interpreted as follows: 42% of patients have a hospitalization in the University Hospital of Dijon for a lung cancer (210780581,C341), where they underwent supplement procedures (coded as YYYY030) for passing the chest radiography (coded as ZBQK002). Then, the same patients go to any university hospital for doing chemotherapy (CHU/CHR,Z511), where they underwent only the chest radiography (coded as ZBQK002). In the second experiments, we study the scalability of our approach. We consider the number of extracted patterns and the running time with respect to two different parameters, the number of the dimension in the multidimensional components and the average length of the itemsets in the data sequences. The first batch of synthetic data generated contains 10 000 sequences defined over (2, 3, 4 and 5) analysis' dimensions. Each sequence contains 30 events and each event is described, in average, by 15 items in the itemset. Hierarchical relations are defined over 5 levels of granularity between elements of each analysis dimension. Figure 7 reports the results according to different values of support threshold for different number of analysis dimension in multidimensional component. We can notice that the running time increases for each newly added analysis dimension. The second batch of synthetic data generated contains 10 000 sequences with varying number of items 5, 10, 15 and 20. The sequences in the four generated data sets have an average cardinality of 30 events,by 3 analysis dimensions. Hierarchical relations are defined over 5 levels of granularity between elements of each analysis dimension. Figure 8 reports the results according to different values of support threshold for different lengths of itemsets.

Conclusion

In this paper, we propose a new approach to mine multi-dimensional itemset sequential patterns. Our approach is able to capture knowledge from dataset represented over both multi-dimensional component and itemsets. We provide formal definitions and propose a new algorithm MMISP to mine this new kind of pattern. We conduct experiments on both real and synthetic datasets. The method was applied on real-world data where the problem was to mine healthcare patients trajectories. According to medical experts, new patterns are easier to understand.

  Example 3. e =((UH,T disease ), {p 1 ,p 2 ,p 3 }) is an event, where (UH,T disease )is a multidimensional component with two dimensions representing hospital and diagnosis. {p 1 ,p 2 ,p 3 } denotes the set of medical procedures. An event e = ((UH Paris ,C 1 ), {p 1 ,p 2 })i si n c l u d e di ne, e ⊆ e e, because (UH Paris ,C 1 ) mdc (UH,T disease ) and {p 1 ,p 2 }⊆{p 1 ,p 2 ,p 3 }.

Fig. 2 . 7 .

 27 Fig. 2. Block partition of the database according to DR={Patient}

  , where supp(a)= supp(a ) and a mdc a.

Frequent

  multi-dimensional comp onent (UH Paris ,C1) (GH Lyon ,R1)

!Fig. 3 .

 3 Fig. 3. Frequent multi-dimensional components generation

Fig. 6 .

 6 Fig. 6. Examples of taxonomies used in multilevel sequential pattern mining

Fig. 7 .

 7 Fig. 7. Running Time (left) and Number of extracted pattern (right) obtained by MMISP with varying in the number of dimension

Fig. 8 .

 8 Fig. 8. Running Time (left) and Number of extracted pattern (right) obtained by MMISP with varying itemsets' cardinalities

Table 1 .

 1 An example of a database of patient trajectoriesLet S DB be a database of multi-dimensional-itemset data sequences. Table1illustrates such a database. Figure1shows two hierarchical taxonomies characterizing both the hospital and diagnosis dimensions. For hospital dimension, UH Paris ∈ down(UH) as UH Paris is a direct descendant of UH.

	Patients					Tra jectories
	P1	((UH Paris ,C1), {p1,p2}), ((UH Paris ,C1), {p1}), ((GH Lyon ,R1), {p2})
	P2	((UH Paris ,C1), {p1}), ((UH Paris ,C1), {p1,p2}), ((GH Lyon ,R1), {p2})
	P3		((UH Paris ,C1), {p1,p2}), ((GH Lyon ,R1), {p2})
	P4	((UH Paris ,C1), {p2}), ((UH Paris ,R2), {p3}), ((GH Lyon ,R2), {p2})
	Definition 1. (Dimensions and specialization down(d)) A dimension (D, ) is
	a partially ordered set. For a given d in D, down(d) (resp. up(d)) denotes the
	set of all specializations {x ∈ D|x d} (resp. generalizations {x ∈ D|d x})o f
	d.
	Example 1. !"#$%&'#()*+,'-#(./*0!+1* 2%"%&./*+,'-#(./*03+1* T >,'-#(./ *	8%'-#&.(,&) *	T C#'%.'% *	*3."6%&*
	!+ 4.&#' *	!+ 5."6) *	2+ 4.&#' *	2+ 7)," *	*4"%9:,"#;'*08 < 1 *	='(>:.*08 ? 1*	**79"A*3."6%&*03 < 1*	*@&%.'(*3."6%&*03 ? 1 *
		+,'-#(./'*						B#'%.'%'*
			Fig. 1. Hospital and diagnoses taxonomies
	By taking into account taxonomies, we define a multi-dimensional component
	as follows:
	Definition 2. (Multi-dimensional component) Given a dimension (D, ),amulti-
	dimensional component over D, denoted (mdc, mdc ),i sat u p l e(d 1 ,...,d m )
	where d i ∈ D, i =1, •••,m. For two given multidimensional components mdc =
	(d 1 ,...,d m ) and mdc =(d	1 ,...,d	m ), mdc mdc mdc denotes that mdc is more
	specific than mdc, if for every i =1,...,m, d	i ∈ down(d i ).
	Example 2. Let (UH Paris , Lung Cancer) and (UH, Cancer)b et w om u l t i d i -
	mensional components. (UH Paris , Lung Cancer) mdc (UH, Cancer) because
	UH Paris ∈ down(UH) and Lung Cancer ∈down( Cancer).

1 ,...,p n }.

Table 2 .

 2 The most specific frequent multi-dimensional components Example 8. Fiqure 3 illustrates the mechanism work of generation all frequent multi-dimensional components on our running example with σ =0 .75. We can notice that the most specific components are (UH Paris ,C 1 ) and (GH Lyon ,R 1 ).

Table 6 .

 6 Some healthcare patients trajectories obtained by MMISP

			Runtime sequences over frequency threshold		#Frequent sequences over Frequency threshold
		250				2-dim 3-dim		4000				2-dim 3-dim
						4-dim						4-dim
						5-dim						5-dim
		200						3500
	Run Time(s)	150					Number of Patterns	2500 3000
		100
								2000
		50
		1000	1200	1400	1600	1800	2000	1000	1200	1400	1600	1800	2000
				Frequency					Frequency

http://illimine.cs.uiuc.edu/

http://apps.who.int/classifications/apps/icd/icd10online/


Nhat Phan

Dino Hai

Pascal Ienco

Maguelonne Poncelet
email: maguelonne.teisseire@teledetection.frpascal.poncelet@lirmm.fr

Teisseire

GeT Move: An Efficient and Unifying Spatio-Temporal Pattern Mining System for Moving Objects

Keywords: Spatio-temporal patterns, trajectories, visualization

Recent improvements in positioning technology has led to a much wider availability of massive moving object data. A crucial task is to find the moving objects that travel together. Usually, these object sets are called spatio-temporal patterns. Analyzing such data has been applied in many real world applications, e.g., in ecological study, vehicle control, mobile communication management, etc. However, few data mining tools are available for flexible and scalable analysis of massive scale moving object data. The main reason is that there is no a unifying approach to manage the patterns while many different kinds of spatio-temporal patterns have been proposed in recent years. Each approach only focuses on mining a specific kind of pattern. Our system, GeT Move, is designed to extract and manage different spatio-temporal patterns concurrently. A user-friendly interface is provided to facilitate interactive exploration of mining results. Since GeT Move is tested on many kinds of real data sets, it will benefit users to carry out versatile analysis on these kinds of data by exhibiting different kinds of patterns at the same time.

Introduction

Nowadays, many electronic devices are used for real world applications. Telemetry attached on wildlife, GPS installed in cars, sensor networks, and mobile phones have enabled the tracking of almost any kind of data and has led to an increasingly large amount of data that contain moving objects and numerical data. Therefore, analysis on such data to find interesting patterns is attracting increasing attention for applications such as movement pattern analysis, animal behavior study, route planning and vehicle control.

Despite the growing demands for diverse applications, there have been few scalable tools available for mining massive and sophisticated moving object data. Even if some tools are available for extracting different kinds of patterns, they suffer the problem that they only extract a specific kind of pattern at a moment.

When considering a dataset, it is quite difficult, for the decision maker, to know in advance the kind of pattern embedded in the data. Additionally, it is difficult to recognize the differences between patterns. Therefore proposing an system able to automatically extract all these different kinds of patterns can be very useful. To tackle the issue, our system, GeT Move is designed to reveal collective movement patterns like convoys [START_REF] Jeung | Discovery of Convoys in Trajectory Databases[END_REF], group patterns [START_REF] Wang | Efficient Mining of Group Patterns from User Movement Data[END_REF], closed swarms [START_REF] Li | Swarm: Mining Relaxed Temporal Moving Object Clusters[END_REF], moving clusters [START_REF] Kalnis | On Discovering Moving Clusters in Spatiotemporal Data[END_REF] and also periodic patterns [START_REF] Mamoulis | Mining, Indexing, and Querying Historical Spatiotemporal Data[END_REF]. GeT Move has been applied on many different datasets published by the Movebank 1 project.

The GeT Move System Architecture

GeT Move general architecture 2 , described in Figure 1-a, has three main layers: (i) collection and cleaning, (ii) mining, and (iii) visualization. The bottom layer is responsible for collecting and preprocessing of moving objects. Object movements are quite complex and therefore data can be noisy. So preprocessing is needed to integrate and clean the raw data and to interpolate missing points.

Unifying approach is then applied on the preprocessed data on the middle layer. Then GeT Move can automatically extract different kinds of patterns such as convoys, closed swarms, group patterns, moving clusters and periodic patterns. The brief definitions of such patterns and an overview of GeT Move are presented in the following.

The top layer shows the visualized results with some statistics. They can be plotted on 2D plane or embedded into other visualization tools (e.g., Google Map 3 and Google Earth 4 ). Our system enable users to visualize all the patterns with some statistics. Moreover, different patterns can be visualized in the same time on the same map. Consequently, users can be more insights into these results.

Informally, a swarm is a group of moving objects O containing at least ε individuals which are closed each other for at least min t timestamps T . To avoid redundant swarms, Zhenhui Li et al. [START_REF] Li | Swarm: Mining Relaxed Temporal Moving Object Clusters[END_REF] propose the notion of closed swarm for grouping together both objects and time. A swarm (O, T ) is a closed swarm if it cannot be enlarged in terms of timestamps T and objects O. Another pattern is convoy which is also a group of objects such that these objects are closed each other during at least min t time points. The main difference between convoy and swarm (or closed swarm) is that convoy lifetimes must be consecutive. Furthermore, moving clusters can be seen as special cases of convoys with the additional condition that they need to share some objects between two consecutive timestamps [START_REF] Han | Mining Moving Object, Trajectory and Traffic Data[END_REF]. We can consider that the main difference between convoys and swarms is about the consecutiveness and non-consecutiveness of clusters during a time interval. In [START_REF] Wang | Efficient Mining of Group Patterns from User Movement Data[END_REF] group pattern is a set of disjointed convoys which are generated by the same group of objects in different time intervals. By considering a convoy as a time point, a group pattern can be seen as a closed swarm of disjointed convoys.

By taking into account the similar characteristics of the patterns, our system is able to extract all the patterns. Additionally, GeT Move provides a platform for users to flexibly tune parameters and supports visualization of the results in different formats. The output can be written in Google Map and Google Earth format to help users better explore the results. Furthermore, users can compare the results in different datasets to understand the differences between object movement behaviors. Moreover, thanks to the system, one can explore that which kinds of patterns are embedded in each dataset. Moreover, users can plot as many patterns as users would like at the same time or clean and redraw the map with the other patterns. Additionally, to due with the large scale moving object databases, GeT Move is implemented by combining different techniques such as clustering technique and pattern mining technique so that GeT Move is very efficient. Therefore, the system can performed online and therefore it is easier for interested one to use and discover interesting patterns.

Let us examine a real example of Swainsoni and Buffalo movements. The raw trajectories of 43 Swainsonies are gathered from July 1995 to June 1998. Buffalo dataset concerns 165 buffalos and the tracking time is from year 2000 to year 2006. We further select the parameters ε and min t and then the system automatically discover different patterns. See Figure 2, there are two closed swarms discovered from Swainsoni and Buffalo. Additionally, with the same parameters there are closed swarms, convoys and group patterns in Swainsoni while there are only closed swarms in Buffalo. Therefore, it is interesting for users to recognize the differences between Swainsonies and Buffalos.

The GeT Move effectiveness and efficiency have been evaluated by using real and synthetic datasets. All the experimental results are also integrated on our demo website.

Outlook

In this paper, we propose a system, GeT Move, which is designed to automatically extract different kinds of spatio-temporal patterns by applying a unifying system. By experimenting with real data sets, one can efficiently observe interesting patterns, and also find the limitations of the current kinds of patterns. This is useful to promote further research into the new challenge issues in moving object mining.

Fig. 1 .

 1 Fig. 1. General System.

  4P h a n N h a t H a i e t a l . (a) One of discovered closed swarms in Swainsoni (b) One of discovered closed swarms in Buffalo

Fig. 2 .

 2 Fig. 2. Screen shots of GeT Move.

  , Hwang et al. propose a general pattern, called a group pattern, which essentially is a combination of both convoys and closed swarms. Basically,

	(a) General system
	1 http://www.movebank.org/
	2 http://www.lirmm.fr/∼phan/index.jsp
	3 http://code.google.com/apis/maps/
	4 http://earth.google.com/


Vincent Russier

Caroline De Montferrand

Yoann Lalatonne

Laurence Motte

Size and polydispersity effect on the magnetization of densely packed magnetic nanoparticles

The magnetic properties of densely packed magnetic nanoparticles (MNP) assemblies are investigated from Monte Carlo simulations. The case of iron oxide nanoparticles is considered as a typical example of MNP. The main focus is put on particle size and size polydispersity influences on the magnetization curve. The particles are modeled as uniformly magnetized spheres isolated one from each other by a non magnetic layer representing the organic coating. A comparison with recent experimental results on γ-Fe2O3 powder samples differing by their size is given.

I. INTRODUCTION

The physics and chemistry of nanoscale magnetic particles (MNP) still gives rise to an important research activity due both to their wide range of potential applications and their own fundamental interest [START_REF] Dormann | Advances in chemical physics[END_REF][START_REF] Skomski | [END_REF][3][4]. Among the large variety of MNP, iron oxide based ones γ-Fe 2 O 3 and Fe 3 O 4 take a particular place in the field of biological and medical applications because of their bio-compatibility and suitable superparamagnetic properties. To translate intrinsic properties of nanoparticles to various applications, there is a need to control nanoparticle dispersions. Consequently nanoparticles are usually coated by an organic surfactant [5,6] in order to prevent aggregation. The influence of this non magnetic layer and then the nanoparticles contact distance play a major role on collective magnetic properties [7]. A complete understanding of the macroscopic magnetic properties of MNP assemblies in terms of their individual intrinsic characteristics on the one hand and of the size distribution and volume concentration on the other hand is of crucial importance. Indeed this is a mean to get informations on the relevant parameters of the distribution and MNP properties from the magnetic measurements. Two key features which strongly influence the macroscopic magnetic properties of these systems are the magnetic structure at the particle scale, where core shell structure and spin canting effect can be invoked [8][9][10][11] and the size distribution generally described through a lognormal law for the diameters distribution.

At temperatures higher than the blocking temperature T b where the MNP are in the superparamagnetic regime [START_REF] Dormann | Advances in chemical physics[END_REF]4] and in case of weak interparticle interactions, namely for both small particles concentrations and in the absence of cluster formation the physical properties of MNP assemblies are well understood. The magnetization curve, M (H) of the whole assembly follows then a Langevin like function weighted by the diameter distribution function and eventually modified in order to take into account a core-shell structure [12][13][14]. Moreover the one-body magnetocrystaline anisotropy energy of the MNP can also be taken into account and this modifies the M (H) curve from the original Langevin function [START_REF] García-Palacios | Advances in chemical physics[END_REF][START_REF] Wiekhorst | [END_REF]. The core shell structure of the MNP may consist simply of the inclusion of a magnetic dead layer at the surface of the MNP [10,11] or of the introduction of an additional paramagnetic component in the MNP [14]. The symmetry breaking at the surface can lead to surface effects on the anisotropy energy of each MNP with noticeable effects on the M (H) curve [17][18][19] . In case of diluted assemblies of spherical MNP when the particles are non or weakly interacting, the non interacting particles type of approach of the magnetization curve leads to a reasonable determination of the characteristics of the individual particles and of the size distribution namely the median diameter d m and the ln(d) standard deviation σ. However, when the NP concentration increases, the interparticles interactions must be taken into account.

These ones which for spherical and well coated MNP include mainly the interactions between the ods are thus available going from mean field approximation, thermodynamic perturbation theory (TPT) [18,20] for weakly interacting systems to numerical simulations for moderate to strongly interacting systems [18,21,22]. The mean field and TPT provide an illustrative physical picture of the relation between the local structure and either the magnetization in terms of the applied field or the susceptibility. For instance the demagnetizing field effect depending on the external shape of the system, is well reproduced by the TPT [20]. As a link between TPT and numerical simulations, the description based on the interaction fields distributions [23] which explains the DDI induced reduction of the magnetization of an isotropic system as a generalization of a similar result obtained using the TPT and suggests that the DDI induced reduction of the magnetization is not related to an antiferromagnetic behavior. However, for strongly interacting systems, as in lyophilized powder samples or high concentration MNP assemblies embedded in a non magnetic matrix the numerical simulations seem more adapted. Although numerical simulations of magnetic properties of MNP assemblies are now many, a systematic study of the mean size and polydispersity effects especially for randomly organized particles with high concentration is still missing.

The aim of this work is to investigate this problem and to interpret recent experimental measurements [24] on powder samples of maghemite MNP assemblies differing by their median size.

We present a Monte Carlo (MC) simulation of the mean particle size and polydispersity effect on the DDI in random and densely packed spherical clusters of coated spherical maghemite MNP.

Our main purpose is to model the case of lyophilized powders or high concentration of particles embedded in a non magnetic matrix. A particular attention is paid to the linear susceptibility χ, and its dependence on the median size of the size distribution. It is found that χ as a function of d m may present a plateau, leading to a quasi independence of the magnetization with respect to d m in the vicinity of the low external fields. The magneto crystalline anisotropy is then shown to play a role for larger values of the field when the particles remain in the superparamagnetic regime in agreement with the findings of Ref. [START_REF] García-Palacios | Advances in chemical physics[END_REF]25] for non interacting particles, in the TPT regime [20] and in preceding MC simulations [18,22,26]. As an application, we focus on the experimental magnetization curves of Ref. [24].

II. MODEL FOR DENSELY PACKED ASSEMBLIES

The model we use is designed to simulate the properties of either lyophilized powders samples or high concentration nanoparticles assemblies embedded in non magnetic matrix. As is usually done to model single domain MNP, the nanoparticles are modeled as non overlapping spheres bearing at their center a permanent point dipole representing the uniform magnetization of the particle (super spin). The moment of each particle is equal to its volume times the bulk magnetization, M s , which means that neither spin canting effect nor magnetic dead layer at the particle surface is considered. We also include the magneto crystalline anisotropy with the same anisotropy constant K 1 on all particles. The particles are supposed to be coated by a non magnetic layer of thickness ∆/2, representing the usual coating by organic surfactant molecules. The layer thickness is taken as ∆/2 for convenience (see below). The particle diameters, {d i } are distributed according to a log-normal law defined by the median diameter d m and the standard deviation σ of ln(d),

f (d) = 1 d √ 2πσ exp - (ln(d/d m )) 2 2σ 2 (1)
d m and σ are related to the mean diameter and the diameter standard deviation σ d through

d 1 = d m e σ 2 /2 and σ d = d m (e σ 2 -1)e σ 2 .
In the following, we use d m as the unit of length, and thus in reduced unit, the distribution function is totally determined by the single parameter σ which characterizes the system polydispersity.

We consider mainly spherical clusters, where owing to the global shape isotropy the demagnetizing effects vanish, with free boundary conditions. This choice of large spherical clusters can be justified on the experimental point of view since upon drying the NP are likely to aggregate in spherical shaped large clusters which has been confirmed from simulations [5]. Our first purpose is to focus on the contribution of the dipolar interactions (DDI) to the magnetization curve, especially in the moderate to strong coupling regime when particles surrounded by their coating layer are at contact. The geometrical configuration of two particles of different sizes at contact with their coating layer is displayed in figure [START_REF] Dormann | Advances in chemical physics[END_REF]. Moreover, we consider temperatures such that the particles of size d m are superparamagnetic; as we shall see later for polydisperse systems due to the presence of large particles in the distribution, this condition may not be strictly fulfilled. When taken into account the magnetocrystalline anisotropy is considered in its simplest form, namely in the uniaxial symmetry and at lowest order [START_REF] Skomski | [END_REF]4]. The total energy thus includes the DDI, the one-body anisotropy term and the Zeeman term corresponding to the interaction with the external applied field H a = H a ĥ. Let { r i }, {v(i)} , { m i } and { n i } denote the particles locations, volumes, moments and easy axes respectively. The total energy of the cluster reads

E = µ 0 4π i<j m i m j mi mj -3( mi rij )( mj rij ) r 3 ij -K 1 i v(i)(n i mi ) 2 -µ 0 H a i m i mi ĥ (2)
where hated letters denote unit vectors, m i are the moment magnitudes, r ij = | r i -r j |. It is worth mentioning that the consideration of the anisotropy term with a fixed easy axes distribution means that the magnetization relax according to a Néel process [START_REF] Du | Trémolet de Lacheisserie[END_REF][START_REF] Coey | Magnetism and magnetic materials[END_REF], namely the particles are considered fixed while their moment relaxes relative to their easy axis. In this work only the case of a random distribution of easy axes is considered. In the following we use reduced quantities; first the energy is written in k B T 0 units, T 0 being a suitable temperature (T 0 = 300K in the present work) and we introduce a reference diameter, d ref .

The reference diameter, d ref is a length unit independent of the size distribution, useful for the energy couplings, and can be chosen from a convenient criterion independently of the actual structure of the MNP assembly. The reduced total energy is given by

β 0 E = -ǫ (0) K d m d ref 3 i d * 3 i (n i mi ) 2 -ǫ (0) d d m d ref 3 i<j d * 3 i d * 3 j mi mj -3( mi rij )( mj rij ) r * 3 ij -h i d * 3 i mi ĥ ǫ (0) d = β 0 µ 0 4π (π/6) 2 M 2 s d 3 ref ǫ (0) K = β 0 K 1 v(d ref ) h = β 0 µ 0 M s (π/6)d 3 m H a ≡ d m d ref 3 H a H ref (3)
where Concerning the structure in position, the nanoparticles surrounded by their coating layer ∆ form an assembly of hard spheres of effective diameters {d i + ∆} (see figure [START_REF] Dormann | Advances in chemical physics[END_REF]) which are arranged in large densely packed clusters with either a random or a well ordered structure (simple cubic or face centered cubic lattice). We build these clusters in two steps. First a large stacking of the coated spheres is made in a parallepipedic box with the desired structure, random or well ordered. In the random case, this first step is made from a sequential random rain plus compression algorithm in such a way to maximize the packing fraction. Doing this we can get a packing fraction ϕ for the effective spheres corresponding to the so-called loose random packing [START_REF] Soppe | [END_REF] (ϕ ≃ 0.60 in the monodisperse case). Once this first step is performed, we cut within the global stacking the cluster we want to study by imposing both the external shape, either spherical or prismatic, and the number of particles N p , with typically N p ≃ 1000. The central part of some of the clusters used in the present work corresponding to different values of the polydispersity, σ is shown in figure [START_REF] Skomski | [END_REF] to illustrate the structures obtained. It is important to note that because of the coating layer of thickness ∆/2 the closest distance of approach between particles i, j is shifted from (d i + d j )/2 to (d i + d j )/2 + ∆ and as a result the sum involved in equation (3) depends on the actual magnetic particles concentration of the cluster through the value chosen for ∆. One can rewrite the DDI sum by using another length scale, namely (d m + ∆) in order to exhibit a contribution independent of ∆. Doing this, the total DDI energy reads

β 0 = (k B T 0 ) -
ǫ (0) d d m d ref 3 d m d m + ∆ 3 i<j d * 3 i d * 3 j mi mj -3( mi rij )( mj rij ) (r ij /(d m + ∆)) 3 (4)
We recall that the distribution of reduced diameters, {d * i } depends only on the value of σ, which is conserved through a scaling operation corresponding to a change of d m . The sum of equation ( 4) is a geometric sum characteristic of the DDI expected, at least for small values of σ, to be independent of ∆ and thus to characterize the reduced DDI sum of the most concentrated cluster (∆ = 0) of the structure (s.c., f.c.c., random) considered. In other words, equation ( 4) allows to explicit the dependence of the dipolar coupling with respect to the particles volume fraction, ϕ v .

For this we note that (d m /(d m + ∆)) 3 can be rewritten as ϕ v /ϕ m where ϕ m (≡ ϕ v (∆ = 0)) is the maximum value of ϕ v for the given configuration, namely the volumic fraction corresponding to the spheres of diameters {d i + ∆}. ϕ m = π/6, 0.74, and ∼ 0.60 for the simple cubic, fcc and the loose random packed structures respectively. Then from (4), we can introduce an effective dipolar coupling constant, say ǫ

(ef f ) d = ǫ (0) d (d m /d ref ) 3 (d m /(d m + ∆)) 3 which is rewritten as ǫ (ef f ) d = (ϕ v /ϕ m )(d m /d ref ) 3 , since ǫ (0) d =1.
0. Now, one can replace both d m and ∆ by say d m2 and ∆ 2 respectively in such a way that the total DDI energy remains constant by imposing

d m d ref d m d m + ∆ = d m2 d ref d m2 d m2 + ∆ 2 (5)
leading to

d m2 = d m 2(1 + ∆/d m ) 1 + 1 + 4 ∆ 2 d m 1 + ∆ d m 1/2 (6)
In the absence of anisotropy energy, namely when only the DDI is taken into account, the two systems characterized by (d m , ∆) and (d m2 , ∆ 2 ) are similar and therefore present the same magnetization curve in terms of the reduced field h. Furthermore this holds also whatever the value of ǫ K in the vicinity of zero external field because for random distribution of easy axes the linear susceptibility χ does not depend on ǫ K in the superparamagnetic regime. Doing the transformation (6), the actual values of { r i } are scaled according to the value of (d m + ∆). Our hypothesis of a value of σ for the reduced diameter distribution to be not (or only negligibly) modified holds rigorously in the quasi monodisperse case (σ << 1). Consequently we shall use in the following the scaling transformation (6) only in quasi monodisperse situations.

In the present work we focus on both the reduced magnetization per unit magnetic volume in the direction of the external applied field,

M r = M (h) M s = < i v(i) mi ĥ > i v(i) (7)
and the linear susceptibility,

χ = ∂M ∂H = M s H ref d m d ref 3 ∂M r ∂h = M s H ref d m d ref 3 χ r with χ r = ∂M r (h) ∂h (8)
where we have used equation ( 3) to introduce the reduced susceptibility, χ r . The susceptibility can also be obtained from the fluctuations :

χ r = β * i v(i) v(d m )    < i v(i) mi ĥ 2 > ( i v(i)) 2 -M 2 r    . (9)
As a rule, we use this second way with the direct derivative merely used as a check of the calculation.

When the anisotropy energy is zero, the magnetization curve can be simulated either starting from h = 0 and increasing the field step by step or from the starurated situation, and decreasing h down to h = 0. When the anisotropy energy is included and since we may get an opening of the hysteresis loop, we start from the saturated case at sufficiently high applied field, and decrease the field beyond -h irr where the irreversible field h irr is defined as the value of h below which the hysteresis cycle opens. In cases where the hysteresis cycle opens, we also define an anhysteretic magnetization curve from the downward and the upward magnetization curves which because of the symmetry of our system reads

M (an) r (h) = 1 2 (M (d) r (h) + M (u) r (h)) = 1 2 (M (d) r (h) -M (d) r (-h)) (10)
The magnetization curves M (h) in terms of the external field are determined from Monte Carlo simulations, by fixing the locations of the particles in the cluster. We consider free boundary conditions, and the clusters includes ca 1000 particles. The dipolar coupling parameter is determined from equation (3). In section III C we consider a given set of experimental results in order to illustrate the model; nevertheless we do not restrict this latter only to this well specified set of samples but instead use the characteristics of maghemite as typical example for MNP assemblies. For the bulk magnetization M s we use a commonly accepted value for maghemite. Using M s varying from 80 to 84 emu/g, or ∼ 75.0 emu/g, if we take into account the temperature dependence, and ρ = 4.870g/cm 3 leads to µ 0 M s from 0.459 T to 0.514 T; at T 0 = 300K we get ǫ found to be much larger when the particle size decreases. A rather wide spectrum of values can be found in the literature for K 1 , corresponding to (K 1 /K b ) lying in between ∼ 4 to 15 for particle diameters of ca 12 nm or smaller [30][31][32][33][34]. In the following we use either (K 1 /K b ) ≃ 4 or 2, since we consider particles with mean diameters larger than 10 nm. With K b = 0.47 10 4 J/m 3 [10,31] this leads to ǫ (0) K = 2.38 and 1.19 respectively. In any case, both M s and K 1 are to be understood merely as realistic instead of truly accurate experimental values given the simplicity of the model.

Our Monte Carlo simulations are based on the usual Metropolis scheme [START_REF] Binder | Monte carlo simulation in statistical physics[END_REF][START_REF] Allen | Computer simulation of liquids[END_REF]; the averages are taken over 10 to 40 independent runs each of which consists in 10 4 to 2.10 4 MC steps (MCS) of equilibration followed by 2.10 4 to 3.10 4 MCS for the averages calculations. Each MCS consists in one trial move per moment in average. The trial move on the unit moment mi consists in moving mi to ( mi + w u)/ | mi + w u| where u is a random vector picked within the unit sphere with uniform probability density. This remains to move mi in a cone of maximum deviation δθ whose value is controlled by the amplitude parameter, w. For δθ << 1, we have δθ ≃ w. The value of w can be either fixed for a time scale mapping of the MCS or determined in a self consistent way in order to optimize the sampling by imposing a value for the acceptance ratio, R. The former version of this scheme corresponds to the time quantified Monte Carlo algorithm (TQMC) [START_REF] Nowak | [END_REF]38] in its first formulation ignoring the precessional step [START_REF] Nowak | [END_REF]. In the absence of anisotropy energy, the time scale mapping is irrelevant for the present purpose since we expect neither a ferromagnetic behavior nor a metastable blocked regime. Thus in this case, w is self consistently determined in such a way that R = 0.5. Conversely, when ǫ (0) K = 0, especially for polydisperse distributions we expect the largest particles to be in blocked state leading to a remanent state all the more that the DDI increase the blocking temperature. Hence, especially in the vicinity of h = 0, we deal with a metastable state whose life time must be comparable to the long scale measuring time τ m .

Strictly speaking one has to perform MC simulations corresponding to τ m and to use the version of the scheme outlined above allowing a mapping of the MC step on the true relaxing time. Since we are interested only in the long time behavior (corresponding to the SQUID measurements time scale), we do not focus on a precise mapping of the MCS scaling time. Instead, we determine w from the behavior of the instantaneous polarization M (t), versus t in MCS along a MC run at h = 0 starting from { mi } = ẑ. In other words, we chose w in order to avoid nonphysical jumps over the anisotropy energy barrier. By varying w we get as expected a w dependent evolution of M (t) before reaching a fluctuating behavior around a well defined plateau; the long time mean value < M (t) > determined beyond some threshold t value and for t up to 2.10 5 MCS is found independent of w at least for w varying in the range w = 0.03 to 0.25 for typical values of the parameters we consider (ǫ d ≃ 2 to 8, ǫ (0) K ≃ 2.3) and the polydispersity deviation σ = 0.28.

Therefore, in the following, we fix w = 0.25 when ǫ (0) K = 0.

III. RESULTS AND DISCUSSION

A. Weak coupling case

Before focusing on the powder like situation characterized by a moderate to strong dipolar coupling, we consider the weak coupling limit of the DDI, ǫ d < 1 with ǫ K = 0 where one can compare the results to the analytical one obtained from the thermodynamic perturbation theory and make the link with the mean field approximation. The important point is that one can deduce at least qualitatively when ǫ d deviates from the limit ǫ d << 1, the general behavior of the magnetization with respect to the DDI. In this framework, we can expand both the magnetization M r (h) and the susceptibility χ r in terms of ǫ d [18,20].

M r (h) = M (0) r (h) + M (1) r (h)ǫ d + 1 2 M (2) r ǫ 2 d χ r = χ (0) r + χ (1) r ǫ d + 1 2 χ (2) r ǫ 2 d ( 11
)
M (0) r and χ (0)
r correspond to the non interacting case, namely

M (0) r (h) = v(d)L(β * (d/d m ) 3 h)f (d)d(d) v(d)f (d)d(d) ( 12
)
where L is the Langevin function;

χ (0) r is directly related to M (0)
r (h) and at h = 0 leads to the linear susceptibility of the non interacting system

χ (0) r (0) = β * (d * 6 /d * 3 )/3 = β * exp(27σ 2 /2)/3 (13)
where d * s is the s-th reduced moment of the distribution f (d). Equation ( 13) explicits the effect of the polydispersity through the factor (d * 6 /d * 3 ), written here in terms of σ for the lognormal distribution. The expansions (11) which have been explicited in the framework of the TPT in [20,26] depend on geometrical sums which can be directly calculated from the structure considered. More-over, the linear part with respect to h of M

H ef f = H a -(D α -1/3)M (0) v (H a ) ; ( 14
)
D α is the demagnetizing factor of the sample in the direction of the external field, ĥ = e α , and M v is the total magnetization per unit volume which is related to either the number of MNP per unit volume, ρ, or the MNP volumic fraction, ϕ v , through

M (0) v = M s M (0) r (H a )ρ v(d)f (d) = M s M (0) r (H a )ϕ v (15)
Using equation (12) for M (0) r and keeping only the first order term with respect to h, we get

H ef f = H a - 1 3 (D α -1/3)M s β * d * 6 d * 3 h (16)
which is then inserted in the mean field expression for M r (H a ); then form an expansion of M (0) r at first order with respect to h and from equation (3) for ǫ d , we get

M (1) r (h) = -β * 2 (D α -1/3) 8 3 ϕ v d * 6 d * 3 2 h (17)
Equation ( 17) can be equivalently rewritten, in terms of ∆M r = M r (ǫ d ) -M r (ǫ d = 0), as

∂(∆M r (h)) ∂ǫ d = -β * 2 (D α -1/3) 8 3 ϕ v d * 6 d * 3 2 h and C 2 = ∂ 2 ∆M r (h) ∂ǫ d ∂h (h = 0, ǫ d = 0) ≡ ∂χ r ∂ǫ d = -β * 2 (D α -1/3) 8 3 ϕ v d * 6 d * 3 2 (18)
A result in agreement with refs. [20] and [26] in the monodisperse case. Here, the important point is that we explicitly write down the effect of the polydispersity through the factor (d * 6 /d * 3 ) 2 which strongly deviates from unity once σ takes non negligible values. It is worth mentioning that the preceding equations hold when either ǫ d << 1 or ϕ v << 1. We have performed MC simulations of the magnetization at small values of the coupling constant for prismatic clusters corresponding to either well ordered (simple cubic, and c.f.c) or random structures with a monodisperse particles distribution, and a random structure with a polydispersity characterized by σ = 0.28. The results for the second derivative of ∆M r with respect to ǫ d and h, C 2 , is displayed in table I. As can be seen, especially for ĥ = ẑ, the mean field approximation or equivalently the linear contribution of the TPT compares well with the MC simulations and in particular for the polydispersity effect.

Since spherical or cubic systems are characterized by D α = 1/3, this first term vanishes in these situations and one is left for the DDI contribution with ∆M r ∝ ǫ 2 d and similarly for χ r . Moreover, still for isotropic systems, we know that the DDI contribution to both M r and χ is negative.

Therefore the magnetization is all the more reduced due to the DDI that the coupling constant ǫ d increases. From the analytical results of the TPT we can calculate the proportionality coefficient relating (χ r (ǫ d ) -χ r (0)) to ǫ 2 d . We have thus compared the MC simulation to the theoretical small ǫ d expansion in the simple cubic structure case and a monodisperse distribution. From this comparison, see figure 3, we can check that the TPT gives an accurate result only for ǫ d < 0.2 as expected. Furthermore, from a description based on the dipolar fields distributions which can be seen as a generalization of the mean field type of approach, ref. [23] have shown also that the dipolar interactions in isotropic systems decrease the magnetization. This decreases is related to the non linearity with respect to the applied field of the non interacting contribution χ (0) r to the susceptibility. Notice that this second type of approach, which remains qualitative in the absence of a theory to deduce the dipolar field distribution, is not restricted to the weak coupling case. Hence, as a general rule, we expect that in an isotropic sample the DDI tend to reduce the magnetization. However, this reasoning does not hold at high fields where the Zeeman term dominates on the DDI and where we expect an approach to saturation, close to what is obtained in the non interacting case deduced from the high field expansion of M (0) r of equation (12), namely

(M r (h) ∼ 1 -1/(β * d * 3 h)).

B. Spherical clusters in the strong coupling case

We now consider, exclusively for spherical clusters, the moderate to strong dipolar coupling case corresponding to the experimental situation of typical coated maghemite NP powders [24],

with ǫ First of all we analyse the linear susceptibility, χ which provides the behavior at low field of the magnetization. Since in our model, with a constant coating layer thickness, ∆, the dipolar coupling constant scales as (d m /d ref ) 3 we expect in the vicinity of h ∼ 0 a reduction of the magnetization higher for large median diameters, where the initial non interacting magnetization M (0) r is higher.

In the quasi monodisperse case, σ = 0.05 we make use of the scaling transformation introduced in equation ( 5) to explicit the effect of the coating layer thickness ∆ on χ by using only one set of simulations for ∆/d ref = 0.20. We checked for ∆ = 0.8 and 2 values of d m the reliability of this scaling transformation (see figure (4)). Therefore, in the quasi monodisperse case we have a rather complete picture of both the effect of the variation of the median diameter, d m and of the distance of closest approach between NP, controlled by the coating layer thickness, ∆. The result for χ determined with σ = 0.05 is displayed and compared to the non interacting case χ (0) in figure (4). As expected, when increasing the particle size and consequently the DDI coupling constant ǫ d , an increasing reduction of χ is obtained. This reduction is of course all the more important that ∆ is small. The important result is that we can get a plateau, which means that χ may becomes particle size independent beyond a threshold value which is, as expected, strongly ∆ dependent. As already mentioned, for random distribution of easy axes, χ does not depend on ǫ K in the superparamagnetic regime, and accordingly this result holds also in the case where the anisotropy is included.

The dependence of χ on both d m and ∆ can be used to deduce the behavior of χ with the NP volumic fraction (or concentration) at fixed value of d m through the relation 3 with ϕ m = ϕ v (∆ = 0). Doing this, in agreement with other MC results [22,39], we get a monotonous decrease of χ with the increase in ϕ, as shown on figure (5). Furthermore this shows that a fit of the NP size on the magnetization curve by using a Langevin function does not hold beyond a critical value ϕ c of the volumic fraction. We can estimate this latter by imposing that χ/χ(ǫ d = 0) is larger than some threshold value say λ, leading the determination of ϕ c through χ(ϕ c )/χ(ǫ d = 0) = λ. The result obtained by using λ = 0.80 is displayed on figure (5). with respect to h, which seems coherent with the interaction fields distribution description [23].

ϕ v = ϕ m (d m /(d m + ∆))
Then we introduce the polydispersity at fixed values of d m . First we consider the case d m /d ref = 1.33, as an example of median diameter located in the plateau region of the χ(d m ) curve. In this case we expect a very weak dependence of the magnetization with respect to the polydispersity in the low field region and this is confirmed by the MC simulations. Indeed, we get only small changes of M r (h) with σ as can be seen in figure (7). The magnetization curves corresponding to σ up to 0.40 are very close to each other for the values of the field for which M r < 0.70; beyond this value, the deviations between the different magnetization curves reflect mainly the approach to saturation where M r (h) ∼

C. Comparison with experiment

We now consider experimental results obtained recently on γ -F e 2 O 3 NP powders samples differing by their size [24]. The experimental protocol for the synthesis is described in [24]. The particles are coated with (5-hydroxy-5,5-bis(phosphono)pentanoic acid) which provides a coating layer of thickness ca 2 nm between particles. As a result of the synthesis method, the standard deviation of the diameter distributions as determined by TEM takes nearly the same value in the 4 samples considered, namely, σ ≃ 0.26. The saturation magnetization M s is found to be in between 61 and 70 emu/g for the distributions characterized by d m = 10, 12, 18 and 21 nm. Although these values are smaller than the bulk value at room temperature (∼ 75 emu/g) the difference is small enough for the spin canting to be neglected in first approximation. It is worth mentioning that the magnetic properties of these NP assemblies have been measured also in diluted solution and, although the possible formation of clusters and/or chains in the presence of the external field cannot be ruled out, this allows for an estimation of the interaction effect. When going from the dispersed samples to the powder ones, we observe both a strong reduction of the magnetization and its weak size dependence in the low field region [24]. According to our simulations, both effects result from the DDI. The DDI induced reduction of M r (H) in the absence of demagnetizing effects is a general simulation result [22,39,40] and a similar trend has been obtained experimentally [41],

and can be deduced from the FC/ZFC measurements in the superparamagnetic regime of either bare or Si coated γ-Fe 2 O 3 NP [33]. Beside its rather weak size dependence the other feature of the experimental reduced magnetization curves, M r (H) in the low field region (see figure (11) is the opening of the hysteresis cycle for the largest sizes beyond d m = 12 nm. These two points are in qualitative agreement with the MC simulations on our model although the opening of the cycle becomes noticeable for larger median diameters (d m ∼ 18 nm) than in MC simulations.

In the present work, we do not compare the experimental magnetization curve in the whole range of field with the results of either a mean field approach or the TPT. In any case the values of the dipolar coupling corresponding to the experimental powders samples (ǫ (ef f ) d ∼ 0.6 to 6.0 when ǫ d ∼ 1 to 8 and the effect of ∆ is taken into account) fall outside of the range of validity of the TPT. Indeed this later is limited to ca ǫ d < 1/6 according to ref. [20], the analytical approach based on TPT of ref. [18] is shown to be very accurate for ǫ d < 0.25 and valid for ǫ d < 0.50 in the monodisperse case and in section III A we found that χ as calculated from the TPT start to deviate from the simulated results at ǫ d ∼ 0.2. Moreover the accounting of the polydispersity is expected to worsen the lack of accuracy of the TPT with the increase of ǫ d .

In figure (12) we compare the experimental and simulated M r (H) for the applied field in the low to intermediate range for d m = 10 nm. The agreement is quite satisfactory up to H = 60 kA/m where M r ≃ 0.7. Then for median diameters d m > 12 nm (d m /d ref > 1.2, we get for both the experimental samples and the MC simulations an opening of the hysteresis cycle. However, as can be deduced from figures (11) and (10) the irreversible field, H irr is found much larger in the MC simulations than in the experiment. Notice that we do not try to map rigorously the MC time scale to the actual measurement time τ m , and this is plays a central role on this point. Therefore concerning the MC simulations we consider in the following the anhysteretic magnetization as defined in equation (10); using this later remains to ignore the hysteresis cycle (i.e. the remanence and the coercive field) or to consider the infinite time scale limit. We compare the experimental M (H) to the simulated ones in figures (13,14) for the median sizes d m 12 and 20 nm respectively.

Given the simplicity of the model which does not include at all any structure in the particles, and the absence of fitting parameter the agreement is qualitatively satisfactory when ǫ Since in this range of fields, the deviation of M r from the non interacting case is governed by the DDI, we can conclude that the strong reduction of the M r (H) variation with respect to H and its relative size independence when compared to the diluted solutions counterpart is the DDI signature. For median sizes larger than d m = 10 nm, the main discrepancy between the simulated and experimental M r (H) curves is the strong non linearity in the very vicinity of H = 0. This is clearly due to the oversimplification of the model in which the particles are uniform single domain ones.

D. Conclusion

In this work we have used MC simulations to investigate both the median size and polydispersity effects on the magnetization curve of densely packed clusters of single domain magnetic NP. An important result is the plateau in the χ(d m ) curve in the quasi monodisperse case for small values of the coating layer ∆, which emphasizes the much reduced size dependence of the M (H) low field dependence in the concentrated systems. Despite of the simplicity of the model, some important features of the experimental M (H) on powder samples are reproduced, especially concerning the DDI signature which occurs principally at low fields and its dependence on the particle size. In order to get a satisfactory agreement with experiments, it appears that the internal polarization structure of the NP should be introduced. Tables struct. cs (a) cfc (b) rand. σ = 0 (c) rand. σ = 0.28 (a) , 1103 (b) , 1054 (c) and 879 (d) particles, characterized by Lx = Ly and Lz = Lx/5. MC simulation with ĥ = x (f ) , ĥ = ẑ (g) ; equation ( 18) with ĥ = x (h) , ĥ = ẑ (i) . (h), (i) The demagnetizing factor entering equation ( 18) is taken from [42] and the moments d * s are taken from the actual diameters distribution of the cluster considered.

Figure captions

1

 1 and the stared lengths are in d m unit. The dimensionless dipolar coupling constant and anisotropy constant are then ǫ d = (d m /d ref ) 3 ǫ (0) d and ǫ K = (d m /d ref ) 3 ǫ (0) K respectively; the reference diameter, d ref can be chosen such that ǫ d (d m = d ref ) ≡ ǫ (0) d = 1 ; the reduced external field h coincides with the usual Langevin variable at temperature T 0 for a monodisperse distribution with d = d m In equation (3), we also introduce the reference external field, H ref for convenience.

  d ref varying from 9.665 nm to 10.422 nm and we use in the following except otherwise mentioned d ref = 10 nm which corresponds to µ 0 M s = 0.488T and H ref = 16.20 kA/m. The anisotropy constant K 1 cannot be taken equal to the bulk effective magnetocrystalline anisotropy constant K b as it is

( 1 )

 1 r (h) can be deduced in the mean field approximation of the magnetization which introduces the DDI contribution to M r from the demagnetizing field and follows from M r (H a ) = M (0) r (H ef f ) where M (0) r corresponds to the non interacting system and H ef f is the effective field

  d ref = 10nm and a coating layer of ∆ = 2.0nm. The median diameter varies from d ref to 2 × d ref and the standard deviation of the distribution ln(d) is taken from σ = 0.05 to represent the quasi monodisperse case to σ = 0.50 to represent a large polydispersity. The importance of σ on the MNP distribution in the clusters is clearly seen on figure[START_REF] Skomski | [END_REF]. Notice that a standard value obtained experimentally is ca. 0.20 ∼ 0.30 which is represented here by σ = 0.28. In the first step we neglect the anisotropy contribution (ǫ K = 0) and focus only on the DDI.

  The magnetization curves in terms of the reduced external field H a /H ref for three values of the median diameter, still for σ = 0.05 is shown in figure (6) and compared to the non interacting diameter distribution weighted Langevin curves. We clearly see the important reduction of the magnetization compared to the non interacting case, and the very weak dependence of the low field behavior with respect to the median diameter which is expected as the considered sizes are either close to the onset of the χ(d m ) curve plateau corresponding to (∆/d ref ) = 0.20 (d m /d ref = 1) or pertain to this later (d m /d ref = 1.33 and 2.00). On the other hand at low external fields the nearly size independence of the magnetization is correlated with a quasi linear behavior of M r (h)

( 1 -

 1 1/(β * d * 3 h)) depends on σ through d * 3 . The deviation from the quasi monodisperse situation over the whole field range occurs for σ ≥ 0.5. Conversely, when the median diameter is taken outside of the χ(d m ) plateau, as is the case for d m /d ref = 1.0 the polydispersity has a noticeable influence on the magnetization as shown in figure (8) for σ ranging from 0.05 to 0.40. In the superparamagnetic regime the MNP anisotropy energy modifies the magnetization curve for intermediate values of the field and leads to a reduction of M r since the moments tend to be pinned in the easy axes directions. Taking into account ǫ K thus reduces further M r for h between the low field region controlled by the DDI and the approach to saturation controlled by the Zeeman energy. In the quasi monodisperse case, the blocking temperature corresponds to that of the median diameter, namely for non interacting particles, k B T b ≃ K 1 v(d m )/25 or equivalently for the reduced blocking temperature 1/β * b ≃ ǫ (0) K (d m /d ref ) 3 /25 leading to 1/β * b ≃ 0.225 for d m /d ref = 1.33. Here we restrict to the room temperature, β * = 1, and we expect the system to be in the superparamagnetic regime even for short times. Indeed for σ = 0.05 our MC simulations confirm the superparamagnetic regime. The result is displayed and compared to the ǫ K = 0 case in figure (9) for d m /d ref = 1.0 and 1.33. As expected, the anisotropy energy does not affect the M (H) curve in the vicinity of H = 0 due to the random distribution of easy axes. Moreover, when d m /d ref = 1, the M (H) curve for intermediate values of the field is only weakly modified by the anisotropy energy while for d m /d ref = 1.33 a noticeable deviation is obtained. The influence of the polydispersity on the magnetization curve when the anisotropy energy is included is shown for d m /d ref = 1.33 on figure (10) for σ ranging from 0.05 to 0.35. Because of the largest particles in the distribution, the system is no more in the superparamagnetic regime for the MC runs considered up to 10 5 MC steps. On the qualitative point of view this is expected since 1/β * b behaves as d * 3 in the absence of DDI and moreover increases with the DDI. As a result an opening of the hysteresis cycle is obtained with remanence magnetization and coercive field increasing with σ as shown in figure (10) in the particular case d m /d ref = 1.33. The magnitude of the hysteresis cycle opening is expected to increase with d m /d ref and is indeed found very weak for d m /d ref = 1.0. The determination of the remanence in terms of d m and the measuring time is beyond the scope of this work; we nevertheless note (see section III C) that the hysteresis cycle opening for large values of d m is in qualitative agreement with experiment.

  particular for the overall variation of M r (H) at low fields (H < 40 kA/m).

Figure 1 Figure 2

 12 Figure 1 Shematic view of the configuration for two particles coated by the layer of thickness ∆/2 at contact. Figure 2 Central part of the clusters corresponding to d m /d ref = 1.33 and σ = 0.05, 0.28 and 0.50 from top to bottom.

Figure 3

 3 Figure 3 Comparison of χ r in terms of ǫ d as calculated from the TPT of Ref. [15], (solid line) and the present MC simulation (symbols) for a spherical cluster of simple cubic structure with N p = 1021 particles and a monodisperse distribution. β * = 1.0.

Figure 4

 4 Figure 4 Linear susceptibility in terms of the median size d m /d ref for different values of the coating layer thickness ∆ in the quasi monodisperse case, σ = 0.05 and β * = 1.0. The two crosses on the ∆/d ref = 0.8 curve correspond to the direct calculation without using the scaling transformation (6). The dotted lines are guides to the eye and the solid line corresponds to the non interacting case.

Figure 5

 5 Figure 5 Susceptibility χ in terms of the reduced volumic fraction, for σ = 0.05, β * = 1.0, and d m /d ref = 2.0 (solid circles); 1.50 (solid squares); 1.25 (upward triangles) and 1.0 (downward triangles). ϕ m = ϕ(∆ = 0). In the present work ϕ m ≃ 0.585. Inset: reduced critical volumic fraction defined as χ(ϕ, ǫ d )/χ(ǫ d = 0) = 0.80 in terms of the median particle size.

Figure 6

 6 Figure 6 Magnetization in terms of the applied field for different values of the median diameter, ∆/d ref = 0.20, β * = 1.0 and N p = 1007 in the quasi monodisperse case, σ = 0.05. The corresponding non interacting curves (diameter distribution weighted Langevin curves) for d m /d ref = 1.0 (long dash), 1.33 (short dash) and 2.0 (dotted line) are displayed for comparison.

Figure 7

 7 Figure 7 Magnetization in terms of the applied field for different values of the ln(d) standard deviation σ and d m /d ref = 1.33, ∆/d ref = 0.20 and β * = 1.0. N p = 1007 (σ = 0.05), 923 (σ = 0.28); 985 (σ = 0.40) and 990 (σ = 0.50). The dotted lines are quides to the eye. The solid lines are the asymptotic limits for σ = 0.05 (bottom) and σ = 0.50 (top).

Figure 8

 8 Figure 8 Magnetization in terms of the applied field for σ = 0.05, 0.20, 0.28 and 0.40 from bottom to top. and d m /d ref = 1.00, ∆/d ref = 0.20, β * = 1.0. The dotted lines are guides to the eyes and the solid lines are the corresponding asymptotic limits for σ = 0.05 and σ = 0.40.

Figure 9

 9 Figure 9 Magnetization curve with the anisotropy energy ǫ (0) K = 2.38, σ = 0.05 and the particles sizes as indicated. The solid lines correspond to ǫ K = 0. β * = 1 and N p = 1007.

Figure 10

 10 Figure 10 Magnetization curve with the anisotropy energy ǫ (0) K = 2.38, d m /d ref = 1.33 and β * = 1 in the polydisperse case with σ = 0.35 (open circles); 0.28 (triangles); 0.20 (squares) compared to the quasi monodisperse case, σ = 0.05 (solid line). Insert: detail of the downward magnetization curve in the vicinity of h = 0, showing the evolution of the remanence and coercivity with σ for σ = 0.35, 0.28, 0.20 and 0.10 (open squares) from top to bottom.

Figure 11

 11 Figure 11 Experimental low field magnetization curves at room temperature for powder samples of γ-Fe 2 O 3 from ref.[24] with polydispersity σ ≃ 0.26 and median sizes as indicated.

Figure 12

 12 Figure 12 Comparison between MC simulations and the experimental magnetization curve at room temperature for d m = 10 nm. The simulation are performed with β * = 1, ǫ (0) K = 2.38 and N p = 923.

Figure 13

 13 Figure 13 Comparison between MC simulations and the experimental magnetization curve at room temperature for d m = 12 nm. The simulations are performed with β * = 1, N p = 980, ǫ (0) K = 2.38 (open circles) or 1.19 (dashed line). The dooted line is a guide to the eye.

Figure 14

 14 Figure 14 Comparison between MC simulations and the experimental magnetization curve at room temperature for d m = 21 nm. The simulations are performed with d m = 20 nm, β * = 1, N p = 998, ǫ (0) K = 2.38 (open circles) or 1.19 (dashed line).

Figure 1 :

 1 Figure 1: Shematic view of the configuration for two particles coated by the layer of thickness ∆/2 at contact.

Figure 2 :Figure 3 :

 23 Figure 2: Central part of the clusters corresponding to dm/d ref = 1.33 and σ = 0.05, 0.28 and 0.50 from top to bottom.

Figure 4 :

 4 Figure 4: Linear susceptibility in terms of the median size dm/d ref for different values of the coating layer thickness ∆ in the quasi monodisperse case, σ = 0.05 and β * = 1.0. The two crosses on the ∆/d ref = 0.8 curve correspond to the direct calculation without using the scaling transformation(6). The dotted lines are guides to the eye and the solid line corresponds to the non interacting case.

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Susceptibility χ in terms of the reduced volumic fraction, for σ = 0.05, β * = 1.0, and dm/d ref = 2.0 (solid circles); 1.50 (solid squares); 1.25 (upward triangles) and 1.0 (downward triangles). ϕm = ϕ(∆ = 0). In the present work ϕm ≃ 0.59. Inset: reduced critical volumic fraction defined as χ(ϕ, ǫ d )/χ(ǫ d = 0) = 0.80 in terms of the median particle size.

Figure 8 :

 8 Figure 8: Magnetization in terms of the applied field for σ = 0.05, 0.20, 0.28 and 0.40 from bottom to top. and dm/d ref = 1.00, ∆/d ref = 0.20, β * = 1.0. The dotted lines are guides to the eyes and the solid lines are the corresponding asymptotic limits for σ = 0.05 and σ = 0.40.

Figure 9 :

 9 Figure 9: Magnetization curve with the anisotropy energy ǫ (0) K = 2.38, σ = 0.05 and the particles sizes as indicated. The solid lines correspond to ǫK = 0. β * = 1 and Np = 1007.

Figure 10 :

 10 Figure 10: Magnetization curve with the anisotropy energy ǫ (0) K = 2.38, dm/d ref = 1.33 and β * = 1 in the polydisperse case with σ = 0.35 (open circles); 0.28 (triangles); 0.20 (squares) compared to the quasi monodisperse case, σ = 0.05 (solid line). Insert: detail of the downward magnetization curve in the vicinity of h = 0, showing the evolution of the remanence and coercivity with σ for σ = 0.35, 0.28, 0.20 and 0.10 (open squares) from top to bottom.

Table I :

 I Second derivative C2 = ∂ 2 (∆Mr)/∂ǫ d ∂h at ǫ d = 0 and h = 0 for a primatic cluster of 1024

					(d)
	ϕ v	π/6	0.74	0.525	0.580
	C 2 (f ) C (g) 2 C (h) 2 C (i) 2	0.2775 0.3367 -0.4745 -0.6109 0.2217 0.3133 -0.4433 -0.6265	0.2534 -0.4390 0.2217 -0.4433	1.7934 -3.1506 1.6173 -3.2354

Acknowledgements

This work was granted access to the HPC resources of CINES under the allocation 2012-c096180 made by GENCI (grand Equipement National de Calcul Intensif).