Chuan Xu

Jianqiang Cheng

Abdel Lisser

A Sampling Method to Chance-constrained Semidefinite Optimization

Keywords: stochastic programming, chance-constrained programming, sample approximation, semidefinite program

Semidefinite programming has been widely studied for the last two decades. Semidefinite programs are linear programs with semidefinite constraint generally studied with deterministic data. In this paper, we deal with a stochastic semidefinte programs with chance constraints, which is a generalization of chance-constrained linear programs. Based on existing theoretical results, we develop a new sampling method to solve these chance constraints semidefinite problems. Numerical experiments are conducted to compare our results with the state-of-the-art and to show the strength of the sampling method.

INTRODUCTION

It is well known that optimization models are used for decision making. In the traditional models, all the parameters are assumed to be known, which conflicts with many real world problems. For instance, in portfolio problems, the return of assets are uncertain. Further, real world problems almost invariably include some unknown parameters. Therefore, the deterministic optimization models are inadequate and a new optimization model is needed to tackle the uncertainty. In this case, stochastic programming is proposed to handle the uncertainty.

As a branch of stochastic programming, chanceconstrained problem (CCP) which is called probabilistic problem as well, was first proposed in [START_REF] Charnes | Cost horizons and certainty equivalents: an approach to stochastic programming of heating oil[END_REF] to deal with an industrial problem. The authors considered a special case of CCP where the probabilistic constraints are imposed individually on each constraint. Latter, [START_REF] Prékopa | On probabilistic constrained programming[END_REF] generalized the model of CCP with joint probabilistic constraints and dependant random variables. See [START_REF] Dentcheva | Concavity and efficient points of discrete distributions in probabilistic programming[END_REF][START_REF] Prékopa | Probabilistic programming[END_REF][START_REF] Henrion | Convexity of chance constraints with independent random variables[END_REF] for a background of CCP and some convexity theorems. In order to circumvent CCP, we usually consider tractable approximation.

For instance, convex approximation (Nemirovski and Shapiro, 2006a;[START_REF] Nemirovski | On safe tractable approximations of chance constraints[END_REF] is a way which analytically generates deterministic convex problems which can be solved efficiently. However, it requires the known structure of the distribution and structural assumptions on the constraints. Another way is simulation-based approach based on Monte-Carlo sampling, for example the well-known scenario approach [START_REF] Calafiore | Uncertain convex programs: randomized solutions and confidence levels[END_REF][START_REF] Calafiore | The scenario approach to robust control design[END_REF][START_REF] Nemirovski | Scenario approximations of chance constraints[END_REF]. As the sampling number N is large enough, we can ensure the feasibility of the solution. In [START_REF] Campi | A sampling-anddiscarding approach to chance-constrained optimization: feasibility and optimality[END_REF], the authors developed a sampling-anddiscarding approach which removes some sampling constraints in the model. They gave theoretical proofs where discarding suitable number of constraints in the sampling model, the result remains feasible and intact. A greedy algorithm to select the constraints to be removed was mentioned and some numerical results are shown in [START_REF] Pagnoncelli | Risk-return trade-off with the scenario approach in practice: a case study in portfolio selection[END_REF]. Recent work of [START_REF] Garatti | Modulating robustness in control design: Principles and algorithms[END_REF] presented a precise procedure of this algorithm on control design.

The probabilistic problem that we work on is the minimum-volume invariant ellipsoid problem in control theory which can be formulated as semidefinite program with chance constrains (CCSDP). In [START_REF] Cheung | Linear matrix inequalities with stochastically dependent perturbations and applications to chance-constrained semidefinite optimization[END_REF], authors proposed a convex safe tractable approximation to solve this problem. In our work, we develop a simulation-based method base on [START_REF] Campi | A sampling-anddiscarding approach to chance-constrained optimization: feasibility and optimality[END_REF]. For the related work to CCSDP, we refer the reader to [START_REF] Yao | LQ Control without Riccati Equations: Stochastic Systems[END_REF][START_REF] Ariyawansa | Chance-constrained semidefinite programming[END_REF][START_REF] Zhu | Semidefinite programming under uncertainty[END_REF].

The paper is organised as follows. In section 2, we present mathematical formulation of the chance constrained semidefinite problem. In section 3, we present simulation-based methods applied on semidefinite program with chance constraints and introduce our method of sampling. In section 4, we show numerical results on the problem in control the-ory. Finally, a conclusion is given in section 5.

CHANCE CONSTRAINED SEMIDEFINITE PROGRAM

Conic optimization problems with chance constraints can be generalized as

(CCP) min{ f (x) : Pr{F(x, ξ) ∈ K} ≥ 1 -ε, x ∈ X} where x ⊆ R n is a vector of decision variables, X is a deterministic feasible region, ξ is a random vec- tor supported by a distribution Ξ ⊆ R d , K ⊂ R l is a closed convex cone, F : R n , R d → R l is
F(x, ξ) = A 0 (x) + m ∑ i=1 ξ i A i (x) + ∑ 1≤ j≤k≤m ξ j ξ k B jk (x)
where A i , B jk are symmetric matrix. Therefore, the chance constrained semidefinite program can be presented as:

(CCSDP) min{ f (x) x∈X : Pr{F(x, ξ) 0} ≥ 1 -ε} 3 SIMULATION-BASED APPROXIMATION 3.1 Scenario Approach
The simplest method of simulation-based approximation is scenario approach. The approximation of CCSDP is:

(CCP -SA) min{ f (x) x∈X : F(x, ξ i ) 0, ∀i = 1, ..., N}
where N is the number of sampling, ξ i is a random sample. CCP -SA yields a feasible solution to CCSDP with probability of at least 1β for

N ≥ 2 ε log( 1
β ) + 2n + 2n ε log( 2 ε ).
( [START_REF] Calafiore | The scenario approach to robust control design[END_REF]

Big-M Semidefinite Sampling Approach

In [START_REF] Luedtke | A sample approximation approach for optimization with probabilistic constraints[END_REF], the authors proposed a simulation-based method which adds a sample average constraint involving expectations of indicator functions. They showed that their simulation-based approximation method yields a feasible solution to the chance constrained problem with high confidence.

If we choose "big-M" function with integer variables to be the indicator function, we have the following tractable approximation of CCSDP:

(CCP -BM) min f (x) s.t F(x, ξ i ) + y i MI 0, ∀i ∈ 1, ..., N N ∑ i=1 y i ≤ ε × N x ∈ X, y ∈ {0, 1} N
where I is an identity matrix, M is a large constant. We see that if y i = 1, the constraint i is satisfied for any candidate solution x including those x ∈ {x|F(x, ξ i ) 0, x ∈ X} discarded by scenario approach (CCP-SA). This "big-M" method is less conservative than CCP -SA, but it introduces the binary variables which increases the computation effort. The advantage of this method is that it gives a less conservative solution.

Combination of Big-M and Constraints Discarding

In order to have a less conservative solution than the scenario approach and reduce the computation effort, our sampling method starts by solving a relaxed CCP -BM model. As we suppose that the relaxed values of y could help select the constraints to be removed in sampling-and-discarding approach proposed by [START_REF] Campi | A sampling-anddiscarding approach to chance-constrained optimization: feasibility and optimality[END_REF].

In our method, we suppose that the relaxed value of y i ∈ [0, 1] obtained by the relaxed CCP -BM indicates the probability of discarding the constraint i. Therefore, we develop a new sampling method which combines the "big-M" approximation and samplingand-discarding method. The main procedure is that we solve the relaxed CCP -BM at first and then according to the sorted value of y i , remove the corresponding constraints in CCP -SA and solve the new reduced problem.

We apply our method to a minimum-volume invariant ellipsoid problem in control theory [START_REF] Cheung | Linear matrix inequalities with stochastically dependent perturbations and applications to chance-constrained semidefinite optimization[END_REF] and compare the performance with scenario approach, sampling-and-discarding approach with greedy procedure [START_REF] Pagnoncelli | Risk-return trade-off with the scenario approach in practice: a case study in portfolio selection[END_REF].

Control System Problem

First of all, we state out the problem and its mathematical model. Supposed that we have the following discrete-time controlled dynamical system:

x(t + 1) = Ax(t) + bu(t) t = 0, 1, ... x(0) = x
where A ∈ R n×n and b ∈ R n are system specifications, t is the index of discrete time, x is the initial state, and u(t) is the control at time t. In order to keep the system stable for any A, b and possible u(t), the safe region for x could be an invariant ellipsoid. An ellipsoid is expressed by:

E(Z) = {x ∈ R n : x T Zx ≤ 1}
where Z is a symmetric positive definite matrix.

An invariant ellipsoid means that if x ∈ E(Z), then A(x) + b ∈ E(Z). [START_REF] Nemirovski | Lectures on modern convex optimization[END_REF]) has shown that the ellipsoid E(Z) is invariant if and only if there exists a λ ≥ 0 such that 1

-b T Zb -λ -b T ZA -A T Zb λZ -A T ZA 0, ||A|| < 1. For
this problem, we prefer to have a smaller safe region for x to ensure the stability. Thus, this control problem is equivalent to finding the minimum-volume of an invariant ellipsoid which could be formulated as a bilinear semidefinite programming problem. If we considered the chance constrained case, the model should be : {CCMV IE(λ), λ ∈ D} [START_REF] Cheung | Linear matrix inequalities with stochastically dependent perturbations and applications to chance-constrained semidefinite optimization[END_REF].

CCMV IE(λ) : max w s.t w ≤ (detZ) 1/n ) Pr 1 -b T Zb -λ -b T ZA -A T Zb λZ -A T ZA 0 ≥ 1 -ε Z 0
where D = {0.00, 0.01, ..., 0.99, 1.00} is a finite set. We assume that the system could be disturbed by some random noise. Like the design of numerical experiment in [START_REF] Cheung | Linear matrix inequalities with stochastically dependent perturbations and applications to chance-constrained semidefinite optimization[END_REF], b is corrupted and b i = bi + ρξ i , ∀i = 1, ..., N where b ∈ R N is the nominal value, ρ ≥ 0 is a fixed parameter to control the level of perturbation, ξ i is a standard Gaussian random variable of sample i.

Sampling Procedure

Scenario Approach

We generate N random samples and solve the following model {CCSC(λ), λ ∈ D} :

CCSC(λ) : max w s.t w ≤ (detZ) 1/n ) 1 -b T i Zb i -λ -b T i ZA -A T Zb i λZ -A T ZA 0, ∀i = {1, ..., N} Z 0

Greedy Procedure for Sampling-and-Discarding Method

For each λ ∈ D, we apply a greedy and randomized constraint removal procedure [START_REF] Pagnoncelli | Risk-return trade-off with the scenario approach in practice: a case study in portfolio selection[END_REF] to the sample counterpart (SP) of CCMV IE(λ) [START_REF] Campi | A sampling-anddiscarding approach to chance-constrained optimization: feasibility and optimality[END_REF].

CCSP(λ) : max w s.t w ≤ (detZ) 1/n ) 1 -b T i Zb i -λ -b T i ZA -A T Zb i λZ -A T ZA 0, ∀i = {1, ..., N} -A Z 0
where A is the set of the indexes of the k removed constraints.

The greedy removal procedure iteratively removes k constraints. At each iteration i, we solve a CCSP(λ) with A i-1 to determine the set of n i active constraints.

Then we randomly choose one of these active constraints such as constraint c to have A i = A i-1 ∪ {c} for following iteration i + 1.

Big-M Procedure for Sampling and Discarding Method

Our sampling method contains two parts. First, we solve a relaxed "big-M" model CCRBM(λ) and obtain the solution of the relaxed binary variable y:

CCRBM(λ) max w s.t w ≤ (detZ) 1/n ) 1 -b T i Zb i -λ -b T i ZA -A T Zb i λZ -A T ZA + y i MI 0, ∀i = 1, ..., N N ∑ i=1 y i ≤ ε × N Z 0 0 ≤ y i ≤ 1, ∀i = 1, ..., N
We sort the elements of y in descending order and take the first k indexes into set A = {i 1 , ..., i k }.

Then, we solve CCSP(λ).

Precise procedure:

1. For each λ ∈ D:

(a) Solve CCRBM(λ) and obtain the relaxation values of y, (b) Determine the set A of removed constraints according to y, (c) Solve CCSP(λ), and let v(λ) be the objective value and Z(λ) be the corresponding solution.

2. Return Z(λ * ) as the optimal solution, where λ * = argmax λ∈D v(λ).

Design of The Experiments

Data

We use the same instances as [START_REF] Cheung | Linear matrix inequalities with stochastically dependent perturbations and applications to chance-constrained semidefinite optimization[END_REF].

We have two group of data.

Data 1 : A = -0.8147 -0.4163 0.8167 -0.1853 , b = 1 0.7071 , ε = 0.05, ρ = 0.01, β = 0.05 Data 2 : A =      0 2 0 0 0 0 0 0.0028 0.0142 0 0 0 0 1 0 0 0 -0.0825 -0.4126 0 1 0 0 0 0      , b =      0 0.0076 0 -0.1676 0      , ε = 0.03, ρ = 0.001, β = 0.05
where β is a confidence parameter which is needed to decide the sample size N and number of removal constraints k.

Selecting the Sample Size and the Number of Constraints to be Removed

For data 1, we consider four sample sizes N ranging from 400 up to 1000. The number of constraints to be removed is calculated as following :

k = εN -d + 1 -2εIn (εN) d-1 β ,
where d is the dimension of variable Z. It has been proven in [START_REF] Campi | A sampling-anddiscarding approach to chance-constrained optimization: feasibility and optimality[END_REF] that with this number of k, the solution obtained by CCSP(λ) (with optimal removal) is feasible to CCMV IE(λ) with high probability 1β.

As choosing the optimal set of constraints to be removed is an NP-hard problem, the solution that we obtain with our procedure can not ensure conservativeness. Therefore, we vary the ratio of k/N from 0.03 to 0.05 to study the influence of k on the result.

For data 2, we consider three sample sizes N ranging from 1000 to 1400 with k calculated as in (4.3.2). In addition, we set the ratio of k/N to be 0.02 and 0.03 for each sample size respectively.

Numerical Results

All experiments are run under MATLAB R2012b on a Windows 7 operating system with i7 CPU 2GHz and 4GB of RAM. The computations are performed using CVX 2.1 with semidefinite program solver SeDuMi.

Tables [1] and[2] provide the computational results of Data 1 and Data 2 respectively. N presents the sampling number. k is the number of removal constraints and k/N is the corresponding ratio. We use the average linear size measure, which is defined as ALS(E(Z)) = (Vol n (E(Z)) 1/n ) [START_REF] Cheung | Linear matrix inequalities with stochastically dependent perturbations and applications to chance-constrained semidefinite optimization[END_REF], to evaluate the volume of ellipsoid. The smaller the volume of ellipsoid is, the smaller the average linear size of ellipsoid is. The columns SC, Greedy, BMSP give the average linear size of ellipsoid obtained by scenario approach (4.2.1), greedy approach (4.2.2) and our method (4.2.3) respectively. 1 -Vio shows the satisfaction rate of each solution estimated under 100000 simulated random samples. Gap presents the gap between the solution of the current method and the solution of the scenario approach. Table [3] shows the CPU time expressed in seconds. The columns SC, Greedy, BMSP show the average CPU time of all tests in Table [1] and[2] when applying scenario approach, greedy approach and our We observe that the real violation is significantly below 5% and 3% respectively in Tables [1] and[2]. It is easy to see that as k increases, we obtain a better solution both with greedy method and with our BMSP method; and the violation of the solution is larger. The reason is that as the more constraints we remove, the larger feasible set of CCSP(λ) we obtain, which involves more violated elements of CCSC(λ).

In Table [1], for each sampling number N, BMSP obtains better solution than Greedy with smaller final value (average linear size of ellipsoid) and a larger violation which is below 5%. For greedy method, the gap is between 0.5‰-10.4‰, compared with scenario approach. While for our method, the gap is between 2.7‰-11.8‰. Figure [1] gives a precise look on the final value obtained by Greedy and BMSP for different values of k for 400 samples. In Figure [2], we compare the violation of Greedy and BMSP. We observe that the increasing rate of violation is nearly the same. Figure [3] shows the local view of ellipsoid for Data 1 obtained by scenario approach, greedy approach and our method with N = 400 and k = 20. We can see that the ellipsoid obtained by our method has the smallest volume.

In Table [2], we obtain a Gap more obvious than the previous one on Data 1. For the case where k is chosen by (4.3.2), our method obtains a gap better than Greedy method with 0.2% to 0.6% improvement. While for other choices of k, their gap are very close to each other.

The advantage of our method compared with Greedy procedure is on the computing time. In the Greedy procedure, we need to solve k + 1 times the semidefinite program CCSP(λ) in order to decide removal constrains, while in our method, we only need to solve 2 semidefinite programs. Therefore, we observe from Table [3] that BMSP consumes much less CPU time than Greedy and almost twice CPU time than scenario approach. But as a counterpart of the CPU time, we obtain better solution than scenario approach.

CONCLUSION

In this paper, we introduce a new simulation-based method to solve stochastic chance constrained program. This method is a combination of Big-M relaxation and a sampling-and-discarding method. We apply this method to semidefinite programming problem in control theory. The numerical results show that our method provides better solutions within a reasonable CPU time.

Figure 1 :

 1 Figure 1: Comparison of average linear size for Data 1

Figure 3 :

 3 Figure 3: Local view of chance-constrained invariant ellipsoid of Data 1 with N = 400, k = 20.


Sahbi El Hog

H T Diep

Henryk Puszkarski

Theory of magnons in spin systems with Dzyaloshinskii-Moriya interaction

Keywords: numbers: 75.25.-j, 75.30.Ds, 75.70.-i :

We study in this paper magnetic properties of a system of quantum Heisenberg spins interacting with each other via a ferromagnetic exchange interaction J and an in-plane Dzyaloshinskii-Moriya interaction D. The non-collinear ground state due to the competition between J and D is determined. We employ a self-consistent Green'function theory to calculate the spin-wave spectrum and the layer magnetizations at finite T in two and three dimensions as well as in a thin film with surface effects. Analytical details and the validity of the method are shown and discussed.

I. INTRODUCTION

The Dzyaloshinskii-Moriya (DM) interaction was proposed to explain the weak ferromagnetism which was observed in antiferromagnetic Mn compounds. The phenomenological Landau-Ginzburg model introduced by I. Dzyaloshinskii [1] was microscopically derived by T. Moriya [2]. The interaction between two spins S i and S j is written as

D i,j • S i ∧ S j (1)
where D i,j is a vector which results from the displacement of non magnetic ions located between S i and S j , for example in Mn-O-Mn bonds. The direction of D i,j depends on the symmetry of the displacement [2]. The DM interaction is antisymmetric with respect to the inversion symmetry.

There has been a large number of investigations on the effect of the DM interaction in various materials, both experimentally and theoretically for weak ferromagnetism in perovskite compounds (see references cited in Refs. 3 and 4, for example). However, the interest in the DM interaction goes beyond the weak ferromagnetism: for example, it has been recently shown in various works that the DM interaction is at the origin of topological skyrmions [5][6][7][8][9][10][11][12][13][14][15][16] and new kinds of magnetic domain walls [START_REF] Heide | Dzyaloshinskii-Moriya interaction accounting for the orientation of magnetic domains in ultrathin films: Fe/W(110)[END_REF][START_REF] Rohart | Skyrmion confinement in ultrathin film nanostructures in the presence of Dzyaloshinskii-Moriya interaction[END_REF]. The increasing interest in skyrmions results from the fact that skyrmions may play an important role in the electronic transport which is at the heart of technological application devices [19].

In this paper, we are interested in the spin-wave (SW) properties of a system of spins interacting with each other via a DM interaction in addition to the symmet-ric isotropic Heisenberg exchange interaction. The competition between these interactions gives rise to a noncollinear spin configuration in the ground state (GS). Unlike helimagnets where the helical GS spin configuration results from the competition between the symmetric nearest-neighbor (NN) and next-nearest neighbor (NNN) interactions [20,21], the DM interaction, as said above, is antisymmetric. This gives rise to a non trivial SW behavior as will be seen below. Note that there has been a number of works dealing with the SW properties in DM systems [START_REF] Puszkarski | Effect of Dzialoshinsky-Moriya Interactions on Propagation of Spin Waves in Ferromagnets: Dynamical Canting[END_REF][START_REF] Kh | Asymmetric Spin-Wave Dispersion on Fe(110): Direct Evidence of the Dzyaloshinskii-Moriya Interaction[END_REF][START_REF] Wang | Magnon-Driven Domain-Wall Motion with the Dzyaloshinskii-Moriya In-teraction[END_REF][START_REF] Stashkevich | Experimental study of spin-wave dispersion in Py/Pt film structures in the presence of an interface Dzyaloshinskii-Moriya interaction[END_REF][START_REF] Moon | Spin-wave propagation in the presence of interfacial Dzyaloshinskii-Moriya interaction[END_REF].

This paper is organized as follows. Section II is devoted to the description of the model and the determination of the GS. Section III shows the formulation of our selfconsistent Green's function (GF) method. Section IV shows results on the SW spectrum and the magnetization in two dimensions (2D) and three dimensions (3D). The case of thin films with free surfaces is shown in section V where layer magnetizations at finite temperature (T ) and the thickness effect are presented. Concluding remarks are given in section VI.

II. MODEL AND GROUND STATE

We consider a thin film of simple cubic (SC) lattice of N layers stacked in the y direction perpendicular to the film surface. For the reason which is shown below, we choose the film surface as a xz plane. The Hamiltonian is given by

H = H e + H DM (2)
H e = - i,j J i,j S i • S j (3)
H DM = i,j D i,j • S i ∧ S j (4)
where J i,j and D i,j are the exchange and DM interactions, respectively, between two Heisenberg spins S i and S j of magnitude S = 1/2 occupying the lattice sites i and j.

For simplicity, let us consider the case where the inplane and inter-plane exchange interactions between NN are both ferromagnetic and denoted by J and J ⊥ , respectively. The DM interaction is supposed to be between NN in the plane with a constant D. Due to the competition between the exchange J term which favors the collinear configuration, and the DM term which favors the perpendicular one, we expect that the spin S i makes an angle θ i,j with its neighbor S j . Therefore, the quantization axis of S i is not the same as that of S j . Let us call ζi the quantization axis of S i and ξi its perpendicular axis in the xz plane. The third axis ηi , perpendicular to the film surface, is chosen in such a way to make ( ξi , ηi , ζi ) an orthogonal direct frame. Writing S i and S j in their respective local coordinates, one has

S i = S x i ξi + S y i ηi + S z i ζi (5)
S j = S x j ξj + S y j ηj + S z j ζj (6)

We choose the vector D i,j perpendicular to the xz plane, namely

D i,j = De i,j ηi (7)
where e i,j =+1 (-1) if j > i (j < i) for NN on the ξi or ζi axis. Note that e j,i = -e i,j .

To determine the GS, the easiest way is to use the steepest descent method: we calculate the local field acting on each spin from its neighbors and we align the spin in its local-field direction to minimize its energy. Repeating this for all spins and iterating many times until the convergence is reached with a desired precision (usually at the 6-th digit, namely at ≃ 10 -6 per cents), we obtain the lowest energy state of the system (see Ref. [START_REF] Ngo | Effects of frustrated surface in Heisenberg thin films[END_REF]. Note that we have used several thousands of different initial conditions to check the convergence to a single GS for each set of parameters. Choosing D i,j lying perpendicular to the spin plane (i. e. xz plane) as indicated in Eq. ( 7), we determine the GS as a function of D. An example is shown in Fig. 1 for θ = π/6 (D = -0.577) with J = J ⊥ = 1. We see that each spin has the same angle with its four NN in the plane (angle between NN in adjacent planes is zero). Let us show the relation between θ and J : the energy of the spin S i is written as where θ = |θ i,j | and care has been taken on the signs of sin θ i,j and e i,j when counting NN, namely two opposite NN have opposite signs. The minimization of E i yields

E i = -4J S 2 cos θ -2J ⊥ S 2 + 4DS 2 sin θ (8)
dE i dθ = 0 ⇒ - D J = tan θ ⇒ θ = arctan(- D J ) (9)
The value of θ for a given D J is precisely what obtained by the steepest descent method.

In the present model, the DM interaction is supposed in the plane, so in the GS the angle between in-plane NN is not zero. We show in Fig. 1 the relative orientation of the two NN spins in the plane.

Note that the perpendicular axes ηi and ηj coincide. Now, expressing the local frame of S j in the local frame of S i , we have ζj = cos θ i,j ζi + sin θ i,j ξi (10) ξj = -sin θ i,j ζi + cos θ i,j ξi (11) ηj = ηi (12) so that S j = S x j (cos θ i,j ζi -sin θ i,j ξi ) +S y j ηi + S z j (cos θ i,j ζi + sin θ i,j ξi )

The DM term of Eq. ( 4) can be rewritten as

S i ∧ S j = (-S z i S y j -S y i S x j sin θ i,j + S y i S z j cos θ i,j ) ξi +(S x i S x j sin θ i,j + S z i S z j sin θ i,j )η i +(S x i S y j -S y i S z j sin θ i,j -S y i S x j cos θ i,j ) ζi (14)
Using Eq. ( 7), we have

H DM = i,j D i,j • S i ∧ S j = D i,j
(S x i S x j e i,j sin θ i,j + S z i S z j e i,j sin θ i,j ) = D 4 i,j

[(S + i + S - i )(S + j + S - j )e i,j sin θ i,j

+4S z i S z j e i,j sin θ i,j ] (15)
where we have replaced S x = (S + + S -)/2. Note that e i,j sin θ i,j is always positive since for a NN on the positive axis direction, e i,j = 1 and sin θ i,j = sin θ where θ is positively defined, while for a NN on the negative axis direction, e i,j = -1 and sin θ i,j = sin(-θ) = -sin θ.

Note that for non collinear spin configurations, the local spin coordinates allow one to use the commutation relations between spin operators of a spin which are valid only when the z spin component is defined on its quantification axis. This method has been applied for helimagnets [28][START_REF] Diep | Low-temperature properties of quantum Heisenberg helimagnets[END_REF][START_REF] Quartu | Phase diagram of bodycentered tetragonal Helimagnets[END_REF].

III. SELF-CONSISTENT GREEN'S FUNCTION METHOD: FORMULATION

The GF method has been developed for non collinear surface spin configurations in thin films [START_REF] Ngo | Effects of frustrated surface in Heisenberg thin films[END_REF][START_REF] Ngo | Frustration effects in antiferrormagnetic face-centered cubic Heisenberg films[END_REF][START_REF] Diep | Quantum Theory of Helimagnetic Thin Films[END_REF][START_REF] Hog | Helimagnetic Thin Films: Surface Reconstruction, Surface Spin-Waves, Magnetization[END_REF]. Let us briefly recall here the principal steps of calculation and give the results for the present model. Expressing the Hamiltonian in the local coordinates, we obtain

H = - <i,j> J i,j 1 4 (cos θ i,j -1) S + i S + j + S - i S - j + 1 4 (cos θ i,j + 1) S + i S - j + S - i S + j + 1 2 sin θ i,j S + i + S - i S z j - 1 2 sin θ i,j S z i S + j + S - j + cos θ i,j S z i S z j + D 4 i,j [(S + i + S - i )(S + j + S - j )e i,j sin θ i,j +4S z i S z j e i,j sin θ i,j ] (16)
As said in the previous section, the spins lie in the xz planes, each on its quantization local z axis (Fig. 2). Note that unlike the sinus term of the DM Hamiltonian, Eq. ( 15), the sinus terms of H e , the 3rd line of Eq. ( 16), are zero when summed up on opposite NN (no e i,j to compensate). The 3rd line disappears therefore in the following.

At this stage it is very important to note that the standard commutation relations between spin operators S z and S ± are defined with z as the spin quantization axis. In non collinear spin configurations, calculations of SW spectrum using commutation relations without paying attention to this are wrong.

It is known that in two dimensions (2D) there is no long-range order at finite temperature (T ) for isotropic spin models with short-range interaction [34]. Thin films have small thickness, therefore to stabilize the ordering at finite T it is useful to add an anisotropic interaction. We use the following anisotropy between S i and S j which stabilizes the angle determined above between their local quantization axes S z i and S z j :

H a = - <i,j> I i,j S z i S z j cos θ i,j (17)
where I i,j is supposed to be positive, small compared to J , and limited to NN. Hereafter we take I i,j = I 1 for NN pair in the xz plane, for simplicity. As it turns out, this anisotropy helps stabilize the ordering at finite T in 2D as discussed. It helps also stabilize the SW spectrum at T = 0 in the case of thin films but it is not necessary for 2D and 3D at T = 0. The total Hamiltonian is finally given by

H = H e + H DM + H a (18)
We define the following two double-time GF's in the real space

G i,j (t, t ′ ) = << S + i (t); S - j (t ′ ) >> = -iθ(t -t ′ ) < S + i (t), S - j (t ′ ) > (19) F i,j (t, t ′ ) = << S - i (t); S - j (t ′ ) >> = -iθ(t -t ′ ) < S - i (t), S - j (t ′ ) > (20)
The equations of motion of these functions read

i dG i,j (t, t ′ ) dt = < S + i (t), S - j (t ′ ) > δ(t -t ′ ) -<< H, S + i ; S - j >> (21)
i dF i,j (t, t ′ ) dt = < S - i (t), S - j (t ′ ) > δ(t -t ′ ) -<< H, S - i ; S - j >> (22)
For the H e and H a parts, the above equations of motion generate terms such as << S z l S ± i ; S - j >> and << S ± l S ± i ; S - j >>. These functions can be approximated by using the Tyablikov decoupling to reduce to the above-defined G and F functions:

<< S z l S ± i ; S - j >>≃< S z l ><< S ± i ; S - j >> (23) << S ± l S ± i ; S - j >>≃< S ± l ><< S ± i ; S - j >>≃ 0 (24)
The last expression is due to the fact that transverse SW motions < S ± l > are zero with time. For the DM term, the commutation relations [H, S ± i ] give rise to the following term:

D l sin θ[∓S z i (S + l + S - l ) ± 2S ± i S z l ] (25)
which leads to the following type of GF's:

<< S z i S ± l ; S - j >>≃< S z i ><< S ± l ; S - j >> (26)
Note that we have replaced e i,j sin θ i,j by sin θ where θ is positive. The above equation is related to G and F functions [see Eq. ( 24)]. The Tyablikov decoupling scheme neglects higher-order functions.

We now introduce the following in-plane Fourier transforms g n,n ′ and f n,n ′ of the G and F Green's functions:

G i,j (t, t ′ , ω) = 1 ∆ BZ d k xz e -iω(t-t ′ ) ×g n,n ′ (ω, k xz )e i kxz.( Ri-Rj ) (27)
F i,j (t, t ′ , ω) = 1 ∆ BZ d k xz e -iω(t-t ′ ) ×f n,n ′ (ω, k xz )e i kxz .( Ri-Rj) (28)
where the integral is performed in the first xz Brillouin zone (BZ) of surface ∆, ω is the spin-wave frequency, n and n ′ are the indices of the layers along the c axis to which R i and R j belong (n = 1 being the surface layer, n = 2 the second layer and so on). We finally obtain the following matrix equation

M (E) h = u, (29)
where M (E) is a square matrix of dimension (2N × 2N ), h and u are the column matrices which are defined as follows

h =              g 1,n ′ f 1,n ′ . . . g n,n ′ f n,n ′ . . . g N,n ′ f N,n ′              , u =       2 S z 1 δ 1,n ′ 0 . . . 2 S z N δ N,n ′ 0       , (30)
where E = ω and M (E) is given by

           E + A 1 B 1 C 1 0 0 0 0 0 0 -B 1 E -A 1 0 -C 1 0 0 0 0 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 C n 0 E + A n B n C n 0 0 • • • 0 0 -C n -B n E -A n 0 -C n 0 • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 0 0 0 0 C N 0 E + A N B N 0 0 0 0 0 0 -C N -B N E -A N            (31)
with

A n = -J [8 < S z n > cos θ(1 + d n ) -4 < S z n > γ(cos θ + 1)] -2J ⊥ (< S z n-1 > + < S z n+1 >) -4D sin θ < S z n > γ +8D sin θ < S z n > (32) B n = 4J < S z n > γ(cos θ -1) -4D sin θ < S z n > γ (33) C n = 2J ⊥ < S z n > (34)
where n = 1, 2, ..., N , d n = I 1 /J , γ = (cos k x a + cos k z a)/2, k x and k z denote the wave-vector components in the xz planes, a the lattice constant. Note that (i) if n = 1 (surface layer) then there are no n -1 terms in the matrix coefficients, (ii) if n = N then there are no n + 1 terms. Besides, we have distinguished the in-plane NN interaction J from the inter-plane NN one J ⊥ .

In the case of a thin film, the SW eigenvalues at a given wave vector k = (k x , k z ) are calculated by diagonalizing the matrix .

The layer magnetization of the layer n is given by (see technical details in Ref. 36):

S z n = 1 2 - 1 ∆ dk x dk z 2N i=1 Q 2n-1 (E i ) e Ei/kB T -1 (35)
where n = 1, ..., N , and Q 2n-1 (E i ) is the determinant obtained by replacing the (2n -1)-th column of M by u at E i .

The layer magnetizations can be calculated at finite temperatures self-consistently using the above formula. The numerical method to carry out this task has been described in details in Refs. [START_REF] Diep | Quantum Theory of Helimagnetic Thin Films[END_REF]. One can summarize here: (i) using a set of trial values (inputs) for S z n (n = 1, ..., N ), one diagonalizes the matrix to find spin-wave energies E i which are used to calculate the outputs S z n (n = 1, ..., N ) by using Eq. ( 35), (ii) using the outputs as inputs to iterate the equations, (iii) if the output values are the same as the inputs within a precision (usually at 0.001%), the iteration is stopped. The method is thus self-consistent.

The value of the spin in the layer n at T = 0 is calculated by [START_REF] Diep | Quantum Theory of Helimagnetic Thin Films[END_REF][START_REF] Diep | Theory of Magnetism: Application to Surface Physics[END_REF]

S z n (T = 0) = 1 2 + 1 ∆ dk x dk z N i=1 Q 2n-1 (E i ) (36)
where the sum is performed over N negative values of E i (for positive values, the Bose-Einstein factor in Eq. ( 35) is equal to 0 at T = 0). The transition temperature T c can be calculated by letting S z n on the left-hand side of Eq. ( 35) to go to zero. The energy E i tends then to zero, so that we can make an expansion of the exponential at T = T c . We have

1 k B T c = 2 ∆ dk x dk z 2N i=1 Q 2n-1 (E i ) E i (37)

IV. TWO AND THREE DIMENSIONS: SPIN-WAVE SPECTRUM AND MAGNETIZATION

Consider just one single xz plane. The above matrix is reduced to two coupled equations

(E + A n )g n,n ′ + B n f n,n ′ = 2 < S z n > δ(n, n ′ ) -B n g n,n ′ + (E -A n )f n,n ′ = 0 (38)
where A n is given by ( 32) but without J ⊥ term for the 2D case considered here. Coefficients B n and C n are given by ( 33) and (34) with C n = 0. The poles of the GF are the eigenvalues of the SW spectrum which are given by the secular equation

(E + A n )(E -A n ) + B 2 n = 0 [E + A n ][E -A n ] + B 2 n = 0 E 2 -A 2 n + B 2 n = 0 E = ± (A n + B n )(A n -B n ) ( 39
)
where ± indicate the left and right SW precessions. Several remarks are in order: (i) if θ = 0, we have B n = 0 and the last three terms of A n are zero. We recover then the ferromagnetic SW dispersion relation

E = 2ZJ < S z n > (1 -γ) (40)
where Z = 4 is the coordination number of the square lattice (taking

d n = 0), (ii) if θ = π, we have A n = 8J < S z n >, B n = -8J < S z n > γ.
We recover then the antiferromagnetic SW dispersion relation

E = 2ZJ < S z n > 1 -γ 2 (41)
(iii) in the presence of a DM interaction, we have 0 < cos θ < 1 (0 < θ < π/2). If d n = 0, the quantity in the square root of Eq. ( 39) is always ≥ 0 for any θ. It is zero at γ = 1. The SW spectrum is therefore stable at the long-wavelength limit. The anisotropy d n gives a gap at γ = 1.

As said earlier, the necessity to include an anisotropy has a double purpose: it permits a gap and stabilizes a long-range ordering at finite T in 2D systems.

Figure 3 shows the SW spectrum calculated from Eq. (39) for θ = 30 degrees (π/6 radian) and 80 degrees (1.396 radian). The spectrum is symmetric for positive and negative wave vectors and for left and right precessions. Note that for small θ (i. e. small D) E is proportional to k 2 at low k (cf. Fig. 3a), as in ferromagnets. However, as θ increases, we observe that E becomes linear in k as seen in 3b. This is similar to antiferromagnets. The change of behavior is progressive with increasing θ, we do not observe a sudden transition from k 2 to k behavior. This feature is also observed in three dimensions (3D) and in thin films as seen below.

It is noted that, thanks to the existence of the anisotropy d, we avoid the logarithmic divergence at k = 0 so that we can observe a long-range ordering at finite T in 2D. We show in Fig. 4 the magnetization M (≡< S z >) calculated by Eq. ( 35) for one layer using d = 0.001. It is interesting to observe that M depends strongly on θ: at high T , larger θ yields stronger M . However, at T = 0 the spin length is smaller for larger θ due to the so-called spin contraction [START_REF] Diep | Theory of Magnetism: Application to Surface Physics[END_REF] calculated by Eq. ( 36). As a consequence there is a cross-over of magnetizations with different θ at low T as shown in Fig. 4.

Let us study the 3D case. The crystal is periodic in three directions. We can use the Fourier transformation in the y direction, namely g n±1 = g n e ±iky a and f n±1 = f n e ±ikya . The matrix [START_REF] Quartu | Phase diagram of bodycentered tetragonal Helimagnets[END_REF] is reduced to two coupled equations of g and f functions, omitting index n,

(E + A ′ )g + Bf = 2 < S z > -Bg + (E -A ′ )f = 0 (42)
where

A ′ = -J [8 < S z > cos θ(1 + d) -4 < S z > γ(cos θ + 1)] +4J ⊥ < S z > cos(k y a) -4D sin θ < S z > γ +8D sin θ < S z > (43) B = 4J < S z > γ(cos θ -1) -4D sin θ < S z > γ ( 44
)
The spectrum is given by Unlike the 2D case where the angle is inside the plane so that the antiferromagnetic case can be recovered by setting cos θ = -1 as seen above, one cannot use the above formula to find the antiferromagnetic case because in the 3D formulation it was supposed a ferromagnetic coupling between planes, namely there is no angle between adjacent planes in the above formulation.

E = ± (A ′ + B)(A ′ -B) (45
The same consideration as in the 2D case treated above shows that for d = 0 the spectrum E ≥ 0 for positive precession and E ≤ 0 for negative precession, for any θ.

The limit E = 0 is at γ = 1 ( k = 0). Thus there is no instability due to the DM interaction. Using Eq. (45), we have calculated the 3D spectrum. This is shown in Fig. 5 for a small and a large value of θ. As in the 2D case, we observe E ∝ k when k → 0 for large θ. Main properties of the system are dominated by the in-plane DM behavior.

Figure 6a displays the magnetization M versus T for several values of θ. As in the 2D case, when θ is not zero, the spins have a contraction at T = 0: a stronger θ yields a stronger contraction. This generates a magnetization cross-over at low T shown in the inset of Fig. 6a. The spin length at T = 0 versus θ is displayed in Fig. 6b. Note that the spin contraction in 3D is smaller than that in 2D. This is expected since quantum fluctuations are stronger at lower dimensions.

V. THE CASE OF A THIN FILM: SPIN-WAVE SPECTRUM, LAYER MAGNETIZATIONS

In the 2D and 3D cases shown above, there is no need at T = 0 to use a small anisotropy d. However in the case of thin films shown below, due to the lack of neighbors at the surface, the introduction of a DM interaction the other half negative (left and right precessions): E i (i = 1, ..., 2N ). Figure 8 shows the case of a film of 8 layers with J = J ⊥ = 1 for a weak and a strong value of D (small and large θ). As in the 2D and 2D cases, for strong D, E is proportional to k at small k (cf. Fig. 8b). It is noted that this behavior concerns only the first mode. The upper modes remain in the k 2 behavior. Figure 9 shows the layer magnetizations of the first four layers in a 8-layer film (the other half is symmetric) for several values of θ. In each case, we see that the surface layer magnetization is smallest. This is a general effect of the lack of neighbors for surface spins even when there is no surface-localized SW as in the present simple-cubic lattice case [START_REF] Diep | Theory of Magnetism: Application to Surface Physics[END_REF]. The spin length at T = 0 for a 8-layer film is shown in Fig. 10 as a function of θ. One observes that the spins are strongly contracted with large θ. Let us touch upon the surface effect in the present model. We know that for the simple cubic lattice, if the interactions are the same everywhere in the film, then there is no surface localized modes, and this is true with DM interaction (see spectrum in Fig. 8) and without DM interaction (see Ref. 37). In order to create surface modes, we have to take the surface exchange interactions different from the bulk ones. Low-lying branches of surface modes which are "detached" from the bulk spectrum are seen in the SW spectrum shown in Fig. 11a with J s = 0.5, J s ⊥ = 0.5. These surface modes strongly affect the surface magnetization as observed in Fig. 11b: the surface magnetization is strongly diminished with increasing T . The role of surface-localized modes on the strong decrease of the surface magnetization as T increases has already been analyzed more than 30 years ago [START_REF] Diep-The-Hung | Effect of surface spin-waves and surface anisotropy in magnetic thin films at finite temperatures[END_REF]. We show now the effect of the film thickness in the present model. The case of thickness N = 12 is shown in Fig. 12a with θ = π/6 where the layer magnetizations versus T are shown in details. The gap at k = 0 due to d is shown in Fig. 12b as a function of the film thickness N for d = 0.1 and θ = π/6, at T = 0. We see that the gap depends not only on d but also on the value of the surface magnetization which is larger for thicker films. The transition temperature T c versus the thickness N is shown in Fig. 12c where one observes that T c tends rapidly to the bulk value (3D) which is ≃ 2.82 for d = 0.1.

VI. CONCLUDING REMARKS

By a self-consistent Green's function theory, we obtain the expression of the spin-wave dispersion relation in 2D and 3D as well as in a thin film. Due to the competition between ferromagnetic interaction J and the perpendicular DM interaction D, the GS is non linear with an angle θ which is shown to explicitly depend on the ratio D/J. The spectrum is shown to depend on θ and the layer magnetization is calculated self-consistently as a function of temperature up to the critical temperature T c .

We have obtained new and interesting results. In particular we have showed that (i) the spin-wave excitation in 2D and 3D crystals is stable at T = 0 with the non collinear spin configuration induced by the DM interaction D without the need of an anisotropy, (ii) in the case of thin films, we need a small anisotropy d to stabilize the spin-wave excitations because of the lack of neighbors at the surface, (iii) the spin-wave energy E depends on D, namely on θ: at the long wave-length limit, E is proportional to k 2 for small D but E is linear in k for strong D, in 2D and 3D as well as in a thin film, (iv) quantum fluctuations are inhomogeneous for layer magnetizations near the surface, (v) unlike in some previous works, spin waves in systems with asymmetric DM interactions are found to be symmetric with respect to opposite propagation directions.

FIG. 1 .

 1 FIG. 1. (a) The ground state is a planar configuration on the xz plane. The figure shows the case where θ = π/6 (D = -0.577), J = J ⊥ = 1 using the steepest descent method ; (b) a zoom is shown around a spin with its nearest neighbors.

FIG. 2 .

 2 FIG. 2. Local coordinates in the xz plane. The spin quantization axes of Si and Sj are ζi and ζj , respectively.

FIG. 3 .FIG. 4 .

 34 FIG. 3. Spin-wave spectrum E(k) versus k ≡ kx = kz for (a) θ = 0.524 radian and (b) θ = 1.393 in two dimensions at T = 0.1. Positive and negative branches correspond to right and left precessions. A small d (= 0.001) has been used to stabilized the ordering at finite T in 2D. See text for comments.

)

  If cos θ = 1 (ferromagnetic), one has B = 0. By regrouping the Fourier transforms in three directions, one obtains the 3D ferromagnetic dispersion relation E = 2Z < S z > (1 -γ 2 ) where γ = [cos(k x a)+ cos(k y a)+ cos(k z a)]/3 and Z = 6, coordination number of the simple cubic lattice.

FIG. 5 .FIG. 6 .FIG. 7 .

 567 FIG. 5. Spin-wave spectrum E(k) versus k ≡ kx = kz for θ = π/6 (red circles) and θ = π/3 (blue circles) in three dimensions at T = 0.1, with d = 0. Note the linear-k behavior at low k for the large value of θ (inset). See text for comments.

FIG. 8 .

 8 FIG. 8. Spin-wave spectrum E(k) versus k ≡ kx = kz for a thin film of 8 layers: (a) θ = π/6 (in radian) (b) θ = π/3, using d = dc for each case (dc = 0.012 and 0.021, respectively). Positive and negative branches correspond to right and left precessions. Note the linear-k behavior at low k for the large θ case. See text for comments.

FIG. 9 . 8 -

 98 FIG. 9. 8-layer film: layer magnetizations M versus temperature T for (a) θ = π/6 (radian), (b) θ = π/3, with d = 0.1. Red circles, blue void circles, green void triangles and magenta squares correspond respectively to the first, second, third and fourth layer.

FIG. 10 .

 10 FIG. 10. Spin length S0 at T = 0 of the first 4 layers as a function of θ, for N = 8, d = 0.1. Red circles, blue void circles, green void triangles and magenta squares correspond respectively to the first, second, third and fourth layer.

FIG. 11 .

 11 FIG. 11. Surface effect: (a) spin-wave spectrum E(k) versus k = kx = kz for a thin film of 8 layers: θ = π/6, d = 0.2, J s = 0.5, J s ⊥ = 0.5, the gap at k = 0 is due to d. The surface-mode branches are detached from the bulk spectrum. (b) Layer magnetizations versus T for the first, second, third and fourth layer (red circles, green void circles, blue void circles and magenta filled squares, respectively). See text for comments.

FIG. 12 .

 12 FIG.12. 12-layer film: (a) Layer magnetizations versus T for θ = π/6 and anisotropy d = 0.1. Red circles, blue squares, green void squares magenta circles, void turquoise triangles and brown triangles correspond respectively to first, second, third, fourth, fifth and sixth layer, (b) Gap at k = 0 as a function of film thickness N for θ = π/6, d = 0.1, at T = 0.1, (c) Critical temperature Tc versus the film thickness N calculated with θ = π/6 and d = 0.1 using Eq.[START_REF] Diep-The-Hung | Effect of surface spin-waves and surface anisotropy in magnetic thin films at finite temperatures[END_REF]. Note that for infinite thickness (namely 3D), Tc ≃ 2.8 for d = 0.1.

ACKNOWLEDGMENTS

This study is a part of a project financed by Narodowe Centrum Nauki (National Science Center of Poland), Grant No. DEC-2013/08/M/ST3/00967. SEH acknowledges a financial support from Agence Universitaire de la Francophonie (AUF).