Dariusz Kajtoch

Krzysztof Pawłowski

Emilia Witkowska

Entanglement storage by classical fixed points in the two-axis countertwisting model

Keywords: numbers: 03.67.Bg, 03.75.Dg, 03.75.Gg

We analyze a scheme for storage of entanglement quantified by the quantum Fisher information in the two-axis countertwisting model. A characteristic feature of the two-axis countertwisting Hamiltonian is the existence of the four stable center and two unstable saddle fixed points in the mean-field phase portrait. The entangled state is generated dynamically from an initial spin coherent state located around an unstable saddle fixed point. At an optimal moment of time the state is shifted to a position around stable center fixed points by a single rotation, where its dynamics and properties are approximately frozen. We also discuss evolution with noise. In some cases the effect of noise turns out to be relatively weak, which is explained by parity conservation.

I. INTRODUCTION

The idea of dynamical generation of entangled states in ultra-cold atomic systems has its origin in the work of Kitagawa and Ueda [1]. The concept of spin-squeezed states was introduced in a collection of qubits. Two different models for dynamical generation of spin-squeezed states were proposed, namely the one-axis twisting (OAT) and the two-axis countertwisting (TACT) Hamiltonians, with the last one giving the strongest level of squeezing. Later on, it was recognized that highly entangled states, including the NOON state, can also be generated by further dynamics. However, due to small nonlinearities and decoherence only atomic states appearing at the beginning of the OAT scenario have been observed [2][3][4][5][6][7] so far, and used to improve precision in prototypes of atomic clocks [7] and magnetometers [8]. On the contrary, the TACT model, being the subject of this paper, has not been realized experimentally yet. Even though, it still attracts attention of many physicists, as a high degree of squeezing and entanglement can be reached at much shorter time scales (or alternatively smaller nonlinearities) [9].

Apart from difficulties in realizing twisting Hamiltonians experimentally, there are other problems related to metrological schemes. Once the desired entangled state is reached, the problem which appears is how to save it in a way robust against decoherence. In many systems, including ultracold atoms, it is not easy to simply switch off the nonlinearity to avoid further evolution into possibly less interesting states. In this paper we propose a scheme, based on the TACT model, in which the entangled state is reached and then its entanglement is preserved for an infinite time even though the nonlinearity is not reduced, and even more surprisingly, also some noise is included.

The main idea of the scheme is very simple, and utilizes the structure of the mean-field phase portrait. In our recent paper [9] we have shown that the mean-field phase portrait of the TACT Hamiltonian consists of two unstable saddle and four stable center fixed points located symmetrically on the Bloch sphere. The strong squeez-ing and entanglement can be reached from an initial spin coherent state located around an unstable saddle fixed point which is reflected in the Heisenberg-like scaling of the quantum Fisher information. For the initial spin coherent state located around a stable center fixed point the quantum dynamics is approximately frozen. The scheme we have in mind joins the qualitatively different dynamics which takes place in a vicinity of a fixed point. The initial spin coherent state located around an unstable fixed point is evolved via the TACT Hamiltonian to an interesting non-Gaussian state. Then, the state is rotated in order to locate its important parts around stable fixed points. As a result, further dynamics is confined to a narrow region of the phase space around two antipodal stable fixed points, and the quantum Fisher information will forever remain at a very high level.

A natural question which arises is how noise, always present in experiments, affects the final results of the scheme? We address this question by investigating the dynamics of the quantum Fisher information with noise inspired by experiments with ultracold two-level atoms. The serious problem here is the noise due to external slowly changing fields, which affects the energy levels and hence causes extra rotation of the state. These changes are assumed to be such slow that during a single experimental run the fields are constant, but they vary from shot to shot. Hence, from one experimental realization to another one obtains an entangled state landing at a different position. The state smears out and the entanglement is reduced while averaging over realizations. The experimental way around is a so-called spin-echo technique. We incorporate such noise by adding stochastic terms into the Hamiltonian. Instead of a single pure state we have a collection of them for each realization of the noise. The final state is obtained as a mixture of them. Such the procedure, at least in the OAT scheme, is also related to other sources of decoherence, i.e. particle losses and finite temperature effects [10]. The effect of noise on the proposed scheme turns out to be relatively weak since one ends up with the Heisenberg-like scaling of the quantum Fisher information. We observe that, depending on details of the noise, the quantum Fisher information ei-ther decreases down to some universal value or remains at a very high level, almost unaffected. We explain the latter finding on the quantum level by showing that it is protected by the conservation laws of energy and parity.

II. THE MODEL

We consider a collection of N qubits e.g. particles in two orthogonal modes. The system is conveniently described using the collective spin operator ˆ S whose components written in the Schwinger representation are

Ŝx = 1 2 â † b + b † â , (1a)
Ŝy = 1 2i â † b -b † â , (1b)
Ŝz = 1 2 â † â -b † b , (1c)
where â † , b † are the bosonic creation operators associated with the two modes.

A representative form of the two-axis countertwisting Hamiltonian was proposed in [1],

Ĥ = χ Ŝx Ŝy + Ŝy Ŝx . (2)
An SU(2) rotation of the Hamiltonian Û Ĥ Û † , where Û is a group element, produces a mathematically equivalent form which appears e.g. in the Lipkin-Meshkov-Glick model [11]. In this paper we will operate on the rotated Hamiltonian ĤTACT = -χ Ŝy Ŝz + Ŝz Ŝy ,

with Û = e -i Ŝyπ/2 , as a matter of convenience simplifying the form of observables of interest.

In what follows, we will consider the effect of noise on the quantum dynamics. A source of the noise is closely related to experimental conditions. In applications to Bose-Einstein condensates the noise was effectively modeled as a linear combination of spin components added to an unperturbed Hamiltonian [4,10,12]. To adopt the same approach we first introduce the family of Hamiltonians

Ĥ γ = ĤTACT + γ • ˆ S, (4)
where γ = (γ x , γ y , γ z ). In this convention the unperturbed Hamiltonian, denoted by Ĥ 0 , is equal ĤTACT . We assume that the parameters γ j are time independent but vary randomly from one experimental realization to another mimicking a stationary random dephasing environment. We assume that they are normally distributed with the Gaussian probability density

P (γ j ) = 1 σ j √ 2π e - γ 2 j 2σ 2 j , (5)
where σ j is the standard deviation.

In the presence of noise a quantum state of the system is described by a density matrix operator ρ(t) which is constructed as a statistical average over stochastic realizations,

ρ(t) = d γ j P (γ j ) |ψ γ (t) ψ γ (t)| . (6)
Here, the pure states |ψ γ are the solutions of the Schrödinger equation

i ∂ t |ψ γ (t) = Ĥ γ |ψ γ (t) . (7)
The initial state for the evolution is the spin coherent state located along the X axis of the Bloch sphere with radius N/2. In a standard way one can expand it in the Fock state basis [9],

|ψ γ (0) = 2 -N/2 N k=0 N k 1/2 |k, N -k , (8)
which is the eigenstate of the Ŝx operator with the eigenvalue N/2. The spin coherent state can be visualized on the Bloch sphere as a disk of diameter √ N /2. In order to quantify the amount of quantum correlations that are useful for atomic interferometry we will concentrate on the quantum Fisher information (QFI). According to the Cramér-Rao bound [13], the precision of a two-mode linear interferometer is bounded by the inverse square root of the QFI [14]. The higher the QFI the better the interferometric precision. It was also recognized that the QFI is a measure of multiparticle entanglement [15,16]. The entanglement is necessary to beat the shot-noise limit with F Q = N which is characteristic for uncorrelated particles [14,[17][18][19][20]. The upper bound for precise measurements with F Q = N 2 is attainable by maximally entangled states, e.g. the NOON state [21]. The value of the QFI is

F Q = 4λ max , (9)
where λ max is the maximal eigenvalue of the covariance matrix [15][16][17] Γ

ij [ρ] = 1 2 l,m (p l -p m ) 2 p l + p m Re l| Ŝi |m m| Ŝj |l , (10)
where p l and |l are eigenvalues and the corresponding eigenvectors of the system density matrix ρ respectively. For a pure state the covariance matrix (10) reduces to

Γ ij = 1 2 Ŝi Ŝj + Ŝj Ŝi -Ŝi Ŝj , (11)
and the QFI is determined by its maximal eigenvalue.

Alternatively, the quantum Fisher information for a pure state can be calculated from

F Q = 4 max n ∆ 2 Ŝ n , (12)
1.0 0.5 0.0 where Ŝ n = n x Ŝx + n y Ŝy + n z Ŝz and the maximization is over all vectors (n x , n y , n z ) satisfying the constraint n 2

-0.5 -1.0 z -π -π/2 0 π/2 π ϕ -π -π/2 0 π/2 π ϕ -π -π/2 0 π/2 π ϕ -π -π/2 0 π/2 π ϕ -π -π/2 0 π/2 π ϕ (a) (b) (c) (e) (d)
x + n 2 y + n 2 z = 1. Time evolution of a quantum system can be traced using a quasi-distribution function in the phase space. It was shown, that the classical Liouvillian flow highly resembles the Husimi function evolution (as long as the curvature of the phase space does not play a significant role) [22]. Let us first analyze the classical phase portrait in order to understand the quantum dynamics and the idea for entanglement storage.

III. MEAN-FIELD PHASE PORTRAITS

In the mean-field description (the limit of large system size N ≫ 1) one replaces bosonic creation and annihilation operators by c-numbers [23],

â → √ N √ ρ a e iϕa , b → √ N √ ρ b e iϕ b , (13)
were ρ a + ρ b = 1 due to conservation of the total particle number. Two canonical variables, namely the population difference z = ρ a -ρ b and relative phase ϕ = ϕ b -ϕ a are sufficient to describe classical dynamics [24]. Conventionally, z ∈ [-1, 1] and ϕ ∈ [-π, π[, reflecting the spherical topology of the phase space. Mean-field counterparts of spin operators take the form

S x = N 2 1 -z 2 cos ϕ, (14a)
S y = N 2 1 -z 2 sin ϕ, (14b)
S z = N 2 z, (14c)
and the quantum Hamiltonian (4) transforms into its classical counterpart

H = -χ N 2 2 z 1 -z 2 sin ϕ + N 2 zγ z + N 2 1 -z 2 (γ x cos ϕ + γ y sin ϕ) . (15)
In principle, one can follow time evolution of the classical system by solving the Hamilton's equations φ =

(2/ N )(∂H/∂z), ż = -(2/ N )(∂H/∂ϕ): φ N χ = - 1 -2z 2 √ 1 -z 2 sin ϕ - z √ 1 -z 2 (γ y sin ϕ + γx cos ϕ) + γz , (16a)
ż N χ =z 1 -z 2 cos ϕ -1 -z 2 (γ y cos ϕ -γx sin ϕ) , (16b)
with γj = γ j /N χ. The dynamics can be qualitatively deduced by looking at the phase portrait which consists of trajectories in the phase space tangent to the velocity field ( φ, ż). In Fig. 1 we show phase portraits for typical cases. A characteristic feature of the TACT Hamiltonian is the existence of unstable saddle and stable center fixed points [9]. When γ = 0, the phase portrait has four stable center fixed points located at (z, ϕ) = (±1/ √ 2, ±π/2) and two unstable saddle fixed points at (z, ϕ) = (0, 0) and (z, ϕ) = (0, -π) (compare with Fig. 1a). Nonzero value of γj shifts locations of fixed points, and a bifurcation occurs when |γ j | = 1 (or γ j = N ). A striking feature is the stability of the phase portrait against weak perturbations (γ j < N ) i.e. if a fixed point was of unstable saddle type, it will remain unstable saddle (the same applies for stable centers). For |γ j | > 1 the mean-field dynamics enters a Rabi-like regime with two stable fixed points remaining in the phase portrait.

Explicit analytical expressions for the positions of particular fixed points are given in Appendix A.

IV. STOCHASTIC EVOLUTION OF THE QUANTUM FISHER INFORMATION

Time evolution of the QFI in the unperturbed case γ = 0 was analyzed in [9] (shown also by the red solid line in Fig. 2). We briefly remind its regular behavior for short times using the Husimi function in the phase space picture. The initial coherent state located around the unstable saddle fixed point (visualized in Fig. 1a by the red disk) is stretched along the prime meridian towards the second unstable fixed point located at the opposite side of the Bloch sphere. The QFI grows up from F Q = N , and reaches the first maximum at χt = ln(2πN )/2N with the Husimi density cumulated around the two poles of the Bloch sphere (see Fig. 3a). At subsequent moments of time, the QFI decreases reaching a minimum at χt = ln(2πN )/N when separated pieces of the Husimi distribution meet together at the opposite side of the Bloch sphere around the second unstable fixed point.

The stretching and separation of the Husimi function takes place once again, but this time along the equator towards the initial fixed point. The QFI grows up reaching the second maximum at χt ≃ 1.75 ln(2πN )/N . Thereafter, the dynamics becomes irregular and overall results depend on the total number of particles.

In the presence of noise the QFI is determined by the system density matrix according to Eq. (10). A numerical procedure for obtaining ρ(t) is based on a discrete limit of Eq. ( 6):

ρ(t) = 1 M M i=1 ρi (t), (17)
where ρi (t) = |ψ i (t) ψ i (t)| and |ψ i (t) is a solution of the Schrödinger equation (7) for given γ i . The optimal ensemble size M can be determined e.g. by looking at the von-Neumann entropy: when its value stabilizes as a function of M then the structure of the density matrix does not change.

Numerical results for the QFI with nonzero noise term along a chosen axis of the Bloch sphere are shown in Fig. 2. Illustrative analysis of the QFI using the phase space picture is hindered in the mixed state case. Nevertheless, the information extracted from each realization of γ i is valuable to the overall behavior of the QFI in a noisy system. Due to convexity of the QFI [14,25,26], a convex mixture of quantum states contains fewer quantum correlations than the ensemble average

F Q [ρ(t)] < 1 M M i=1 F Q [|ψ i (t) ] . (18)
If the QFIs for particular realizations of γ i lie below the unperturbed value then the QFI of the noisy system decreases. This is what we observe in the early stage of the evolution for χt < ln(2πN )/2N . On the other hand, if the QFIs for particular realizations of γ i exceed the value for the unperturbed case in some time interval, then the inequality (18) does not exclude the possibility that the QFI value may be larger as compared to the noiseless case. This effect is observed for noise along the X axis and time χt around ln(2πN )/N (see Fig. 2c).

As a general principle, the stronger the noise the greater the reduction of the QFI value. Suppression of the QFI is observed for nonzero γ z (Fig. 2a) and nonzero γ y (Fig. 2b). When |γ y, z | ≪ χ √ N then the initial state remains within the range of the unstable fixed point, and the evolution of the QFI is close to the noiseless case. When |γ y, z | is of the order of or larger than χ √ N then the position of the unstable fixed point shifts, and effectively the initial state is no more within its range, see Figs. 1d and1e. This prevents the first spreading and separation which leads to a decreased value of the QFI. A different situation takes place for a nonzero γ x and γ y, z = 0. The positions of unstable fixed points do not change, and the time evolution is always initialized at the unstable fixed point, see Fig. 1b. As a result, the first part of the evolution is very close to the noiseless case. The impact of the linear term enhances at later times due to modification of the classical trajectories, see Fig. 1b. The noise-enhanced effect is spotted around χt = ln(2πN )/N , as it can be observed in Fig. 2c. In order to give ground for comparison we show the QFI at χt = ln(2πN )/N for the noisy system (left panel) and a particular realization of noise (right panel) in insets of Fig. 2c. A sufficiently large γ x opens a new path for the evolution. The separated pieces of the distribution function have no occasion to localize around the second unstable fixed point, they turn around the stable fixed points and meet earlier at the initial saddle fixed point. Consequently, the state remains in the two elongated pieces for a longer time increasing the value of the QFI around χt ≃ ln(2πN )/N . When γ j > N unstable fixed points disappear in the phase portrait and the QFI is never larger than N in the limit γ j → ∞.

V. THE STORAGE SCHEME

The regular part of the evolution and exceptional behavior of the Husimi function give a possibility of a stabilization scheme with nearly stationary value of the QFI at a relatively high level. The scheme consists of three steps. The quantum system is evolved until the QFI reaches the first maximum at τ = ln(2πN )/2N χ. Then an instantaneous pulse, e.g. of laser light, rotates the state through π/4 around the X axis,

ψ 0 (τ + ) = e i Ŝxπ/4 ψ 0 (τ -) , (19)
where τ + and τ -denote the time just after and before the rotation respectively. Later on (t ≥ τ + ), the dynamics is governed by the TACT Hamiltonian without any manipulations. The idea of the scheme can be understood as illustrated in Fig. 3a. As we know from previous considerations, shortly before the pulse the Husimi function is highly stretched with the two maxima localized around the two poles of the Bloch sphere. Rotation of the state throws the two maxima into stable regions A small rotation of the state results in a small shifts of the distribution of quasiprobability, here the SU(2) Wigner function [27]. Thus the appearance of small thin structures in the Wigner function leads to a high value of the QFI, as rotating the state perpendicularly to them would strongly affect the state and the quasiprobability distribution. The high narrow fringes, so also the high value of the QFI, remain even in the case with noise along the X axis, but they are mostly smeared out by noise along Y or Z axis.

of the phase space. Thereafter, the dynamics is trapped around the two antipodal stable fixed points. An animation for time evolution of the Husimi function is in [28]. Time evolution of the QFI is shown in Fig. 3b. A roughly stationary value was obtained owing to the rotation. This result does not depend much on the number of particles, and the average value of the QFI remains at the level F Q ≃ 0.75N 2 . We emphasize that a deviation from the optimal time τ or angle of rotation up to 20% does not spoil the "freezing" scheme, but rather lowers the value of the QFI, at worst to

F Q ∼ 0.7N 2 .
In fact, within the scheme not the state itself is saved but rather its useful entanglement quantified by the QFI. Before the rotation the QFI given by Eq. ( 12) is maximized for n = (0, 0, 1), which leads to

F Q (τ -) = 4∆ 2
Ŝz (τ -). The optimal direction after the pulse is

n = 0, 1/ √ 2, 1/ √ 2 .
For t ≥ τ + , the QFI given by the definition (12) fulfills inequality

F Q (t) ≥ 4 Ŝy Ŝz + Ŝz Ŝy (t),
where we have used the fact that Ŝy (t) = Ŝz (t) = 0 (from the symmetry arguments) and the variance ∆ 2

Ŝy -Ŝz (t) is nonnegative. The operator Ŝy Ŝz + Ŝz Ŝy is equal to -H TACT , so it is a constant of motion. Conservation of energy implies Ŝy Ŝz + Ŝz Ŝy (t) = Ŝy Ŝz + Ŝz Ŝy (τ + ). On the other hand, 4 Ŝy Ŝz + Ŝz Ŝy (τ + ) =

Ŝ2

z (τ -) -4 Ŝ2 y (τ -). Just before the pulse, the state is still squeezed along the Y axis, hence the variance ∆ 2

Ŝy (τ -) is relatively small. Finally we have

F Q (t) ≥ F Q (τ -) -4∆ 2 Ŝy (τ -) ≈ F Q (τ -). (20)
It follows that after the rotation the QFI needs to stay at least at the level from the moment of pulse.

A. The effect of noise

In the following we analyze the effect of weak noise on the scheme. We limit the regime of parameters in such a way that most of γ i 's in an ensemble leave the initial state within the range of the unstable fixed point. We also demand that the fidelity function between the unperturbed state and the density matrix coming from the stochastic dynamics

F σj (t) = ψ 0 (t) ρσj (t) ψ 0 (t) = ∞ -∞ dγ j P (γ j ) ψ 0 (t) ψ γj (t) 2 ( 21
)
at t = τ will stay close to 1. The latter defines weak noise. The conditions for the σ i parameters are the following:

(i) σ y, z /χ ≪ √ N /2 and σ x /χ < N, (22a)
(ii) σ x,y,z /χ ≪ 2 √ 2N / ln(2πN ). (22b)
For large N the second condition is stronger than the first one. For example, N = 50 requires σ x,y,z /χ ≪ 3.5. Time evolution of the QFI with the rotation at t = τ for the system with noise is shown in Fig. 3c. Two qualitatively different cases can be distinguished. For γ x = 0 and γ y, z = 0, stabilization of the QFI around F Q x /N 2 ≃ 0.75 is observed, where . x denotes here average in time. The overall shape of the curve stays very close to the result for the unperturbed system even for a very long period of time. This behavior does not depend on the parity of N . In the second case, for γ y = 0 and γ x, z = 0 (or γ z = 0 and γ x, y = 0), initial decaying followed by stabilization of the QFI is observed. The decay time depends both on N and σ y, z , while the average level of the QFI stabilization is typically ≈ 0.2N 2 , irrespective of σ y, z . We noticed that the stable value of the QFI oscillates in the same way even for different values of σ y,z .

These observations agree with the geometrical consideration based on the Wigner function, as illustrated and discussed in Fig. 4.

B. Noise along the Z or Y axis

A distinguishing feature of the QFI is the decaying character observed for γ z = 0 and γ x, y = 0 (or γ y = 0 and γ x, z = 0). The half-decay time behaves as τ 1/2 ∝ 1/(N σ y, z ) according to our numerical simulations, see Fig. 5a. Understanding of the decay process goes handin-hand with the eigensystem of the density matrix since formula (10) for the QFI with noise requires its knowledge. As shown in the inset of Fig. 5a by the red solid line, during the initial decay of the QFI the spectrum is dominated by a single eigenvalue p l ≈ 1. For such spectrum the covariance matrix (10) can be simplified to

Γ ij ≃ p l 2 Ŝi Ŝj + Ŝj Ŝi -2 Ŝi Ŝj , (23)
where the quantum average is calculated over the eigenstate |l corresponding to p l . Moreover, we observe that the dominant eigenvalue p l coincides with the fidelity (21), marked by the black dashed line in the inset of Fig. 5a. This means that the eigenstate |l is nothing else but just the unperturbed state ψ 0 (t) . If this is the case, then the covariance matrix ( 23) is equal to the covariance matrix for the unperturbed evolution multiplied by the decaying weight p l . The initial decay of the QFI in the case with noise is then determined by the p l , being equal to the fidelity (21). In order to compute the decay rate of the QFI we move on to the analysis of the fidelity function. The state after the rotation is

|ψ γ (t) r = e -i Ĥ γ T e i Ŝxπ/4 ψ γ (τ -) , (24)
where T = t -τ . For noise along the Z axis one has

r ψ γz (t)| ψ 0 (t) r 2 ≃ k, k ′ |c k | 2 |c k ′ | 2 e -iγzT (k-k ′ ) , (25)
where c k are decomposition coefficients in the Fock state basis of the state (24) at T = 0. Formula (25) can be found by splitting the exponential e -i Ĥ γ T into the product e -i Ĥ 0 T e -i γz ŜzT . This can be done for T < T * , where T * ≃ (6) 1/3 ( 3 χγ 2 z N ) -1/3 is determined using the Zassenhause formula and BBGKY hierarchy for expectation values of operator products [9,29]. Integration over γ z gives [30]

F σz (T ) = 2 N n=1 e -T 2 σ 2 z n 2 /2 × N k=n |c k | 2 |c k-n | 2 + N k=0 |c k | 4 .
(26) It is clear that the fidelity function is a monotonically decaying function from

F σz (0) = 1 to F σz (+∞) = k |c k | 4 .
Moreover, it is a dimensionless quantity and therefore a function of an independent dimensionless combination that we can form from σ z , T and N . It is evident from Eq. ( 26) that the half-decay time τ 1/2 should be proportional to inverse of σ z , but the question is about the dependence with respect to N . In order to get some insight to our case we first analyze a simpler case with |c k | 2 = 1/(N + 1). It can be shown that the fidelity function for relatively large N is well described by

F σz (T ) = 1 λ √ π Erf(λ) - 1 λ 1 -e -λ 2 , (27)
with λ = N σ z T / √ 2, see Appendix B for analytical calculations. Our case is harder to analyze because the weights |c k | 2 are distributed in a numerically known way. Nevertheless, we performed numerical calculations for different N and σ z and extracted the same scaling of the fidelity function ( 26) for short times. The same holds for the exact fidelity function defined by Eq. ( 21).

When t ≫ τ , then other eigenvalues of the system density matrix become relevant and the analysis spoils out. One cannot simplify the covariance matrix Γ by (23) any more. In the limit t → ∞, all eigenvalues of the system density matrix stabilize at some nonzero level which leads to a nonzero QFI at long times.

Analogous reasoning can be applied for noise along Y axis, but this time one should decompose the state in the eigenbasis of Ŝy .

C. Noise along the X axis and almost optimal interferometric scheme

The evolution with noise along the X axis turned out to be very peculiar: even for quite strong noise, i.e. σ x ≈ χ, still the QFI only slightly oscillates at the average level F Q x /N 2 = 0.75 (see in Fig. 5b) and it is not damped even for the longest time that we were able to study numerically (up to t = 2000τ ). In order to present our understanding of the robustness of the QFI against noise along the X axis, we have to remind the basics of linear interferometry. In our system linear interferometry can be viewed as rotation around a given axis n through an unknown angle θ, followed by measurement of some observable Ô. The task of interferometry is to estimate θ having outcomes of measurements of Ô. The Cramér-Rao bound states that the uncertainty of an estimation is limited by the inequality ∆θ ≥ [F ( Ô, n)] -1/2 , where F ( Ô, n) is the Fisher information associated with the chosen interferometric direction n and the chosen observable Ô. Possible outcomes of measurements of Ô are the eigenvalues µ. The probability of measuring µ is equal to p(µ) = Tr ρ Pµ , where Pµ is the operator projecting onto a subspace of eigenvectors corresponding to the eigenvalue µ. The conditional probability of measuring µ in the state rotated around n through an angle θ reads p(µ|θ) = Tr e -iθ Ŝ n ρ e iθ Ŝ n Pµ .

For θ = 0, the Fisher information is

F ( Ô, n) ≡ lim θ→0 µ 1 p(µ|θ) dp(µ|θ) dθ 2 . ( 28
)
The QFI we were presenting so far is the largest possible value of F ( Ô, n) [31], namely for any n and Ô it is known that

F Q ≥ F ( Ô, n). (29)
In this subsection, we use the inequality (29) to bound F Q from below, by computing analytically F ( Ô, n) for well chosen n and Ô.

As the rotation axis n we choose, as in the case with unitary dynamics, the direction n = 0,

1 √ 2 , 1 √ 2
after the pulse, and n = (0, 0, 1) before the pulse. The main problem is to find a good observable Ô, such that the statistics of the measurements would be both, sensitive to rotation around n and robust against noise along the X axis. We suggest

Ô := Π := P+ -P-, (30)
where P+ = n=N,N -2,... |n x x n|, P-= n=N -1,N -3,... |n x x n|, and |n x is the eigenstate of Ŝx with the eigenvalue n -N/2. The two eigenvalues of the operator Π are equal to µ ± = ±1. The operator Π can be written as (-1) Ŝx -N/2 , so it is proportional to the parity operator. This operator have been discussed many times in the context of interferometry as its measurements can saturate the Cramér-Rao bound [32][33][34], which have been partially demonstrated in the experiments [35,36]. Here, this choice is dictated by symmetries of the Hamiltonian. The operator Π commutes not only with the Hamiltonian ĤTACT , but also with its versions used to simulate noise along the X axis, i.e. ĤTACT + γ x Ŝx , and with the operator of the pulse e iπ Ŝx/4 . Moreover, it commutes with the density matrix averaged over the distribution of γ x . Furthermore, the equation H (γx,0,0) , Π = 0 combined with the definition of the parity Π = P+ -Pand the decomposition of the identity I = P+ + Pimplies that the projectors P± are also constants of motion. The initial state (8) can be written as |N, 0 x , so it is entirely contained in the "+" subspace, namely P+ |N, 0 x = |N, 0 x . As the projectors are constants of motion, even for the non zero γ x , one has ρ(t) = P+ ρ(t) P+ . This implies that p(µ + |θ = 0) = 1 and p(µ -|θ = 0) = 0 because of the identity P+ P-= 0.

According to the definition (28), to compute Using these observations we simplify the Fisher information to

F ( Π,
F ( Π, n) = lim θ→0 1 p(µ -|θ) dp(µ -|θ) dθ 2 = 2 lim θ→0 d 2 p(µ -|θ) dθ 2 , (31)
where the last equality follows from the d'Hôpital rule for the 0/0 type of expression. The remaining part is the second derivative, which can be written as

lim θ→0 d 2 p(µ -|θ) dθ 2 = -Tr Ŝ n , Ŝ n , ρ P-= 2 Ŝ n P-Ŝ n . (32)
Finally, rewriting Pas I -P+ and using the fact that for our n one has Pµ Ŝ n Pµ = 0, we find that

F ( Π, n) = 2 lim θ→0 d 2 p(µ -|θ) dθ 2 = 4∆ 2 Ŝ n , (33)
which due to Eq. ( 29) gives the final inequality

F Q ≥ 4∆ 2 Ŝ n . (34)
This result, valid in the case with noise along the X direction, has important consequences. First of all in the case without noise we have

F Q = 4∆ 2 Ŝ n = F ( Π, n), which
shows that the parity is the optimal measurement. As shown in Fig. 6, also in the case with noise the measurements of the parity is almost the optimal one, even though we deal with mixed states. After the pulse, in the case with noise ∆ 2 Ŝ n is of the order of the squared extension of the state on the Bloch sphere. As the dynamics is frozen in a proximity of the stable fixed points, and the distance between them is of the order of N , we expect that F Q has to retain the order of magnitude of N 2 . In deriving the final formulas ( 33) and ( 34) we used only arguments based on symmetries, so these results are valid for any distribution of γ x and also for other types of noise along X, for instance with γ x being a stochastic time-dependent function.

VI. CONCLUSIONS

The short time quantum dynamics generated by the TACT Hamiltonian is determined by its mean-field phase portrait. The mean-field phase portrait consists of four stable center and two unstable saddle fixed points located symmetrically on the Bloch sphere. In the evolution of the initial spin coherent state, which is located around an unstable fixed point, states with the Heisenberg-like scaling of the QFI appear. Weak stochastic time-independent noise terms in the Hamiltonian, which mimic stationary random dephasing environment, do not modify the topology of the phase portrait but shift positions of fixed points and change shapes of classical trajectories. The stochastic evolution reduces the value of the QFI in general, but noise-enhanced effects are also possible in some time interval.

The unique configuration of fixed points in the phase portrait of the TACT model allows for the freezing of entanglement quantified by the QFI. In the scheme we have proposed, a high stationary value of the QFI have been achieved by a single rotation of the state which locates it around stable fixed points. After the rotation, the QFI have to stay at least at the level of the rotation time due to conservation of energy. The effect of noise on the scheme turned out to depend on the direction of the noise. When the noise is along the Y or Z axis, the QFI decreases in time as 1/(N σ y, z ). The scaling law for the time-dependent fidelity to the unperturbed state determines the decay rate of the QFI. On the contrary, the QFI does not decay when the noise is along the X axis, and the level of its stabilization becomes closer to the unperturbed case while increasing the number of particles. We have shown that the parity conservation makes the scheme robust against noise along the X axis. The last finding is quite general since it concerns any noise along the X axis, including a time-dependent γ x .

A similar storage scheme may be applied to another models in vicinity of symmetrically located stable and unstable fixed points. In the one-axis twisting model for example, a single rotation of the NOON state may locate it around stable fixed points leading to a high value of the QFI. A qualitative outcome of the scheme might be the same whenever an additional conservation law is not compatible with the decay mechanism. Nevertheless, generation of desired NOON-like states is strongly hindered by noise. From this point of view the TACT Hamiltonian with the shortened by N timescale is quite promising, although remains an experimental challenge.

After some manipulations, for the first sum we get 1 + e -t 2 σ 2 (N +1) 2 /2 + R 1 (t, σ, N ),

(B5)
where R 1 (t, σ, N ) is an error terms which grows in time and saturates at the level of 0.5, irrespective of N . In fact, one can calculate it exactly in the form of an infinite sum. Similarly, for the second sum we get N n=1 ne -t 2 σ 2 n 2 /2 = 1 -e -t 2 σ 2 (N +1) 2 /2 t 2 σ 2 --1 2 (N + 1)e -t 2 σ 2 (N +1) 2 /2 + R 2 (t, σ, N ), (B6) with the error term R 2 (t, σ, N ) which tends to 0 for growing time. The dominant part of (B1) gives us

F N (t) = 1 λ √ π Erf(λ) - 1 λ 1 -e -λ 2 , ( B7
)
with λ = N σt/ √ 2. Because the error terms depend on time t and function (B7) decays to 0, there exists a time t * after which error terms start to play a significant role. However, they are additionally damped by the size N , so the time t * lengthens once we increase the size N .

FIG. 1 .

 1 FIG. 1. (Color online) Mean-field phase portraits for (a) γx = γy = γz = 0; (b) γx = 0.4, γy = γz = 0; (c) γx = 1, γy = γz = 0; (d) γy = 0.4, γx = γz = 0; (e) γz = 0.4, γx = γy = 0. Red disks visualize the initial spin coherent state.

8 FIG. 2 .

 82 FIG. 2. (Color online) Time evolution of the QFI for γ = 0 (solid red line) and for γ = 0 (gray lines), for N = 50 atoms and M = 2 × 10 3 trajectories. The noise term is nonzero along the Z (a), Y (b), or X (c) axis, with σj = χ (dotted lines), σj = 5χ (dot-dashed lines), σj = 15χ (dashed lines), σj = 50χ (dot-dot-dashed lines). Insets in (c) show the QFI at χt = ln(2πN )/N as a function of σx (left) and γx (right).

FIG. 3 .

 3 FIG. 3. (Color online) (a) Illustration of the storage scheme by classical stable fixed points. The initial state for the evolution is the spin coherent state located around the unstable fixed point. At the optimal moment of time τ = ln(2πN )/2N χ (the first maximum of the QFI) position of the state is shifted by the single rotation e i Ŝxπ/4 to regions around the two stable fixed points. Thereafter, the dynamics is limited to a narrow region of the phase space around the two antipodal stable fixed points. (b) The evolution of the QFI: red dashed line by the TACT Hamiltonian, black solid line by the TACT Hamiltonian with the rotation at t = τ . Both for the noiseless system and N = 50. (c) The QFI with noise and the rotation at τ in time for: σx = χ/10, σy, z = 0 (blue dashed line), σy = χ/10, σx, z = 0 (red solid line), σz = χ/10, σx, y = 0 (green dotted line) and N = 50, M = 2 × 10 3 .

FIG. 4 .

 4 FIG. 4. The Wigner function with noise and the rotation at τ for (a) σx = χ/10, σy, z = 0 and (b) σy = χ/10, σx, z = 0, both at t = 30 ln(2πN )/N χ and N = 50, M = 10 4 . The QFI has the meaning of the speed of change of statistics due to rotation on the Bloch sphere.A small rotation of the state results in a small shifts of the distribution of quasiprobability, here the SU(2) Wigner function[27]. Thus the appearance of small thin structures in the Wigner function leads to a high value of the QFI, as rotating the state perpendicularly to them would strongly affect the state and the quasiprobability distribution. The high narrow fringes, so also the high value of the QFI, remain even in the case with noise along the X axis, but they are mostly smeared out by noise along Y or Z axis.

FIG. 5 .

 5 FIG. 5. (Color online) (a) The half-decay time τ 1/2 of the QFI for γy = 0 and γx, z = 0 (open symbols), and γz = 0 and γx, y = 0 (filled symbols). Dashed line is the best fit τ 1/2 = 0.977(N σy, z ) -1.096 to the numerical data. (b) Averaged values of the QFI after the rotation and noise along X axis. In insets: solid lines are few largest eigenvalues of the system density matrix while dashed line is fidelity function (21) in time. Parameters are the same as in Fig. 3c.

FIG. 6 .

 6 FIG. 6. (Color online) Comparison between the QFI (solid lines) and the Fisher Information F ( Π, n) associated with specific choice of the interferometric direction and measurements of the parity operator (dot-dashed lines). From top to bottom: σx = χ, σx = 10χ and σx = 15χ.

= 0 .

 0 n) one has to find the derivatives dp(µ|θ) dθ θ=0 = -iTr Ŝ n , ρ Pµ . Using the cyclic property of trace we conclude that Tr Ŝ n , ρ P-= 0, which implies dp(µ-|θ) dθ θ=0 As the sum of the probabilities is equal to 1, we have dp(µ+|θ)

ACKNOWLEDGMENTS

We thank O. Hul for a careful reading of the manuscript.

This work was supported by the (Polish) National Science Center Grants DEC-2011/03/D/ST2/01938 and 2014/13/D/ST2/01883.

Appendix A: Fixed points of the phase portrait Below we give explicit analytical expressions for the positions of particular fixed points: (i) |γ x | < 1 and γy = γz = 0 (Fig. 1b). The locations of the four stable center fixed points are (z, ϕ) = (±z x , ϕ x ) and (z, ϕ) = (±z x , π -ϕ x ), where

, while the positions of unstable saddle fixed points do not change as compared to γ = 0 case.

(ii) |γ x | ≥ 1 and γy = γz = 0 (Fig. 1c). The two stable fixed points are located at (z, ϕ) = (0, 0) and (z, ϕ) = (0, -π).

(iii) |γ y | < 1 and γx = γz = 0 (Fig. 1d). The locations of the four stable fixed points are (z y± , ±π/2) and the positions of the two unstable fixed points are (z, ϕ) = (γ y , 0) and (z, ϕ) = (γ y , -π). The locations of unstable fixed points move vertically while changing γy .

(iv) |γ y | ≥ 1 and γx = γz = 0. The two stable fixed points are located at (z, ϕ) = (z y-, ±π/2).

(v) |γ z | < 1 and γx = γy = 0 (Fig. 1e). The locations of the four stable fixed points are (±z z+ , -π/2)

and (±z z-, π/2) while the positions of the two unstable fixed points are (z, ϕ) = (0, arcsinγ z ) and (z, ϕ) = (0, π -arcsinγ z ). The positions of unstable fixed points move horizontally.

(vi) |γ z | ≥ 1 and γx = γz = 0. The two stable fixed points are located at (z, ϕ) = (±z z+ , -sign[γ z ]π/2).

where z y± = (γ y /4) ± γ2 y + 8/4, (A1)

Appendix B: Scaling of the fidelity function with noise along the Z axis

We analyze properties of a general function defined as

where σ is a real constant and p k are positive weights satisfying

There are two discrete sums over the Gaussian functions which can be evaluated using the Euler-Maclaurin integration formula [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. In general one has

where B 2m are the Bernoulli numbers and f (m) (x) denotes the mth derivative of f calculated at x.


Simon Gorecki

Youssef Bouanan

Judicael Ribault

Gregory Zacharewicz

Nicolas Perry
email: nicolas.perry@ensam.eu

INCLUDING CO-SIMULATION IN MODELING AND SIMULATION TOOL FOR SUPPORTING RISK MANAGEMENT IN INDUSTRIAL CONTEXT (a)

Keywords: Model Driven Architecture, Co-Simulation, Functional Mockup-Interface, Risk Management

Due to increasing complexity of engineered system of systems, development of software to design and support them must tend to be more and more concurrent and distributed. To more easily tackle these systems design, global problem is decomposed into several sub-ones where each sub problem is allocated and solved by different contributors. Each participant develops a fragment of the global solution that need after to be integrated with other ones. In this paper we present an extension to the UML/BPMN modeling and simulation tool: Papyrus. This module allows to factor complex tasks during the modelling step and simulation execution process. In detail, we propose to add risk management and other potential interruptions features to BPMN models and Simulation. This is made possible according to Functional Mock-up Interface standard, a cosimulation standard that define how to orchestrate components while simulation execution process.

INTRODUCTION

The Modeling & Simulation (M&S) concept is now a required step in any design of complex systems. It allows to early represent its behavior and interaction. The modeling phase describes a process and allows the development of an executable simulation that virtually designs our subject and anticipates its study. As technologies are growing, systems complexity increases, and makes system more difficult to model and simulate. Along with this growing complexity, risks, hazards, and threat must also be considered during the modeling process. Several research has been done about project risk management [START_REF] References Altuhhova | An Extension of BPMN for Security Risk Management[END_REF][START_REF] Better | Simulation Optimization Applications in Risk Management[END_REF]. In this paper we propose a way to integrate them during the modeling and simulation phases without overloading models. Modeling are one of the primary and most important steps in a project development process. In our case, we need a modeling language able to represent risk management. For that, the Unified Modeling Language (UML) is a general purpose, widely used for describing the different aspects of software and complex systems. It offers several modeling notations to express not only the structure but the behavior of the modeled system. However, this language is also criticized for his low capacity to give precision to the system description. For this reason, the Object Management Group (OMG) has created the concept of UML profile which allows users to create their own UML specification. This UML specification will allow user to define his own language depending on his semantic subject, and provide it to the community: Systems Modeling Language (SysML), Use Case Diagram, etc. This language is a viable solution for representing risks in a system description. Models will be designed and executed with Papyrus tool.

Papyrus is an open source UML/SysML/BPMN modeler of the Eclipse foundation that provides to users and developers a powerful tool for modeling UML models. Another interesting part of Papyrus tool is the execution of UML models due to fUML standards (Semantics of a Foundational Subset for Executable UML Models) allows by the MOKA engine [START_REF] Guermazi | Executable Modeling with FUML and Alf in Papyrus: Tooling and Experiments[END_REF]. Our proposition in this paper consist in defining an extension to MOKA engine able to interrupt the simulation execution and do a request to an external simulation component through Functional Mock-up Interface (FMI). FMI Co-Simulation standard will allow us to relocate and factor complex tasks without overloading models. In the case studied in this paper, an M&S tool is used in a semi-academic, semi-professional context: a French company has launched an innovative project to set up a solar power plant. This project deals with different domains and several simulations tools. All those domains will imply many constraints that must be taken into account during modeling and simulation phases. One of the main aspect we are working on is risk management. Indeed, in the renewable energy domain, many risks and issues must be taken in account by engineer at the modeling phase of the project such as weather issues impact, etc. Our contribution consists in adding a risk management module into an open source modeling and simulation environment: Papyrus.

BACKGROUND

In this section, recent contributions in co-simulation approach are briefly discussed first, efforts in risk management in process modeling are explained in 2.2, Papyrus, a UML modeler is described in 2.3, and finally, Functional Mock-up Interface standard is explained in 2.4 section.

Complex Systems simulation and co-simulation

The emergence of complex engineered systems that integrate both physical, software and network aspects are posing challenges in their design, operation, and maintenance. The current business climate and market pressure are forcing the design of systems to be concurrent, interoperable, distributed and reusable. This is done in order to be divided between different teams and/or external suppliers, each with its own expertise domain and each with its own tools. Here comes the role of Distributed Simulation (DS): one simulation is divided into multiple sub functions (or models) from a large system. Each function is executed on a different computer possibly geographically distributed from others. From a general point of view, this solution divides complex problems into simpler modular sub problems, but also rises interoperability issues. Modeling and Simulation (M&S) of complex systems requires the simultaneous consideration of several points of view. The system behavior has to be considered at different levels and scales. In addition, the study of these systems involves skills from different scientific, business and technical fields. The challenge is then to reconcile these heterogeneous points of view, and to integrate each domain models and tools (or subsystems) within a unified framework, orchestrated by an M&S process. Two of the most popular efforts going in these directions are FMI (Functional Mock-up Interface) and HLA (High Level Architecture). HLA is an IEEE standard (IEEE Computer Society 2010) for distributed computer simulation systems (IEEE Computer Society 2003). In the HLA standard, a distributed simulation is called Federation (see Figure 1). A Federation is composed of several HLA simulation entities, called Federate, which can interact among them by using the Run-Time Infrastructure (RTI). The RTI represents a Federation execution backbone and provides a set of services to manage the communication and data exchange among Federates. FMI (Functional Mock-Up Interface) [START_REF] Blochwitz | Functional Mockup Interface 2.0: The Standard for Tool Independent Exchange of Simulation Models[END_REF]) establishes itself as a standard for model exchange and co-simulation of equational models. The FMI functions are used (called) by a simulation environment to create one or more instances of the FMU (Functional Mockup Unit) and to simulate them, typically together with other models. An FMU may either have its own solvers (FMI for Co-Simulation) or require the simulation environment to perform numerical integration (FMI for Model Exchange). It enforces some generic rules and a software interface to manipulate equational models and their numerical solver using a combination of XML-files and compiled C-code. On that interface, any equational component can be embedded into an FMU (Functional Mock-up Unit) helping to solve the interoperability problem for the co-simulation of equational models. Then, the numerical resolution of a system can be performed by defining a set of communication points between the FMUs according to a trade-off between the accuracy of the simulation results and the performances of the co-simulation process [START_REF] Camus | Hybrid Co-Simulation of FMUs Using DEV&DESS in MECSYCO[END_REF]. The FMI standard defines two interfaces: FMI for Model Exchange and FMI for Co-Simulation [START_REF] Blochwitz | Functional Mockup Interface 2.0: The Standard for Tool Independent Exchange of Simulation Models[END_REF]. The FMU CS contains its own solver that will be built when generating the tool. The advantage of this model is to combine two or more simulation tools in a co-simulation environment. The exchange of data between the subsystems is limited to "Communication Points". Between two Communication Points, the subsystems are solved independently from each other by their individual solver. Within a master-slave view, slaves simulate subproblems while the master is responsible for the coordination of the overall simulation and data transfer. Several tools are compatible with the FMI interface at Export/Import for both components FMU, Model-Exchange (ME) and co-simulation (CS). Example: JModelica, Dymola, LMS AMESim, EnergyPlus, CATIA, NI LabVIEW, Ptolemy II, etc.

Risk management in process modeling

The concept of risk is highly polysemous and supports a large number of definitions. In the context of risk management, we can introduce a number of concepts revolving around risk and conditioned by the environment and the components of the project. The project risk is related to the occurrence of events, from internal or external origin, which may affect the achievement of the initial target. The risk qualify the effect of these events on the achievement of project's objectives. The anticipation of these events via the factor's identification, internal or external, which are the cause, the evaluation of their impact on the project progress and the proposal of appropriate treatment actions are the purpose of risk management. In literature, we observed most recurrent major steps in risk management that we can cite: identification, analysis, evaluation and treatment of risks. These keywords can be used in tools such as brainstorming guidelines in order to anticipate, minimize risks in a project. However, these methods are not much structured, mostly handle qualitative information and are frequently limited to user experience and point of view. In the following, we cite some relevant methods and tools that were proposed in literature to manage risks.

 Information

Systems Security Risk Management (ISSRM) methods and standards (according to [START_REF] Dubois | A Systematic Approach to Define the Domain of Information System Security Risk Management[END_REF])) mainly consist of process guidelines that help identify vulnerable assets, determine security objectives, and assess risks as well as define and implement security requirements to treat the risks. By using these methods one reduces the losses that might result from security problems. However, these methods generally offer quite poor modelling support. Instead, they usually resort to informal documentation in natural language and ad hoc diagrams. According to Figure 1, this analysis method allow user to identify and classify risks into three categories: -asset-related concepts are used to identify in a systems, or in a company, skills of the organization, and security risks that must be avoid. -risk-related concepts identify risk events and threat related to assets defined previously. -Risk-treatment concepts are defined depending on risks possibilities. They will represent the decision of how to treat the identified risks. A treatment satisfies a security need, expressed in generic and functional terms, and can lead to security requirements.

 Business Process Modeling Notation Extension

for Risk Handling was proposed in [START_REF] Marcinkowski | A Business Process Modeling Notation Extension for Risk Handling[END_REF]. The paper identify three different risk types: -Businessdriven risks, strategic in nature aimed at protecting the business and keeping it accessible whenever and whoever in support of continuous business operations. -Data-driven risks, dealing with the availability of data and information in all of its different forms as used by the organization. -Event-driven risks, focusing on actual events that create risk to business continuity and viability. It propose to extend BPMN standard in order to modeling the several risks and handle it in three different ways: reduce it, retain it, avoid it, or transfer it.

Papyrus

Papyrus is the UML/SysML modeler of the Eclipse foundation. It provide tools for executing and debugging UML models. The execution part is handled by MOKA, a fUML execution engine. Papyrus is based on Eclipse and is open source. In accordance with its primary goal to implement the complete standard specification of UML2, Papyrus provides an extensive support for UML profiles. It includes hence all the facilities for defining and applying UML profiles in a very rich and efficient manner. But, it also provides powerful tool customization capabilities similar to DSML-like metatools. This way, Papyrus is a tool enabling to gather the advantages of using a general purpose language such as UML2, but also those of DSML-based approaches.

FMI

The Functional Mock-up Interface (FMI) for Co-Simulation interface is designed both for the coupling of simulation tools (simulator coupling, tool coupling), and coupling with subsystem models, which have been exported by their simulators together with its solvers as runnable code. It is an interface standard for the solution of time dependent coupled systems consisting of subsystems that are continuous in time or time-discrete [START_REF] Bastian | Master for Co-Simulation Using FMI[END_REF][START_REF] Blochwitz | Functional Mock-up Interface for Model Exchange and Co-Simulation[END_REF][START_REF] Sievert | Modelica Models in a Distributed Environment Using FMI and HLA[END_REF]. It provides interfaces between master and slaves and addresses both data exchange and algorithmic issues. FMI for Co-Simulation consists of two parts (Figure 2):

 Co-Simulation Interface: a set of C functions for controlling the slaves and for data exchange of input and output values as well as status information.

 Co-Simulation Description Schema: defines the structure and content of an XML file. This slave specific XML file contains "static" information about the model (input and output variables, parameters …) and the solver/simulator (capabilities …). The capability flags in the XML file characterize the ability of the slave to support advanced master algorithms which use variable communication step sizes, higher order signal extrapolation etc.

Integrating risk management in modeling and simulation process

In order to integrate concepts described in section 2.2 one of the first solution was to implement risks treatment directly on models. However, risk management, particularly in renewable energy production domain imply different and complex methodologies and make the modeling too heavy and not systematic for the user. Modeling and simulation a global solar energy production site with this method of risk management is too complex for user. According to our subject, we must manage risks depending on several levels:  A risk can block a task, or it can slow it down (depending on a degree) as described in Figure 4.

Figure 4. Risk impact on a task  A risk can influence another one (increase or decrease its impact), it can form a causality chain such as described in Figure 5. Those risk management mechanisms appear too complex and heavy for being implemented directly on process model, also the rules can be dependent on the domain studied. It forms the reason why we propose to outsource them from papyrus model designer.

Proposition

We propose to define an extension to a modeling and simulation tool to include risk treatment separated from the main model. In our case, Papyrus, based on Eclipse, allow user to model UML diagram and, with the mechanism of UML-profiles, enable UML-based diagram modeling like SysML or BPMN. One of the advantages of using Papyrus as a modeling tool is the included MOKA execution engine which make possible to simulate fUML models.

Our goal is to implement an ad-hoc system (Figure 7) to Papyrus in order to manage constraints, risks, and hazards during the simulation execution. This module is able to generate issues in the main system described by the process model. All the potentials issues are generated according to equations defined outside of the process model. This risk management extension is disconnected from the initial model in order to keep it clean (see section 4.1). The global simulation execution will react depending on hazards and constraints generated and referenced by the risk management module (as described in Figure 7). The objective of this proposition is to relocate and factorize complex tasks (risk management) during the simulation execution process. This may involve to use a co-simulation standard for insure communication between different external modules. Another way to communicate with the outside is to request web services API. According to this proposition, every type of complex systems can be relocated outside of the model. FMI standard, propose two mechanism: the first is "simulation exchange", the one that we are using in this paper. But the standard propose also "model exchange" mode which could be also implemented in this context. Using FMI as model exchange would means to modeling process into an external FMU and import it in the simulation. However, the MOKA engine could not execute fUML models described in external models.

BUSINESS CONTEXT

In our context, a company designing solar power plants has special needs. This project consists in installing solar panels fields to provide electricity and heat in a large area which is not powered so far. However, the transport of solar panels fields is extremely expensive. To reduce this cost/blow, they are designing a mobile factory which will be able to manufacture and assemble the solar panels on site. Rather than transporting finished products, only the mobile plant and raw materials would be carried out. The main challenges of this project are : the factory miniaturization to fit in the least transport containers (around 20) [START_REF] Benama | Supporting Make or Buy Decision for Reconfigurable Manufacturing System, in Multi-Site Context[END_REF], risks management [START_REF] Rodney | Integrating Risks in Project Management[END_REF] caused by low knowledge, and designing resistant structure depending on the environment of the power plant [START_REF] Piegay | Multi-Objective Optimization and Decision Aid for Spread Footing Design in Uncertain Environment[END_REF], (El Amine 2016).

To guide the project and model the conception process, the company use Papyrus tool. The objective is to create several models for modeling and simulation of all workflow from the deployment of the mobile factory, to the management in real time of the production site. With the implementation of risk management module, the company aim at running simulations for every case possible in order to being able to anticipate and avoid problems.

REALISATION

According to [START_REF] Guermazi | Executable Modeling with FUML and Alf in Papyrus: Tooling and Experiments[END_REF], when the MOKA engine execute a fUML model, a mechanism named visitors is executed at each steps of the simulation process. This visitor can be surcharged and can execute java code at each steps of the simulation process. This is our entry point to the MOKA engine. During the execution process, at each task parsed by the engine, a visitor is called, we can test the name of the task and have access to his parameters (that is the "Reference" link between "Process Model" and "Risk & Hazard management" in Figure 7). During a global simulation process, the risk management module must know the name of each task of the model for allowing interactions with it.

The entry point of the MOKA engine is a java function.

From it, we can integrate a co-simulation environment with Java-FMI. The library allow us to load an FMU file in order to interact with another simulation. In our context, we build a simple FMU which can generate errors in the process model according to a normal law (see Figure 8). In the same time, we are requesting weather information to a web service (OpenWeatherMap). Location is given in input of the API and it return to the visitor weather information. In our context, for a solar system, electricity production yield will be impacted depending on clouds and sunshine. We also add the implementation of the two advices (advice for FailureResource and advice for FailureAction). The first advice is associated to tasks which use Resources with FailableResource stereotype. It will ask to the failure manager if a task can start or finish and it =will also execute additional start/finish actions for a given task. Concerning the second advice, it is associated to Tasks with the FailableAction stereotype. It verify also on the failure manager if a task can start or finish and is executes additional start/finish actions for a given task. The part that allows to import an FMU and to execute it is implemented within these advices.

CONCLUSION

In modeling and simulation domain, many efforts are done to increase reusability and creating bridges between technologies. In this paper we presented a contribution to extend Papyrus execution engine for FMI co-simulation in the context of risk management. The paper demonstrate also the link between BPMN models to risk description. Then these risks characteristics are used at the simulation step. However, many aspects must be improved for our contest. It is necessary to declare risks and hazards equations linked to our model. Efforts must be done to increase usability of this extension for enable model exchange aspects of FMI.

Figure 1 :

 1 Figure 1: ISSRM Risk Data Model

Figure

  Figure 2: FMU / FMI Concepts 3. CONTRIBUTION

Figure 3 :

 3 Figure 3: Risks Modeling in BPMN We can see on the top of the above Figure 3 a BPMN model of a simple database login action between a user and a system. The BPMN diagram at the bottom of Figure 3 represent the same upper model, with risk treatment taken in account: model can be very complex.Modeling and simulation a global solar energy production site with this method of risk management is too complex for user. According to our subject, we must manage risks depending on several levels:

Figure 5 .

 5 Figure 5. Risk causality chain  Several risks can be composed to generate or implicate a new hazard as described in Figure 6.

Figure

  Figure 6: Composed Risk

Figure 7

 7 Figure7represents the proposed architecture, we can observe the process model editor and the MOKA execution engine combined into Papyrus, and a connection to an external risk management module. This new module is connected to both Process Modeler and Moka Execution module to interact with them. The objective of this proposition is to relocate and factorize complex tasks (risk management) during the simulation execution process. This may involve to use a co-simulation standard for insure communication between different external modules. Another way to communicate with the outside is to request web services API. According to this proposition, every type of complex systems can be relocated outside of the model.

Figure 7 :

 7 Figure 7: Defining and Using Risks in M&S

Figure 8 :

 8 Figure 8: Requesting FMU and Web Service from Papyrus VisitorUsing MOKA engine to execute the models allows us the implementation of preconditions for tasks and resources.

Figure 9 :

 9 Figure 9: UML Profile describing the new concepts  FailableAction stereotype: Applied on Tasks (UML actions) with two attributes: (i) Mean Time Between Failure defines the average time (in hours) during an action can be executed before a failure occurs over a specified time period. (ii) Mean Time To Repair represents the average time required to repair a failed action (during this time the task can't be executed). (iii) Notification of Default suspends the execution of the task depending on external variables and indicators.  FailableResource stereotype: Applied on Resources (UML Class) with three attributes: (i) Mean Time Between Failure defines the average time (in hours) during the resource can be used before a failure occurs over a specified time period. (ii) Mean Time To Repair represents the average time required to repair a failed resource (during this time the resource can't be allowed or used by tasks independently from Availability constraints). (iii) Notification of Default suspends the availability of the resource depending on external variables and indicators.The second step consists in generating the source code of the profile using the Eclipse Modeling Framework (EMF). The code generator for EMF models can be adjusted and in its default setting. It provides change notification functionality to the model in case of model

Figure 10 :

 10 Figure 10: Code generated using EMF Then, we customize the execution engine by adding the implementation of the new execution engine which associate new Advices to the execution visitors, interacting with a FailureManager. The FailureManager centralized class counting the number of FailableElement executions and activating failures.We also add the implementation of the two advices (advice for FailureResource and advice for FailureAction). The first advice is associated to tasks which use Resources with FailableResource stereotype. It will ask to the failure manager if a task can start or finish and it =will also execute additional start/finish actions for a given task. Concerning the second advice, it is associated to Tasks with the FailableAction stereotype. It verify also on the failure manager if a task can start or finish and is executes additional start/finish actions for a given task. The part that allows to import an FMU and to execute it is implemented within these advices.