Chao Chen
email: cschaochen@cqu.edu.cn

Shenle Pan
email: shenle.pan@mines-paristech.fr

Using the Crowd of Taxis to Last Mile Delivery in Ecommerce: a Methodological Research

Keywords: Last Mile Delivery, Crowdsourcing, Taxi Trajectory Data Mining, Freight Transport, City Logistics

Using the Crowd of Taxis to Last Mile Delivery in E-commerce: a Methodological Research

Introduction

In E-commerce environment, Last Mile Delivery (hereafter LMD) is the problem of transport planning for delivering goods from e-retailers' hub to the final destination in the area, for example the end consumers' home, see [START_REF] Lee | Winning the last mile of e-commerce[END_REF] and [START_REF] Punakivi | Solving the last mile issue: reception box or delivery box[END_REF]. Speed and cost are the two crucial success factors to LMD. Faster shipping while with lower cost is the major challenge; nevertheless, it is also a paradox to a certain extend. Indeed, when customers are given a choice between fast and cheap delivery, most of them choose the cheap one, observed by a recent report [START_REF] Gibson | The state of the retail supply chain: Essential Findings of the Fifth Annual Report[END_REF]. The report also infers that that lowcost, speedy two-day delivery corresponds to most customers' expectation, opposite to the one-day delivery policy pursued by giant e-retailers such as Amazon and Alibaba etc. This fact may open up new opportunities to innovative freight transport models [START_REF] Meyer-Larsen | LogisticsArena-A Platform Promoting Innovation in Logistics[END_REF] for LMD aiming at reducing delivery cost while respecting shipping time, nevertheless, not necessarily aiming at minimizing shipping time. Being our topic here, crowdsourcing is one of such solutions getting more and more attention [START_REF] Paloheimo | Transport reduction by crowdsourced deliveries -a library case in Finland[END_REF].

In the literature, crowdsourcing has been usually seen as "an interesting idea" for freight transport before seriously moving to real applications. Despite the existence of some internet-based services, scientific research, especially methodology for application is rare in the literature [START_REF] Paloheimo | Transport reduction by crowdsourced deliveries -a library case in Finland[END_REF]. This paper aims to fill this gap by providing a methodological approach of applying the crowdsourcing solution proposed and to assess its performance. In this paper, the crowd studied is taxi fleet in city, supported by a transport network composed by road network and customer self-pickup facilities such as 24 hours shops in city, named as TaxiCrowdShipping system. The system relies on a two-phase decision model, first offline taxi trajectory mining and second online package routing -taxi scheduling. As the first stage of our research, this paper introduces the objective and the framework of the TaxiCrowdShipping system, as well as to define the function of system.

The reminder of the paper is organised as follows. Section 2 consists of a relevant literature review. Then, Section 3 presents the TaxiCrowdShipping system. After introducing some basic concepts and assumptions, we focus on the two steps approach being the decision support tool for the system. Some expected results are also discussed. Finally, Section 4 concludes this paper by giving some research perspectives for future works.

Related Works

Recently some innovative solutions have been studied for city logistics and LMD in E-commerce environment, for example those involved in our study like interconnected city logistics enabled by Physical Internet [START_REF] Sarraj | Analogies between Internet network and logistics service networks: challenges involved in the interconnection[END_REF][START_REF] Crainic | Physical Internet Enabled Interconnected City Logistics[END_REF], self-service parcel station (e.g., DHL PackStation, LaPoste Pickup Station etc.), new tools for LMD (bicycle, motor, electric vehicle etc.), Smart city logistics [START_REF] Neirotti | Current trends in Smart City initiatives: Some stylised facts[END_REF], and crowdsourced delivery [START_REF] Paloheimo | Transport reduction by crowdsourced deliveries -a library case in Finland[END_REF]. Due to the space limitation, here we focus on the works regarding crowdsourcing in Freight Transport.

Being firstly discussed in [START_REF] Howe | The rise of crowdsourcing[END_REF], crowdsourcing has been increasingly studied as a solution to freight transport. It can be simply defined as "outsourcing a task to the crowd via open call" [START_REF] Howe | The rise of crowdsourcing[END_REF]. One the practice side, it occurs mainly in the form of internet-based services for example imoveit.co.uk and zipments.com, where the crowd is undefined. Thus, both professional (e.g., carriers) and non-professional (e.g., inhabitants) service providers may answer the calls. In 2014 Amazon has launched a project to explore taxi deliveries in San Francisco and Los Angeles 1 . The idea is similar to our study, though, their methodology and results are not yet published to our knowledge. Moreover, the package deliveries are completed by ordering free taxis, while our proposed solution leverages the hitchhiking rides provided by occupied taxis when they are sending passengers, thus our solution is more green and economic. On the scientific side, only few relevant works in the area of logistics can be found. A case study of applying crowdsourcing to library deliveries in Finland is conducted in [START_REF] Paloheimo | Transport reduction by crowdsourced deliveries -a library case in Finland[END_REF]. They study a system called PiggyBaggy to assess the sustainability and adaptability of such solution. A taxi-based solution for the waste-collecting or product return problem (i.e., reverse logistics) in metropolitan area is discussed in [START_REF] Pan | A crowdsourcing solution to collect e-commerce reverse flows in metropolitan areas[END_REF], without considering goods delivery. Some other relevant works can be also found in the area of data science. Data scientists are mainly interested at mining the taxi trajectory data to understand the city dynamics, and developing various smart services for taxi drivers, passengers, as well as the city planners [START_REF] Castro | From taxi GPS traces to social and community dynamics: A survey[END_REF][START_REF] Zheng | Trajectory data mining: an overview[END_REF]. However, almost all the current research related to taxi data mining focuses on the people or public transport [START_REF] Chao | B-Planner: Planning Bidirectional Night Bus Routes Using Large-Scale Taxi GPS Traces[END_REF][START_REF] Liu | Exploiting Heterogeneous Human Mobility Patterns for Intelligent Bus Routing[END_REF], little attention has been paid to freight transport.

From the literature we can see that crowdsourcing in freight transport usually occurs in the form of internet-based services in practice, and it is usually investigated via case study in the literature. Methodology for application is not well addressed. Besides, no attention has been paid to crowd selection or definition. People in city are often regarded as eligible crowd. Following the previous work [START_REF] Pan | A crowdsourcing solution to collect e-commerce reverse flows in metropolitan areas[END_REF] dealing with reverse flows, this paper focus on a methodology approach for the LMD problem, where logistics constraints and decision model are different.

3

TaxiCrowdShipping System

Basic Concepts and Assumptions

To ease the description, we define the related concepts based on Fig. 1, and also make some assumptions.

Definition 2. (Taxi Trajectory)

A taxi trajectory is a sequence of time-stamped GPS points. Each GPS point p i =(t i , x i , y i , ind i ) consists of a time-stamp t i , a longitude x i , a latitude y i , and an indicator ind i showing whether the taxi is occupied or not. A pick-up point is a special GPS point, with the indicator changing from 0 to 1 (the redcoloured circle in Fig. 1); while a drop-off point is the one with the indicator changing from 1 to 0 (the green-coloured circle in Fig. 1). Thus we can further define a passenger-delivery trajectory is the GPS sequence from the pick-up point to the followed drop-off one (the red dashed line in Fig. 1); a passenger-hunting trajectory is the GPS sequence from the drop-off point to the followed pick-up point (the green dashed line in Fig. 1).

Definition 3. (Package Pickup Station)

A package pickup station is a Point of Interest (POI hereafter) near roadside that is responsible for storing packages waiting for consumer pickup, (the star in Fig. 1). Here we select 24-hour opening convenience stores near roadside as the package pickup stations.

Definition 4. (Package Delivery Request)

A package delivery request is defined as a triple <o p , d p , t p >, where o p and d p refer to the origin and the destination of the package respectively, and t p refers to the time when the user submits the request. The request is generated by users who need the package express delivery service.

Definition 5. (Real-time Taxi Ordering Request)

A real-time taxi ordering request is defined as a triple <o t , d t , t t >, where o t and d t refer to the passenger's origin and the destination respectively; t t refers to the time when the passenger submits the request. The request is made by passengers who need taxi service. Assumption 1. All selected POI is open 24/7, without capacity issue and with good accessibility to taxi. Assumption 2. The taxi drivers are willing to accept the assigned package delivery tasks.

Assumption 3. The package can be trackable. Since the birth time, the package is either stored at the pickup station or carried by the scheduled taxi. Each pickup station is authorized and has a unique ID; each taxi is registered in taxi management department and also has a unique ID.

Problem Description

To help understand how our proposed solution works to handle the LMD, we intentionally design a simple running example. Suppose in Fig. 1 the leftmost and rightmost star are the origin and the destination of the package respectively. After the generation of package delivery request, there happens to be a passenger who makes a real-time taxi ordering request, intending to go to the same destination. At that time, we can assign the package delivery task to the taxi which has responded the passenger's re-quest. Finally, the package will be also delivered, with a hitchhiking ride provided by the taxi while sending the passenger. The solution can be featured as economic and eco-friendly one since it almost does not incur extra labour cost and energy as well as CO 2 emissions.

Accordingly, the taxi-based crowdsourcing solution to LMD consists of the shortest path finding problem for packages and the scheduling problem for taxis; and it can be described as follows.

Given:

-A road network and a set of package pickup stations in the studied city; -A set of taxi trajectory data in the studied city in history (e.g., last month); -A set of package delivery requests, and a set of real-time generated taxi ordering requests. Note that these requests come in stream.

Objective:

For a given package delivery request, find its optimal package delivery path which can minimize the total package delivery time (i.e. maximize the delivery speed). Once determined the path, we can schedule the next coming taxi to delivery the package having the same destination. Note that in one of the scenarios in this study the path can be re-planned according to the Real-time Taxi Ordering Request.

Constraints:

Only taxis which response the taxi ordering requests after the package delivery request can be scheduled.

Once a taxi is involved into a delivery task, it can be available again to be scheduled to participate only after completing the current task (i.e. sending the package to the predefined pickup station). In another word, a taxi can carry at most one package when sending passengers.

3.3

The TaxiCrowdShipping System

The TaxiCrowdShipping system contains two components, i.e., the Offline-Trajectory Mining and the Online-Package Routing respectively, which will be detailed as follows and shown in Fig 2.

Offline-Trajectory Mining

The objective of the Offline-Trajectory Mining is to estimate the direct package delivery time from one pickup station to another one, by taking a single hitchhiking ride. The time cost mainly includes two parts: the time cost on waiting for the hitchhiking rides, which is related to the frequency of taxi rides, and the time spent on driving on the roads. Here, we propose a two-step procedure to estimate the time cost.

Step 1: From Trajectory Data to Passenger Flow. From the given taxi trajectory data, it is not difficult to compute the passenger flow between any two pickup stations. Specifically, to compute the passenger flow from cs i to cs j , the trajectories meeting Eqns. 1∼2 will be counted. Then, the passenger flow from cs i to cs j is just the number of trajectories satisfying the requirements. Note that there may be no passenger flow between some pickup station pairs, and the passenger flow is different at different time slots. For the purpose of future research, we divide the time into three time slots for a day in advance, i.e., night-time hours, day-time hours and rush hours. Step 2: From Passenger Flow to Time Cost. To estimate the time cost, we need to estimate two parts, i.e. the waiting time and the driving time. The waiting time is defined as the time cost on waiting for the suitable hitchhiking ride event of passenger taking taxis, to help deliver a package from cs i to cs j directly (with no transhipment). Here, we employ the Non-Homogeneous Poisson Process (NHPP) to model the behaviour of passenger taking taxis [START_REF] Qi | How Long a Passenger Waits for a Vacant Taxi--Large-Scale Taxi Trace Mining for Smart Cities[END_REF]. According to the passenger flow, we can estimate the waiting time of packages at different time slots at the pickup stations. Under the Poisson hypothesis within a time slot, we could derive the probability distribution of the waiting time for the next suitable hitchhiking ride event (i.e. t next , the event of a passenger taking taxi from cs i to cs j ), which can be expressed in Eq. 3:

Here N(t) represents the number of event occurring within t, and . Then the probability density function (pdf) of t next is just the derived function of P{}, as can be seen in Eq. 4.

Thus, we can deduce the expectation of t next (i.e. the waiting time for the hitchhiking ride event occurring):

Note that in the model is the frequency of passenger taking taxis from cs i to cs j (i.e. the passenger flow from cs i to cs j ), which can be easily estimated by the Eq. 6.

where is the average number of passengers taking taxis from cs i to cs j during the studied time slot in the observed days; ∆T is the time duration of the that time shot.

Therefore, the waiting time from cs i to cs j is:

For each passenger-delivery ride from cs i to cs j , it is easy to derive its time spent on driving on the roads. The driving time is simply the average one of all such rides, as shown in Eq. 8.

where N is the number of passenger-delivery rides during the studied time slot in the observed days. t e -t s is the time cost of the corresponding taxi ride.

Finally, the time cost is just the sum of waiting time and driving time, as shown in Eq. 9.

Note that the time cost will be +∞ if there was no passenger flow on the respected pickup station pair.

Online-Package Routing and Taxi Scheduling

The objective of the Online-Package Routing and Taxi Scheduling is to schedule the specific taxis to help delivery the packages with the determined optimal path, according to the real-time coming taxi ordering requests. Here, we also propose a two-step procedure to complete, detailed as follows.

Step 1: Find the Optimal Pickup Station Sequence. For a package delivery request, with the estimated time cost values in the last component, it is trivial to find the best pickup station sequence from the origin to the destination of the package, in terms of the total time cost, by applying the classical shortest path finding algorithms.

Step 2: Schedule the Taxis. After obtaining the optimal pickup station sequence for a package delivery request, we schedule the taxis according to the real-time taxi ordering requests. In more detail, from the origin of the package to the followed pickup station in the optimal pickup station sequence, we wait for the taxi which will pick up a passenger at the origin, heading to the followed station (the information is included in the real-time taxi ordering requests), and assign the package delivery task to that taxi. After that, the origin of the package will be also updated. The procedure will be repeated until the package arrives at its destination.

Expected Results

Following the framework proposed here some results are expected in the next steps. First, we will conduct study to assess the implementability of the TaxiCrowdShipping system proposed. A large city in China, namely Hangzhou city is selected to be the test field, thanks to some available data sets there such as Open data of taxi trajectory, may of city shops' and road network etc. However, the data of package delivery request is still to be completed. Second, a set of algorithms for package routing and taxi scheduling problem will be developed and examined. Then a set of scenarios is expected to the study in order to assess the performance of the system as well as its sensibility of setting.

Conclusion

In this paper we aim to propose a methodological approach of applying crowdsourcing solution to Last Mile Delivery (LMD) problem in city logistics. To this end, a system called TaxiCrowdShipping is proposed, whose objective is to use the taxi fleet and shops in city for LMD purpose. The framework of such system is discussed in this paper, as well as the two-phase decision model. Except the expected results discussed above, this study opens some research perspectives. First, the Physical Internet-based active container described in [START_REF] Sallez | On the Activeness of Physical Internet containers[END_REF] may be adapted to this study. Suppose that such containers are able to actively publish a delivery request on web (as a passenger calls a taxi), the taxi scheduling could be more efficient and responsive. It could also provide a good real time tracking and tracing technique to the system. Second, the study can be extended to the automated self-service parcel station implantation problem. Coupling such station placement and crowdsourcing solution is still rarely studied in the literature.

Fig. 1 .

 1 Fig. 1. Illustration of some basic concepts.

Fig. 2 .

 2 Fig. 2. The TaxiCrowdShipping system architecture.

where

  Tr i .o and Tr i .d are the original and destination points of Tr i , respectively; loc(•) gets the latitude and longitude location of the given pickup station; δ is a userspecified parameter. Ddist(a • b) calculates the driving distance from point a to b.

http://www.engadget.com/2014/11/05/amazon-is-exploring-taxi-deliveries-in-san-franciscoand-los-ang/


Julien Tranchida

Pascal Thibaudeau

Stam Nicolis

Colored-noise magnetization dynamics: from weakly to strongly correlated noise

Keywords: Stochastic magnetization dynamics, Non-Markovian effects, Multiscale models

Statistical averaging theorems allow us to derive a set of equations for the averaged magnetization dynamics in the presence of colored (non-Markovian) noise. The non-Markovian character of the noise is described by a finite auto-correlation time, τ , that can be identified with the finite response time of the thermal bath to the system of interest. Hitherto, this model was only tested for the case of weakly correlated noise (when τ is equivalent or smaller than the integration timestep). In order to probe its validity for a broader range of auto-correlation times, a non-Markovian integration model, based on the stochastic Landau-Lifshitz-Gilbert equation is presented. Comparisons between the two models are discussed, and these provide evidence that both formalisms remain equivalent, even for strongly correlated noise (i.e. τ much larger than the integration timestep).

I. INTRODUCTION

The study of the effects of thermal fluctuations on the magnetization dynamics, by way of stochastic simulations, is a topic that has received much attention [START_REF] Brown | Thermal fluctuations of a single-domain particle[END_REF], [START_REF] Coffey | Thermal fluctuations of magnetic nanoparticles: Fifty years after Brown[END_REF]. What has also been the subject of numerous studies, in the field of coarsegrained simulations [START_REF] Thibaudeau | Frequency-dependent effective permeability tensor of unsaturated polycrystalline ferrites[END_REF], is the development of multiscale models, able to simulate the average magnetic behaviour of a spin assembly under thermal stresses [START_REF] Garanin | Dynamics of an ensemble of single-domain magnetic particles[END_REF], [START_REF] Kazantseva | Towards multiscale modeling of magnetic materials: Simulations of FePt[END_REF]. But in both contexts, the coupling between the thermal bath and the spin assembly has mainly been taken into account within the framework of the white-noise limit, and few studies [START_REF] Bose | Correlation effects in the stochastic Landau-Lifshitz-Gilbert equation[END_REF], [START_REF] Atxitia | Ultrafast spin dynamics: the effect of colored noise[END_REF] have considered the effects of a finite auto-correlation time for the response of the thermal bath, i.e. the case of colored noise.

However, the simulation of ultra-fast magnetic experiments [START_REF] Kirilyuk | Ultrafast optical manipulation of magnetic order[END_REF], [START_REF] Vahaplar | All-optical magnetization reversal by circularly polarized laser pulses: Experiment and multiscale modeling[END_REF] is now able to probe time-scales of the dynamics of comparable order of magnitude [START_REF] Rossi | Dynamics of magnetization coupled to a thermal bath of elastic modes[END_REF], [START_REF] Vaterlaus | Spin-lattice relaxation time of ferromagnetic gadolinium determined with time-resolved spin-polarized photoemission[END_REF] to the autocorrelation time of the external random field. In these situations, a colored form for the noise, introducing a finite auto-correlation time, can provide a better description for the random variables of the bath [START_REF] Bose | Correlation effects in the stochastic Landau-Lifshitz-Gilbert equation[END_REF], [START_REF] Atxitia | Ultrafast spin dynamics: the effect of colored noise[END_REF].

The bath variables can be defined by random vectors ω whose components ωi are drawn from a Gaussian distribution with zero mean and second-order autocorrelation functions:

ωi (t)ω j (t ′ ) = D τ δ ij exp - |t -t ′ | τ (1)
with τ the auto-correlation time, D the amplitude of the noise. Wick's theorem allows us to express all other autocorrelation functions in terms of these one-and two-point functions [START_REF] Zinn-Justin | Quantum field theory and critical phenomena[END_REF]. When the auto-correlation time vanishes (τ → 0), the well-known white-noise limit is recovered:

ωi (t)ω j (t ′ ) = 2Dδ ij δ (t -t ′ ).
Corresponding author: J. Tranchida (email: julien.tranchida@cea.fr).

In previous work [START_REF] Tranchida | Closing the hierarchy for non-markovian magnetization dynamics[END_REF], statistical averaging theorems [START_REF] Shapiro | Formulae of differentiation and their use for solving stochastic equations[END_REF] were applied to a non-Markovian form of the stochastic Landau-Lifshitz-Gilbert equation (sLLG), and they allowed us to derive a new hierarchy of equations for the simulation of averaged magnetization dynamics. Under Gaussian assumptions, the hierarchy of equations was closed, and a new model, presenting both transverse and longitudinal forms of damping, was obtained. This model is denoted as the dynamical Landau-Lifschitz-Bloch model (dLLB). A first benchmark of tests was performed, and the validity of our closure assumptions was checked to hold for weakly correlated noise (i.e. τ → 0), and with a constant applied magnetic field. However, the case of strongly correlated noise was left for future investigations-in the present paper we discuss our first results along this path.

In section II, the non-Markovian sLLG model is presented, and compared with other works based on non-Markovian stochastic magnetization dynamics. Then, some numerical experiments are presented in order to assess some aspects of this model, and to probe the impact of the auto-correlation time τ on the magnetization dynamics.

Section III introduces our dLLB model. The results of numerical experiments, performed to assess the range of validity of our closure assumptions, for different domains of values of τ , i.e. beyond Markovian situations are discussed.

We finally conclude, in section IV, on the relevance of our dLLB model, when the auto-correlation time for the bath variables is no longer negligible, and present some ideas about the situations where our closure assumptions may break down and, therefore, when the validity of our model attains its limits.

II. NON-MARKOVIAN SLLG MODELLING

The Einstein summation convention is used throughout, with Latin indices standing for vector components, ǫ ijk is the anti-symmetric Levi-Civita pseudo-tensor and s i denotes the component of the normalized magnetization vector (| s| = 1).

Up to a renormalization over the noise [START_REF] Mayergoyz | Nonlinear magnetization dynamics in nanosystems[END_REF], the stochastic Landau-Lifshitz-Gilbert (sLLG) equation can be written as follows:

∂s i ∂t = 1 1 + λ 2 ǫ ijk s k [ω j + ωj -λǫ jlm ω l s m ] (2)
with ω j the precession frequency assumed independent of the spin state and ωj describes the (fluctuating) contribution of the thermal bath.

When the ω become non-Markovian random variables, drawn from a Gaussian distribution, defined by the autocorrelation functions (1), the sLLG equation (eq. 2) can be recast into the following set of equations [START_REF] Fox | Fast, accurate algorithm for numerical simulation of exponentially correlated colored noise[END_REF], [START_REF] Berdichevsky | Stochastic resonance in linear systems subject to multiplicative and additive noise[END_REF]:

∂s i ∂t = 1 1 + λ 2 ǫ ijk s k [ω j + ωj -λǫ jlm ω l s m ] (3)
∂ ωj ∂t = - 1 τ (ω j -ξ j ) (4)
and where ξ are new Gaussian random variables, such that:

ξ j (t) = 0 and ξ i (t)ξ j (t ′ ) = 2Dδ ij δ (t -t ′ ).
The bath variables are now, themselves, solutions of a stochastic differential equation. This, indeed, is what introduces inertial effects, that are described by the auto-correlation time τ .

An interesting feature of this new set of equations is that, for the case of no transverse damping (λ = 0 in eq. 3), we immediately recover the formalism developed by Myasaki and Seki (MS) [START_REF] Miyazaki | Brownian motion of spins revisited[END_REF]. In this way, the set of equations [START_REF] Thibaudeau | Frequency-dependent effective permeability tensor of unsaturated polycrystalline ferrites[END_REF][START_REF] Garanin | Dynamics of an ensemble of single-domain magnetic particles[END_REF] can be seen as a generalization of their work in the case of finite transverse damping.

However, their equation on the bath variables takes a slightly different form than in our work, because they introduced a phenomenological form of damping, that depends on the magnetization state (and with a spin coupling constant χ). The dynamics of ω is indeed damped by a factor that depends on the value of the magnetization, and therefore, the corresponding distribution (for the bath variable) has no reason to remain Gaussian (especially in the case of a non-trivial Hamiltonian for the spin system). In our case, by construction, it is ensured that the bath variables ωi are, both, always drawn from a Gaussian distribution, and remain sensitive to the inertial effects of the memory kernel, defined by the color of the noise.

This set of equations ( 3) and ( 4) is then integrated using a third-order Omelyan's algorithm [START_REF] Omelyan | Symplectic analytically integrable decomposition algorithms: classification, derivation, and application to molecular dynamics, quantum and celestial mechanics simulations[END_REF], [START_REF] Ma | Large-scale simulation of the spin-lattice dynamics in ferromagnetic iron[END_REF] that respects the symplectic properties of equation [START_REF] Thibaudeau | Frequency-dependent effective permeability tensor of unsaturated polycrystalline ferrites[END_REF]. Regarding the noise, this system is formally integrated as ωi (t) ≡ e tLω ωi (0) = t 0 e -(t-u)/τ ξ i (u)du with L ω is the Liouville operator corresponding to Eq.( 4). Because the Liouville operators of the spin and the noise commute, an integration of the colored system is performed by the combination of the two operators as e tL (s i (0), ωi (0)) = e tLs e tLω (s i (0), ωi (0)) = (s i (t), ωi (t)).

It should be stressed that this scheme indeed incorporates the Stratonovich Ansatz, whose theoretical foundations have been studied by Aron et al. [START_REF] Aron | Magnetization dynamics: path-integral formalism for the stochastic Landau-Lifshitz-Gilbert equation[END_REF] for the case of the sLLG.

An effective averaging model is then constructed: the integration of the set of stochastic equations (3,4) is performed one thousand times, each time for a different realization of the random variables ξ i . The results for the spin realizations are stored, and statistical averages thereof are computed. This methodology presents, of course, some disadvantages: first of all, that it is purely computational; also, that it requires significant storage for the configurations. However, these are more than compensated for by the fact that it provides an exact formulation and resolution of the Langevin-like stochastic equation (our non-Markovian sLLG equation). Fig. 1 presents some results directly obtained with the effective averaging model. Fig. 1 highlights the following interesting features of our system: A first, very important, result is that, for the case N = 1, that plots two stochastic simulations, without any averaging effects, increasing the auto-correlation time, τ , leads to stronger inertial effects, that can be directly seen from the curves of s x . Indeed, from fig. 1, we can deduce that this variation does not have a strong impact on the amplitude of the fluctuations of s x , but tends to smoothen them out. Some other features can be observed, concerning the square of the averaged magnetization norm. First, when N = 1, the norm of the magnetization remains perfectly constant and equal to unity (without any renormalization necessary over the time steps). This is a good test of the validity of the integration scheme that has been applied, and therefore make all our procedures (eqs. ( 3) and ( 4), and the integration scheme) a good framework for non-Markovian atomistic simulations. Moreover, when the number N , of realizations increases (which is equivalent to a better approximation of the statistically averaged quantities), the square of the averaged magnetization norm stabilizes below unity. This constitutes a simple and very direct recovery of the equivalence between the presence of a longitudinal form of damping, and statistical averaging in a thermalized magnetic system. Therefore, the behavior in time of the average magnetization dynamics (that our dLLB model is expected to track) can be used to realize direct comparisons between the two models, for the same set of parameters.

III. TESTING THE LIMITS OF THE AUTO-CORRELATION TIME EFFECTS

The set of statistically averaged equations, called dLLB model, is now introduced. With • denoting the statistical average over the noise, one has the following set of 21 equations:

∂ s i ∂t = 1 1 + λ 2 [ǫ ijk ω j s k + ǫ ijk ωj s k -λǫ ijk ǫ jlm ω l s k s m ] (5) ∂ ωi s j ∂t = - 1 τ ωi s j + 1 1 + λ 2 ǫ jkl ω k ωi s l + D τ ǫ jil s l -λǫ jkl ǫ kmn ω m ( s l ωi s n + s n ωi s l ) (6)
∂ s i s j ∂t = 1 1 + λ 2 ǫ jkl ω k s i s l +ǫ jkl ( s i ωk s l + s l ωk s i ) -λǫ jkl ǫ kmn ω m ( s i s l s n + s l s i s n + s n s i s l -2 s i s l s n ) + (i ↔ j) (7)
Equations ( 5) and ( 7) are strongly related to those derived by Garanin et al. [START_REF] Garanin | Dynamics of an ensemble of single-domain magnetic particles[END_REF], and obtained under slightly different assumptions. Equation ( 6) is deduced from the application of the Shapiro-Loginov formula [START_REF] Shapiro | Formulae of differentiation and their use for solving stochastic equations[END_REF], and defines the time evolution of the mixed moments of the thermal bath variables and the spin system variables. This model was only shown to be valid for the simulation of averaged magnetization dynamics in interaction with a constant magnetic field, and when the auto-correlation time can be considered to be "small" (i.e. close to the Markovian limit) [START_REF] Tranchida | Closing the hierarchy for non-markovian magnetization dynamics[END_REF]. Under those conditions, the Gaussian assumptions that have been applied to the moments (and that allowed us to close the hierarchy of moments and to integrate our dLLB model) hold.

We now want to probe the range of validity of our closure assumptions for other values of the auto-correlation time τ . This is possible via a direct comparison between our dLLB set of equations and the effective averaging model. Figures 2 and3 plot these comparisons for two very different values of the auto-correlation time τ .

For weakly correlated noise (τ = ∆t), fig. 2 shows an excellent agreement between the two models. The value of the autocorrelation time, τ was then increased, while keeping the value of the amplitude of the noise, D, fixed. As shown in fig. 3, the results of the two models only begin to diverge one from the other for "strongly" auto-correlated noise (τ = 10 -2 ns).

Several orders of magnitude separate the two values of the auto-correlation time τ ; nonetheless, we remark that our results provide evidence that our dLLB model and the related closure assumptions remain valid, also, for strongly non-Markovian noise. Both stochastic and deterministic curves of fig. 2 and3 are strongly smoothed by the averaging. Therefore, the difference in terms of inertial effects, between the strongly and slightly Markovian cases (when τ = 10 -2 , or τ = 10 -7 , respectively) is barely observable (unlike it was on fig. 1, for the situations of small realization numbers, N=1 and 10). But an interesting feature, that is recovered from the comparison of the two graphs, is that increasing the auto-correlation time leads, on average, to an enhancement of longitudinal damping effects. An equivalent result was already observed by Atxitia et al. [START_REF] Atxitia | Ultrafast spin dynamics: the effect of colored noise[END_REF] in atomistic non-Markovian simulations, performed with the Myasaki and Seki (MS) model.

IV. CONCLUSION

An integration method for the non-Markovian sLLG equation was presented. This scheme allows us to construct an effective averaging model that solves the stochastic Langevinlike equation, for a fixed realization of the non-Markovian noise process and construct the magnetization correlation functions from the many realizations that were drawn-the only limitation is storage for the configurations.

A first set of numerical experiments was performed with this effective averaging model (results plotted on fig. 1). It allowed us to test the consistency of our integration scheme, and to probe the effects of the auto-correlation time τ on the magnetization dynamics. Inertial effects were clearly observed with the variations of τ , and the influence of statistical averaging on the magnetization norm was also recovered (equivalence between statistical averaging and longitudinal form of damping in a magnetic system connected to a thermal bath).

Then, comparisons between the numerical results obtained from this model and our previously derived dLLB set of equations were carried out. These simulations lead us to two main conclusions. First, it was shown that, even across large intervals of values of the auto-correlation time τ (from weakly correlated noise, τ = 10 -7 ns, to clearly non-Markovian cases, 10 -2 ns), the closure assumptions of our hierarchy of equations hold, and our dLLB model remains valid. The application of these techniques also allowed us to recover the result that increasing the auto-correlation time for the bath variables not only enhances their inertial effects, but can also lead to stronger longitudinal damping effects, that can be measured in the average magnetization dynamics.

In future work [START_REF] Thibaudeau | Non-markovian magnetization dynamics for uniaxial nanomagnets[END_REF], it will be shown that what does affect the closure assumptions and the validity of our model are the magnetic interactions that are taken into account, for example, the contribution of a uniaxial anisotropy, or an exchange interaction.

In this case, it was observed that cumulants of the spin and bath variables, assumed to be equal to zero under Gaussian assumptions, are no longer taking negligible values (like the third order cumulants of the spin, for example). Therefore, non-Gaussian closure assumptions are required, and a particular attention has to be taken for the treatment of the mixed (noise and spin) moments. The stochastic approach to chaotic dynamical systems [START_REF] Nicolis | Closing the hierarchy of moment equations in nonlinear dynamical systems[END_REF] may provide valuable insights to this end.

Figure 1 .

 1 Figure 1. Comparisons between simulations performed with the effective stochastic model, for two different values of the auto-correlation time: τ = 10 -7 above, and τ = 10 -2 below. The red curves plot the dynamics of a component of the averaged magnetization vector ( sx ) in a constant magnetic field in the z-direction, and the blue curves plot the square of the averaged magnetization norm, for different repetition values N (N=1, 10 and 1000). Parameters of the simulations: {D = 1 rad.GHz; λ = 0.1; ω = (0, 0, ωz), with ωz = 176 rad.GHz; timestep ∆t = 10 -7 ns}. Initial conditions: s(0) = (1, 0, 0).

Figure 2 .Figure 3 .

 23 Figure2. Relaxation of the averaged magnetization dynamics in a constant magnetic field in the z-direction for "weakly" correlated noise, τ = 10 -7 ns. The solid lines plot the stochastic case (10 3 repetitions), dots are for the dLLB model. Parameters of the simulations: {D = 1 rad.GHz; λ = 0.1; ω = (0, 0, ωz), with ωz = 176 rad.GHz; timestep ∆t = 10 -7 ns}. Initial conditions: s(0) = (1, 0, 0).

ACKNOWLEDGMENTS

JT acknowledges financial support through a joint doctoral fellowship "Région Centre-CEA" under the grant agreement number 00086667.


Pascal Thibaudeau

Julien Tranchida

Stam Nicolis

Non-Markovian magnetization dynamics for uniaxial nanomagnets

A stochastic approach for the description of the time evolution of the magnetization of nanomagnets is proposed, that interpolates between the Landau-Lifshitz-Gilbert and the Landau-Lifshitz-Bloch approximations, by varying the strength of the noise. Its finite autocorrelation time, i.e. when it may be described as colored, rather than white, is, also, taken into account and the consequences, on the scale of the response of the magnetization are investigated. It is shown that the hierarchy for the moments of the magnetization can be closed, by introducing a suitable truncation scheme, whose validity is tested by direct numerical solution of the moment equations and compared to the averages obtained from a numerical solution of the corresponding colored stochastic Langevin equation. This comparison is performed on magnetic systems subject to both an external uniform magnetic field and an internal one-site uniaxial anisotropy.

I. INTRODUCTION

Idealizations are often employed, when modeling real magnetization dynamics. Far from equilibrium, magnetic fluctuations, modeled by white noise are examples of such an idealization and the magnetization processes predicted in this way have been observed in ultrafast reversal magnetization experiments. However real fluctuations always have finite auto-correlation time and the corresponding spin system dynamics is, rather, described by non-Markovian stochastic processes [START_REF] Atxitia | Ultrafast spin dynamics: the effect of colored noise[END_REF] . Despite the practical relevance of non-Markovian processes, however, these have received much less attention than Markovian processes, chiefly, because of the simplicity of the theory and the familiarity of the corresponding mathematical tools [START_REF] Gardiner | Handbook of Stochastic Methods in Physics[END_REF] . The study of the magnetization dynamics of super-paramagnetic systems, however, is a suitable framework and provides motivation to describe magnetic fluctuations, induced by contact to a heat bath 3 , that are sensitive to non-Markovian effects.

The effective magnetization dynamics of a spin system and a noise is described by a stochastic differential equation of Langevin type, also known as the stochastic Landau-Lifshitz-Gilbert (sLLG) equation :

(1 + λ 2 ) ∂s i ∂t = ǫ ijk (ω j + Ω j )s k + λ(ω i s j s j -ω j s j s i ). (1)
The Einstein summation convention is adopted, with Latin indices standing for vector components, ǫ ijk is the anti-symmetric Levi-Civita pseudo-tensor and s is a magnetization 3-vector, whose norm is fixed and can be set to unity. Here ω is the precession frequency which derives from a spin Hamiltonian and is such connected to the a) Electronic mail: pascal.thibaudeau@cea.fr b) Electronic mail: julien.tranchida@cea.fr c) Electronic mail: stam.nicolis@lmpt.univ-tours.fr spin vector. λ is the damping coefficient. The stochastic properties of Ω describe the fluctuations of external and/or internal degrees of freedom, random perturbations and noise. Its correlation functions are assumed to satisfy the relations

∂ ∂t Ω i (t) . . . Ω in (t n ) = - 1 τ Ω i (t) . . . Ω in (t n ) (2)
with and t ≥ t 1 ≥ . . . ≥ t n is the time-ordered statistical average product of stochastic variables 4 and τ is a single characteristic constant time. With proper conditions Gaussian process, Poisson process and two-state (dichotomy) Markov process enjoy this "time derivative property" [START_REF] Shapiro | Formulae of differentiation and their use for solving stochastic equations[END_REF] . In general a detailed description of the stochastic processes in the extended phase space (s, Ω) is not easy to establish [START_REF] Miyazaki | Brownian motion of spins revisited[END_REF][START_REF] Thibaudeau | Thermostatting the atomic spin dynamics from controlled demons[END_REF] .

Stochastic processes are completely described by an infinite set of multivariate probability distributions. In practice, not all multivariate distributions are not desirable but only a few, mainly the one-variable distribution P (s, t) and the two-variable distribution, P (s 1 , t 1 ; s 2 , t 2 ). The former is used to determine the average f (s, t) of any functional f of the process s(t), e.g. its statistical moments. The latter is used to determine two-time statistics, e.g., the correlation function s i (t 1 )s j (t 2 ) . In the white noise limit, a Fokker-Planck equation on P (s, t) has been established only for a thermally interacting spin system subject to a constant external field [START_REF] García-Palacios | Langevin-dynamics study of the dynamical properties of small magnetic particles[END_REF][START_REF] Brown | Thermal fluctuations of a single-domain particle[END_REF] . Beyond the mean-field approximation, the explicit spin dependance of the effective field, including a magnetic anisotropy, on the probability distribution P , may be addressed.

II. THEORY

In contact with a thermostat, described by the stochastic vector, Ω(t), the spin system precesses about the axis of magnetic effective field, defined by both a uniform external magnetic field with ω (0) = γµ 0 H and an anisotropic single-axis. This can be observed in superparamagnetic nano-magnets which behave as effective, two-state, anisotropic spin systems 3 . The deterministic frequency vector is given by ω = ω (0) + κn (n • s), where n is the anisotropy single-axis and κ the anisotropy frequency. The noise vector, Ω(t), is assumed to be described by a Gaussian stochastic process, that is completely defined by its 1-point function Ω(t) and 2-point function Ω i (t)Ω j (t ′ ) . Therefore the only non-trivial information is contained in the these, that we assume to define an Ornstein-Uhlenbeck process 9 , which means that Ω i (t) = 0 and

Ω i (t)Ω j (t ′ ) = D τ exp(- |t -t ′ | τ )δ ij . (3)
It can be shown that this process does converge to the centered white noise process, when τ → 0. Moreover, it is assumed that the thermostat can be described by a fluctuation-dissipation theorem that relates the strength of the noise, D, with the temperature T by the relationship D = λk b T /h, where λ is the parameter that appeared in Eq.( 1).

These microscopic degrees of freedom give rise to average quantities, that probe collective properties. Such quantities can be operationally defined by taking statistical averages over the noise of Eq.( 1) at equal times. In the limit of low damping (i.e. λ ≪ 1), the equation for the magnetization, s i , is :

∂ s i ∂t = ǫ ijk ( ω j s k + Ω j s k ) +λ( ω i s j s j -ω j s j s i ), (4)
and implicates higher order correlation functions of noise and spin that have to be determined in turn. By applying the Shapiro-Loginov formulae of differentiation [START_REF] Shapiro | Formulae of differentiation and their use for solving stochastic equations[END_REF] to the components of the Ω i s j at equal times, we obtain the equation

∂ Ω i s j ∂t = Ω i ∂s j ∂t - 1 τ Ω i s j . (5)
Introducing the right-hand side of Eq.(1) into Eq.( 5), evolution equations for the mixed moments Ω i s j are obtained. Due to the non-linearity of the sLLG equation, an infinite hierarchy arises 10 : the equations for the second-order moments depend on third-order moments, and so on. In order to solve the system, closure relationships need to be found. So the equations for the moments s i s j are now required and are provided by ∂ s i s j /∂t = s i ṡj + ṡi s j . Assuming that this identity holds, the time derivative of s i is replaced by the right-hand side of Eq.( 1) inside the statistical average. This leads to an additional set of equations for the second-order moments of the microscopic spin degrees of freedom, that are closely related to those written by Garanin et al (see Eq.( 6) [START_REF] Garanin | Dynamics of an ensemble of single-domain magnetic particles[END_REF] ) and Garcia-Palacios et al.

(see Eq.(2.10) 3 ). As long as the stochastic variables follow individual Gaussian processes, their connected moments (labelled as . ), beyond the second moment are zero [START_REF] Van Kampen | Stochastic processes in physics and chemistry[END_REF] . The case of the mixed connected moments of spin and noise is more subtle. If a linear coupling between the spin and the noise is assumed, the mixed cumulants, above the second will, also, vanish. If the Gaussian hypothesis is dropped, the general relations, at equal times, between the desired moments are :

Ω i s j s k = Ω i s j s k + Ω i s k s j + Ω i s j s k , (6) Ω i s j s k s l = Ω i s j s k s l + Ω i s j s k s l + Ω i s k s j s l + Ω i s l s j s k + s l Ω i s j s k + s k Ω i s j s l + s j Ω i s k s l , (7)
s i s j s k = s i s j s k + s j s k s i + s i s k s j -2 s i s j s k + s i s j s k , (8)
s i s j s k s l = s i s j s k s l + s i s k s j s l + s i s l s j s k -2 s i s j s k s l + s k s i s j s l + s j s i s k s l + s l s i s j s k + s i s j s k s l + s i s j s k s l . (9)
These considerations lead to the following system of equations :

∂ s i ∂t = ǫ ijk (ω (0) j s k + κn j n l s l s k + Ω j s k ) +λ ω (0) i s j s j -ω (0) j s j s i +κn l (n i s l s j s j -n j s l s j s i ) , (10)
∂ Ω j s k ∂t = - 1 τ Ω j s k + ǫ klm ( Ω j ω l s m + Ω j Ω l s m ) , +λ ( Ω j ω k s l s l -Ω j ω l s l s k ) , (11)
∂ s i s j ∂t = ǫ ikl ( ω k s l s j + Ω k s l s j ) +λ ( ω i s k s k s j -ω k s k s i s j ) +(i ↔ j). (12)
These expressions are closed by considering first the independence of Gaussian processes, the form of the effective frequency vector and the expression of the autocorrelation function of the noise. One has explicitly : tion is explicitly treated at the third-order.

Ω j ω l s m = ω (0) l Ω j s m +κn l n u ( Ω j s u s m + Ω j s m s u ) Ω j Ω l s m = D τ δ jl s m Ω j ω k s l s l = 2ω (0) k Ω j s l s l +κn k n u ( Ω j s u s l s l + 2 Ω j s l s u s l ) Ω j ω l s l s k = ω (0) l ( Ω j s l s k + Ω j s k s l ) +κn l n u ( Ω j s u s l s k + Ω j s l s u s k + Ω j s k s u s l ) ω k s l s j = ω (0) k s l s j + κn k n u ( s u s l s j + s l s u s j + s j s u s l -2 s u s l s j ) Ω k s l s j = Ω k s l s j + Ω k s j s l ω i s k s k s j = ω (0) i ( s j s k s k + 2 s k s k s j
The third-order correlation functions for the dynamical variables s satisfying a stochastic Langevin-like equation are generally determined by first solving this equation in the presence of the noise and then averaging the product of these variables over the ensemble of the random noise. A path integral formalism can be established for deducing identities between the correlation functions. This strategy has been tested in various approximation schemes, including the quenched approximation (i.e. D = 0). A shortcut might be to look for an alternate description of an effective system of equations in a closed time path formalism of Schwinger allowing to describe a general nonequilibrium field theory at finite temperature [START_REF] Kleinert | Path integrals in quantum mechanics, statistics, polymer physics, and financial Markets[END_REF] . This is why a closure approximation was proposed by Martin, Siggia and Rose 17 which represents the classical limit of the closed time path formalism and will be the subject of further investigation.

JT acknowledges financial support through a joint doctoral fellowship "Région Centre-CEA" under grant number 00086667.

Figure 2 .

 2 Figure 2. Stochastic average magnetization (left) and autocorrelation spin function (right) around a single-axis in the z-direction and a field in the x-direction. The dot lines are for the sLLG for 10 4 repetitions and the solid lines stand for our dLLB model. D = 1 × 10 -3 GHz; λ = 0.05; ω (0) = (π/5, 0, 0)GHz ; κ = 1/10GHz; n = (0, 0, 1); τ = 10 -1 ns are taken with the same initial conditions as previously.

Figure 3 .

 3 Figure 3. Stochastic average over for 10 5 repetitions of the independent sisjs k terms for various values of D. λ = 0.1; ω (0) = (π/5, 0, 0)GHz ; κ = 0.5GHz; n = (0, 0, 1); τ = 10 -3 ns are taken with the same initial conditions as previously.

-2 s k s k s j ) + κn i n u ( s k s k s u s j + 2 s u s k s k s j -2 s u s k 2 s j ω k s k s i s j = ω (0) k ( s j s k s i + s k s i s j + s k s j s i -2 s k s i s j ) +κn k n u ( s u s k s i s j + s u s i s k s j + s u s j s k s i -2 s u s k s i s j )

Combined with Eqs. (10,11,12), the derived system of equations is called the dynamical Landau-Lifshitz-Bloch model (dLLB). In the limit of τ → 0, Eq. (11) gives a direct expression for the mixed moment of the noise and the magnetization as Ω i s j = -Dǫ ijk s k , which makes the Ω i s j matrix antisymmetric. Once introduced into Eq.( 10), a longitudinal contribution for the spin dynamics appears in which the amplitude of the noise is proportional to the square root of the temperature. For finite values of τ , Eq.( 10) exhibits a "memory effect" which shows that the value of s i (t) for the time t depends on all the previous values of s i (t ′ ) before the time t.

III. NUMERICAL EXPERIMENTS

Numerical experiments were performed to test the consistency of the closure hypothesis. For "low" values of τ , the dynamics is close to Markovian and Eqs. (10,11,12) are expected to match the results obtained from an averaged Markovian stochastic dynamics. In order to realize these comparisons, a reference model is constructed from Eq.( 1) with a white noise random vector Ω. Multiple spin trajectories with the same initial conditions are generated and followed in time according to a previously detailed symplectic integration scheme [START_REF] Thibaudeau | Thermostatting the atomic spin dynamics from controlled demons[END_REF][START_REF] Lewis | Geometric integration on spheres and some interesting applications[END_REF][START_REF] Beaujouan | Anisotropic magnetic molecular dynamics of cobalt nanowires[END_REF] . A statistical average is taken over all the realizations. Eqs. (10,11,12) are integrated with an eight-order Runge-Kutta integration scheme, with adaptive steps, constructed from the Jacobian of the whole system.

The figure (1-left) (resp. 1-right) exhibits the average magnetization (resp. spin-spin autocorrelation function) from the stochastic LLG equation and our dLLB model for the same value of the amplitude D in the white noise limit. For low value of D (i.e. in the normal phase), the average stochastic dynamics follows our dLLB model. A similar situation has been already investigated for an external constant field only [START_REF] Tranchida | Closing the hierarchy for non-markovian magnetization dynamics[END_REF] . When τ increases, the noise is no longer white and does not follow a Markov's process. A complete stochastic system of equations for both the spin and the noise is built by supplementing Eq.( 1) with

where ξ i is a Gaussian centered stochastic process [START_REF] Fox | Fast, accurate algorithm for numerical simulation of exponentially correlated colored noise[END_REF] , but with ξ i (t) = 0 and ξ i (t) ξ j (t ′ ) = 2Dδ ij δ(t -t ′ ). Formally, this system is integrated as Ω i (t) ≡ e tLΩ Ω i (0) = t 0 e -(t-u)/τ ξ i (u)du with L Ω is the Liouville operator corresponding to Eq.( 13). Because the Liouville operators of the spin and the noise commute, an integration of the colored system is performed by the combination of the two operators as e tL (s i (0), Ω i (0)) = e tLs e tLΩ (s i (0), Ω i (0)) = (s i (t), Ω i (t)). This is performed for larger values of τ for which the figure (2) exhibits s i (t) and the diagonal part of s i (t)s j (t) . It is observed that, even for "large" values of τ , our dLLB model continues to track the dynamics of the colored stochastic LLG system. The influence of the higher than second-order cumulants on the dynamics of the average magnetization and spin autocorrelation function is mitigated by a large damping factor λ which prevents detailed descriptions of the magnetization dynamics. Moreover, because of the link between the amplitude D and the damping factor λ, through a tobe-derived fluctuation-dissipation theorem in the special case of an explicit dependance of the spin in the effective field, the domain of validity of our dLLB model in temperature could be significantly controlled.

When the intensity of the noise increases, the spin dynamics is no longer Gaussian. It is observed numerically, in the stochastic dynamics, by calculating the cumulant function s i (t)s j (t)s k (t) . The ten independent spin cumulants are depicted, for several values of the amplitude D in the figure (3). As shown, for this particular situation, some terms go to zero in time, whereas some remain finite, depending on the underlying symmetries. At the expense of higher order expressions, the Eqs. (10,11,12) can be supplemented by an expansion of d sisj s k dt and integrated together. This represents an option where the deviation off the Gaussian character of the spin distribu-