Kaouther Moussa
email: kaouther.moussa@uphf.fr

Mirko Fiacchini
email: mirko.fiacchini@grenoble-inp.fr

Mazen Alamir
email: mazen.alamir@grenoble-inp.fr

Probabilistically Certified Region of Attraction of a Tumor Growth Model with Combined Chemo-and Immunotherapy

The recent progress in immunology lead to a considerable interest in modeling cancer dynamics in order to better understand and analyze such complex systems. Many works have been carried out in order to design cancer treatment protocols using mathematical models. One of the main complexities of such models is the presence of different types of uncertainties, which remains less considered in the literature. This paper deals with the estimation of regions of attraction (RoAs) under parametric uncertainties for a cancer growth model with combined therapies. We propose a framework of probabilistic certification, based on the randomized methods, in order to derive probabilistically certified RoAs of a cancer growth model. The model considered in this paper describes the interaction between a tumor and the immune system in presence of a combined chemo-and immunotherapy treatment, with considerations on pharmacokinetics and pharmacodynamics (PK/PD) of both treatments.

The last decades witnessed a considerable progress in experimental and clinical immunology as well as in modeling the immune system dynamics [START_REF] Eftimie | Mathematical models for immunology: current state of the art and future research directions[END_REF]. The progress in cancer dynamics modeling motivated researchers to apply control approaches in order to schedule cancer treatments using optimal control strategies. We can find in the literature many works regarding the application of optimal control approaches on cancer treatment problems. For instance [START_REF] Onofrio | A bi-parametric model for the tumour angiogenesis and antiangiogenesis therapy[END_REF], where optimal protocols for anti-angiogenic therapy were investigated, or [START_REF] De Pillis | Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls[END_REF] where linear controls were designed for a tumor-immune interactions model with chemotherapy delivery. However, only few works addressed the problem of handling parametric uncertainties. One can cite for example, [START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF] where a robust feedback scheme is proposed to schedule anti-angiogenic treatment combined with chemotherapy, [START_REF] Kovacs | Model-based angiogenic inhibition of tumor growth using modern robust control method[END_REF] where an H ∞ -based robust control was applied to the same model and [START_REF] Alamir | On probabilistic certification of combined cancer therapies using strongly uncertain models[END_REF] where a general framework for probabilistic certification of cancer therapies was proposed.

The estimation of the region of attraction for cancer models is an interesting problem since it provides a set of possible initial conditions (tumor volume and immune density for example) that can be driven to a desired target set (benign region). This problem becomes complex when dealing with nonlinear systems and even more challenging for uncertain systems. There are some works which dealt with the problem of estimating the RoA for cancer models, we cite for example [START_REF] Doban | Domain of attraction computation for tumor dynamics[END_REF] and [START_REF] Zarei | Perturbed tumor immunotherapy domain of attraction estimation via the arc-length function[END_REF], where the authors proposed different Lyapunov functions based approaches, to estimate the domain of attraction of the tumor free equilibrium point corresponding to autonomous cancer growth models, where no therapies are considered, see also [START_REF] Merola | An insight into tumor dormancy equilibrium via the analysis of its domain of attraction[END_REF] and references therein. However, only few works considered model uncertainties, in particular in [START_REF] Riah | Iterative method for estimating the robust domains of attraction of non-linear systems: Application to cancer chemotherapy model with parametric uncertainties[END_REF], an iterative method to estimate the robust RoA was presented. Furthermore, the work in [START_REF] Doban | A switching control law approach for cancer immunotherapy of an evolutionary tumor growth model[END_REF] dealt with the estimation of regions of attraction for cancer models, using rational Lyapunov functions and switching control strategies. The latter work has been further investigated in [START_REF] Doban | Stability domains computationand stabilization of nonlinear systems: implications for biological systems[END_REF] in order to take into account model disturbances, based on input to state (ISS) Lyapunov functions. Moreover, the model proposed in [START_REF] Doban | A switching control law approach for cancer immunotherapy of an evolutionary tumor growth model[END_REF] has been used for data fitting in [START_REF] Guillerey | Chemotherapy followed by anti-CD137 mAb immunotherapy improves disease control in a mouse myeloma model[END_REF] for real-life experiments. Finally, it is important to highlight the fact that the estimation of robust RoAs is based on the worst-case scenario analysis leading to potentially pessimistic design, this because the worst-case is considered no matter how small its probability of occurrence is.

The presence of parametric uncertainties can drastically affect the efficiency of a nominal controller as well as the size of the estimated RoAs. In this work, we propose a framework to probabilistically certify the existence of a control structure that drives the states corresponding to tumor cells and immune density from an initial state set to a certified target set. This probabilistic certification framework is based on the randomized methods proposed in [START_REF] Alamo | Randomized strategies for probabilistic solutions of uncertain feasibility and optimization problems[END_REF] and [START_REF] Alamo | Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms[END_REF], which, unlike the robust classical design, avoids focusing on few unlikely extremely bad scenarios allowing to overcome the conservatism of the robust RoA design. The methodology that we propose in this paper consists mainly of two steps. Firstly, we derive an ordered sequence of sets and their control strategy such that the states can be driven from a set to the previous one with a certain probabilistic guarantee. The appropriate choice of the first set allows to insure that the union of the sets is a probabilistically certified approximation of the RoA. The second step consists in providing a global certification on the probability of convergence to the initial certified target set. The randomized methods have been already used in [START_REF] Hokayem | Nonlinear systems stability via random and quasi-random methods[END_REF] in order to determine the stability region for nonlinear deterministic systems, without taking into consideration parametric uncertainties.

The model that we investigate here is a modified version of the classical Stepanova one [START_REF] Stepanova | Course of the immune reaction during the development of a malignant tumour[END_REF] that has been extensively used in the literature, we cite for example [START_REF] Onofrio | On the Dynamics of Tumor-Immune System Interactions and Combined Chemo-and Immunotherapy[END_REF], [START_REF] Ledzewicz | On optimal protocols for combinations of chemoand immunotherapy[END_REF] and [START_REF] Ledzewicz | On the role of the objective in the optimization of compartmental models for biomedical therapies[END_REF] where optimal control methodologies were proposed to schedule chemo-and immunotherapy administration profiles. Furthermore, [START_REF] Sharifi | Multiple model predictive control for optimal drug administration of mixed immunotherapy and chemotherapy of tumours[END_REF] proposed a model predictive control scheme to design chemo-and immunotherapy administration schedules. In [START_REF] Sharifi | Overcoming channel uncertainties in touchable molecular communication for direct drug targeting assisted immunochemotherapy[END_REF], the authors proposed a robust model predictive control scheme, in order to consider direct drug targeting pharmacokinetic uncertainties as well as system model mismatches. In this paper, we model the concentration of the chemotherapy agent in the plasma and the tumor site via a pharmacokinetics compartmental model, we also model the pharmacokinetics of immunotherapy, as well as the pharmacodynamics of both drugs. Although the classical Stepanova model has been widely used, it has never been investigated in the literature to estimate the controlled region of attraction of its corresponding tumor free equilibrium. Therefore, we aim at pointing out the importance of uncertainties considerations in the RoA estimation for such models. This paper is organized as follows: In Section 1, the dynamical cancer model and the problem of RoA probabilistic certification are introduced. Section 2 recalls the randomized algorithms approach for probabilistic certification. In Section 3, a framework for RoA probabilistic certification is proposed, based on the randomized methods presented in [START_REF] Alamo | Randomized strategies for probabilistic solutions of uncertain feasibility and optimization problems[END_REF] and [START_REF] Alamo | Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms[END_REF]. In Section 4, the proposed RoA probabilistic certification framework is applied to the considered cancer model. Finally, Section 5 summarizes the contribution that we present in this paper.

Problem statement

The following nonlinear dynamical system describes the interaction between a tumor and the immune system in presence of chemotherapy and immunotherapy treatments:

ẋ1 = µ C 1 -x 1 x∞ x 1 -γ X x 1 x 2 -κ X x γc 4 x γc 4 +C γc 50c x 1 ẋ2 = µ I (1 -βx 1 ) x 1 x 2 + α Y -δ Y x 2 -η Y x 2 x 3 + κ Y x γ i 5 x γ i 5 +C γ i 50i x 2 ẋ3 = -(k 1 + k 2 ) x 3 + s 1 u 1 (t) V 1 ẋ4 = k 12 V 1 V 2 x 3 -k 2 x 4 ẋ5 = -c i x 5 + s 2 u 2 (t), x(0) = (x 1 (0), x 2 (0), x 3 (0), x 4 (0), x 5 (0)) = x 0 , (1)
where the variables are defined as follows:

x 1 The number of tumor cells (10 6 •cells)

x 2 The density of effector immune cells (ECs) (dimensionless)

x 3 The concentration of chemotherapy in the plasma (µg•mL -1 )

x 4 The concentration of chemotherapy in the tumor effect site (µg•mL -1 )

x 5 The concentration of immunotherapy in the immune cells site (mg•mL This model is an extension of the model presented in [START_REF] Onofrio | On the Dynamics of Tumor-Immune System Interactions and Combined Chemo-and Immunotherapy[END_REF] that has the advantage of being a low dimensional system that nevertheless includes the main aspects of cancerimmune interactions, and it has been widely used in the literature for cancer drug scheduling. In many models it is assumed that the drug concentration is equal to its dosage which is an oversimplification [START_REF] Ledzewicz | Optimal controls for a model with pharmacokinetics maximizing bone marrow in cancer chemotherapy[END_REF]. Therefore, we revisited the model proposed in [START_REF] Onofrio | On the Dynamics of Tumor-Immune System Interactions and Combined Chemo-and Immunotherapy[END_REF] by adding a pharmacokinetic (PK) compartment (involving x 3 and x 4 ) that allows to model the concentration of chemotherapy in the plasma and the tumor effect site. Furthermore, we incorporated to this model the pharmacodynamics of chemotherapy using a Hill function, the equations and the parameters values have been taken from [START_REF] Bukkuri | Optimal control analysis of combined chemotherapy-immunotherapy treatment regimens in a PKPD cancer evolution model[END_REF] and [START_REF] Iliadis | Optimizing drug regimens in cancer chemotherapy by an efficacy-toxicity mathematical model[END_REF] for the Etopside drug. Other types of chemotherapy can be considered such as Pegylated Liposomal Doxorubicin, whose main corresponding medical parameters have been presented in [START_REF] Drexler | Modeling of tumor growth incorporating the effect of pegylated liposomal doxorubicin[END_REF]. Moreover, we extended the model with the pharmacokinetics of immunotherapy which contains only one compartment represented by x 5 which is the concentration of immunotherapy in the immune cells site, since we consider that the immunotherapy does not have a direct inhibition effect on the tumor. The pharmacodynamics effects have also been incorporated using a Hill function. The parameters values related to the PK/PD of immunotherapy have been taken from [START_REF] Ionescu | A minimal PKPD interaction model for evaluating synergy effects of combined NSCLC therapies[END_REF].

Table 1 summarizes the definitions of the model parameters and their nominal values. Furthermore, the parameters s 1 and s 2 allow to scale the drugs effects. Note that we consider examples of possible treatment protocols, nevertheless, it is worth emphasizing that in this paper, we focus on the assessment of a methodology that remains applicable for different parameters values, treatment strategies and equations.

Let's denote by x = (x 1 , x 2 , x 3 , x 4 , x 5 ) and u = (u 1 , u 2 ) respectively, the state and the control input vectors. In this paper, we consider a cycle-based treatment, where the drugs are injected following N C therapeutic cycles. Each cycle having two phases, a hospitalization period lasting 5 days, where the patient receives one infusion per day, and a rest period where the patient recovers. Figure 1 shows a typical temporal combined control structure, the different notations in this figure are defined as follows: Concentration of chemotherapy infusion T = 5

Hospitalization duration (days)

T c = 15
Cycle duration (days)

For a given treatment cycle, the therapeutic profile considered in this paper is completely defined by the following control parametrization θ:

θ = (ν C , σ C , d C , σ I , d I ) . (2)
In cancer treatment design, we usually have some constraints to satisfy, they can be defined either on the states or on the control inputs. These constraints allow to prevent from drug toxicity and excessive immune weakening. In this paper, we consider the following constraints for all t ∈ [0, T ], with T ∈ R + :

x 2 (t) ≥ c, with c ∈ R + , (3)
0 ≤ u 1 (t) ≤ 1, (4)
0 ≤ u 2 (t) ≤ 1, (5)
where ( 3) is a health constraint on the minimal density of immune cells. The constraints on u 1 (t) and u 2 (t) for all t, see ( 4) and ( 5), are normalized drug toxicity constraints, they can be satisfied by properly choosing the parametrization θ (namely, d C and d I ) of the control input u. Therefore, we will consider only the constraint (3), since the satisfaction of the other constraints can be monitored by a proper choice of θ.

The uncontrolled model (1) (for u = (0, 0)) has two locally asymptotically stable equilibriums points. The macroscopic malignant equilibrium is x m = (766.4, 0.018, 0, 0, 0) and the benign one is x b = (41.45, 0.954, 0, 0, 0). In general, the objective of the treatment is to drive the state initial conditions to the region of attraction of the benign equilibrium (safe region), without constraints violation. We are specifically interested in characterizing the set of initial conditions (tumor volume and immune density) from which the trajectories of (1) can be driven to the safe region under parametric uncertainties.

In [START_REF] Moussa | Robust domain of attraction estimation for a tumor growth model[END_REF], we proposed a methodology to characterize the controlled region of attraction of model ( 1) with bang-bang controls (without pharmacokinetics). Then, we used this approach to derive and estimate the robust region of attraction. In this paper, we propose to derive a probabilistically certified RoA for model [START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF], that is based on chance-constrained problems, tolerating some constraints violations provided that their corresponding probability is small enough. Definition 1 We denote by Ω 0 ⊆ R 5 + a set of initial conditions x 0 , such that the state trajectories (x 1 , x 2 ) corresponding to model [START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF], when no control is applied (i.e. u = (0, 0)), converge to their corresponding benign equilibrium in spite of all uncertainties realizations, with a given confidence probability. This set can be seen as a targeted safe set for each treatment. Note that since u = (0, 0) and the initial conditions x 3 (0), x 4 (0), x 5 (0) are all equal to zero, the set Ω 0 is technically determined in R 2 + . The region of uncertainties is defined such that the parameters are normally distributed in the interval [0.9p nom , 1.1p nom ], where p nom stands for the vector of nominal parameters that will be properly defined in Section 4.

Problem 1 (Estimation of a probabilistically certified RoA) We aim at computing a sequence of sets {Ω k } N C k=1 , for N C therapeutic cycles. Those sets are determined in the space of (x 1 , x 2 ), representing respectively, the cancer burden (defined by the number of cancer cells) and the ECs density, such that, in the family of control parametrizations that we consider, there exists a therapeutic protocol that drives, with a desired probability, the states from Ω k+1 to k j=0 Ω j without safety constraints violations.

We denote by Ωpnom 0 an estimation of the region of attraction of the benign equilibrium for u = (0, 0), when nominal model parameters (in Table 1) are considered.

Overview on probabilistic certification

The randomized algorithms were presented in [START_REF] Tempo | Randomized algorithms for analysis and control of uncertain systems, with applications[END_REF], [START_REF] Alamo | Randomized strategies for probabilistic solutions of uncertain feasibility and optimization problems[END_REF] and [START_REF] Alamo | Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms[END_REF] in order to solve optimization problems with probabilistic constraints satisfaction. In contrast to the standard robust control design, which is based on the worst-case scenario analysis leading hence to potentially pessimistic design, the randomized methods provide the possibility to avoid focusing on the worst scenarios if their probability of occurrence is small. Therefore, this framework is very interesting from the cancer treatment point of view, since the latter involves many uncertainties that have to be considered.

This section aims at briefly recalling the main key-points of the randomized methods that are important for the assessment of the approach that we propose in this paper, in which we present a framework of estimation of probabilistically certified region of attraction for a cancer therapies dynamical model.

Consider the following optimization problem :

min θ∈Θ J(θ) s.t. ∀p g c (θ, p) = 0, (6)
where θ ∈ Θ ⊂ R n θ is the decision variable (which can be a parametrization of a control law) and p is the uncertainties vector following the probability measure P defined in the set P (the vector p can contain for example model parameters that are considered to be uncertain), J is the cost to be minimized. In terms of control design for dynamical systems, the cost J can involve the states, the input variables, their respective integrals with respect to time or any combination of these indicators. Finally, g c is an indicator function on the violation of some given constraints and is defined as follows: g c (θ, p) := 0 if all the constraints are satisfied 1 otherwise

The randomized method consists in replacing the original hard problem in ( 6) by the following relaxed problem:

min θ∈Θ J(θ) s.t. Pr P {g c (θ, p) = 1} ≤ η, (7)
where the constraint is on the probability of constraints violation, giving therefore a soft constraint in the sense that we can accept a value of θ which minimizes the cost J, even if the constraints are violated for some realizations of p, provided that the probability of these realizations is less than or equal to η (small enough). Even though the constraint in [START_REF] Doban | Stability domains computationand stabilization of nonlinear systems: implications for biological systems[END_REF] simplifies the previous constraint in [START_REF] De Pillis | Chemotherapy for tumors: An analysis of the dynamics and a study of quadratic and linear optimal controls[END_REF], the computation of the violation probability remains expensive. Authors in [START_REF] Alamo | Randomized strategies for probabilistic solutions of uncertain feasibility and optimization problems[END_REF] and [START_REF] Alamo | Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms[END_REF] proposed a simplification which consists in replacing the probability by the mean value of g c over N p drawn independent identically distributed (i.i.d.) samples of p in P according to the probability distribution P. Therefore, the simplified optimization problem is the following:

min θ∈Θ J(θ) s.t. Np i=1 g c θ, p (i) N p ≤ m N p , (8)
where m is the number of allowed constraints violations. In [START_REF] Alamo | Randomized strategies for probabilistic solutions of uncertain feasibility and optimization problems[END_REF] and [START_REF] Alamo | Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms[END_REF], several bounds on N p are given such that the fulfillment of the constraint in [START_REF] Doban | Domain of attraction computation for tumor dynamics[END_REF] implies that the probability condition in ( 7) is satisfied with a confidence probability greater than or equal to 1δ.

The bounds that are derived on N p involve the precision η and the confidence of fulfillment δ.

In this paper, we are interested in specific control structures, since cancer treatment schedules are often defined by cycles with a hospitalization period where the patient receives several drug infusions and a rest period for recovery. Therefore, it is more adequate in this case to consider that the controls are parametrized by a discrete variable θ with cardinality n Θ ∈ N. This choice of θ simplifies the optimization problem [START_REF] Doban | Stability domains computationand stabilization of nonlinear systems: implications for biological systems[END_REF], since it can be solved by a simple enumeration. In this case, the following proposition from [START_REF] Alamo | Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms[END_REF] holds:

Proposition 1 Let m ∈ N be any integer representing the number of accepted failures.

Let δ ∈ (0, 1) be a targeted confidence parameter. Take N p satisfying

N p ≥ 1 η m + ln n Θ δ + 2m ln n Θ δ 1 2 (9)
then any solution of [START_REF] Doban | Domain of attraction computation for tumor dynamics[END_REF] in which {p (j) } N j=1 are i.i.d. following the probability distribution P satisfies the constraint in [START_REF] Doban | Stability domains computationand stabilization of nonlinear systems: implications for biological systems[END_REF] with a probability greater than or equal to 1δ The inequality ( 9) is mathematically based on the binomial distribution. In this section, we presented a concise overview of the basic theoretical aspects of this methodology, the readers interested in further details might refer to [START_REF] Alamo | Randomized methods for design of uncertain systems: Sample complexity and sequential algorithms[END_REF].

It is interesting to notice that the bound on N p provided by Proposition 1 does not depend on the dimension of p which is useful when having many uncertain parameters in the certification problem. Furthermore, as we can see in Table 3, since the confidence parameter δ affects the bound logarithmically, we can have a highly confident certification with a tractable number of random samples.

Table 3: The evolution of the number of samples N p required to achieve the certification, with respect to the confidence design parameter δ and the number of control parametrizations n Θ , for η = 10 -2 and m = 1. Furthermore, for a specific desired confidence parameter δ = 10 -3 , Table 4 provides an idea on the evolution of the number of trials N p that should be performed for each possible control law θ, with respect to the precision parameter η and the number of control parametrizations n Θ . Therefore, the total number of simulations is

n Θ δ = 0.1 δ = 0.01 δ = 0.
N sim = N p • n Θ .
Table 4: The evolution of the number of samples N p required to achieve the certification, with respect to the precision design parameter η and the number of control parametrizations n Θ , for δ = 10 -3 and m = 1. This approach provides a powerful pragmatic tool allowing to certify control strategies. In [START_REF] Alamir | On probabilistic certification of combined cancer therapies using strongly uncertain models[END_REF], a randomized methods based framework for probabilistic certification of feedback control strategies has been proposed for a combined cancer therapy model.

n Θ η = 0.1 η = 0.01 η = 0.

Probabilistic certification of ROA

In this section, we will establish a framework of RoA probabilistic certification, based on the randomized methods presented in the previous section. We propose to use this general framework in order to probabilistically certify the existence of a control structure which allows to drive initial states from a given set to a target set under parametric uncertainties. Let us rewrite system (1) into the following form:

ẋ = F (x, u, p), ( 10
)
where p is the vector of parameters that model (1) involves. Furthermore, we consider that the variables of system (10) are subject to the following constraints:

x ∈ X, x(T ) ∈ Ω, u ∈ U, (11)
where the sets X and U are defined as follows:

X = {x ∈ R 5 + | x 2 ≥ c}, (12)
U = {(u 1 , u 2 ) ∈ R 2 + | 0 ≤ u 1 , u 2 ≤ 1}, (13)
and Ω ⊆ R 5 + represents the target set that will certified to be safe for each treatment cycle.

Remark 1 Note that since the states x 3 , x 4 and x 5 vanish after the treatment period, when u = (0, 0), the certified safe set Ω will be determined technically in the space (x 1 , x 2 ) which is a subset of R 2 + .

As previously mentioned, we consider that the control inputs are parametrized by a vector θ which lies in a discrete set Θ with cardinality n Θ ∈ N. This choice of θ fits particularly to the case of cancer therapy design, since some of the parameters involved in the treatment scheduling are naturally quantified. Suppose that the parameters vector p is a random variable following the probability distribution P that we denote p ∼ P. Given a set Γ ⊆ R 5 + , containing initial conditions x 0 (Remark (1) holds for Γ since the initial conditions x 3 (0), x 4 (0), x 5 (0) are null), and a parameterization of the input θ ∈ Θ, let's consider the following optimization problem: min θ∈Θ J(θ) s.t. ∀ (x 0 , p) ∈ (Γ × P) g c (θ, x 0 , p, Ω) = 0, [START_REF] Guillerey | Chemotherapy followed by anti-CD137 mAb immunotherapy improves disease control in a mouse myeloma model[END_REF] where J(θ) is a cost function to be minimized. In terms of cancer treatment design, this function can be a combination of many objectives that one seeks to achieve, for example reducing the quantity of injected drugs, to prevent from toxicity, or reducing the duty cycle in order to reduce the hospitalization duration. g c is the failure indicator function, defined on the state trajectories of [START_REF] Onofrio | A bi-parametric model for the tumour angiogenesis and antiangiogenesis therapy[END_REF]. Note that g c depends on Ω in order to emphasize the fact that the terminal constraint, expressing that the health indicators belong to the certified target set, is a part of the failure indicator that will be properly defined in the sequel. The function g c is deterministic such that, for a given initial state x 0 , an input parametrization θ and a model parameters vector p ∈ P, it is equal to one if the constraints [START_REF] Onofrio | On the Dynamics of Tumor-Immune System Interactions and Combined Chemo-and Immunotherapy[END_REF] are violated, zero otherwise. Problem ( 14) then aims at selecting the optimal control strategy such that no constraints violation occurs.

As previously explained, the randomized method consists in replacing the original problem in ( 14) by the following chance-constrained problem tolerating some violations:

min θ∈Θ J(θ) s.t. Pr X 0 (Γ)×P {g c (θ, x 0 , p, Ω) = 1} ≤ η, (15)
where the constraint is on the probability of violation, with respect to the distribution of x 0 on Γ, that we denote X 0 (Γ), and p ∼ P. This problem gives therefore a chanceconstrained formulation in the sense that we can accept a vector θ which minimizes the cost J, even if the constraints are violated for some realizations of (x 0 , p), provided that the probability of these violations is lower than η, hence small enough.

Since problem ( 15) is hard to solve, it can be simplified into the following problem, employing the empirical mean instead of the probability of the constraints violation, Given Γ ⊆ R 5 :

min θ∈Θ J(θ) s.t. N i=1 g c θ, x (i) 0 , p (i) , Ω ≤ m, (x 0 , p) (i) ∼ (X 0 (Γ) × P) , ∀i = 1, . . . , N, (16)
where m is the maximum number of allowed constraints violation. If N satisfies the condition in [START_REF] Doban | A switching control law approach for cancer immunotherapy of an evolutionary tumor growth model[END_REF], then the solution of ( 16) satisfies the constraint in problem (15) with a probability higher than 1δ.

As previously explained, problem ( 16) can be solved by simple enumeration which means that, for each possible treatment protocol defined by θ, we simulate the model [START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF] for N samples of the uncertain parameter vector p, in order to compute N i=1 g c θ, x (i) 0 , p (i) , Ω and thereby to select the treatment protocols satisfying the latter constraint. The solution is the treatment protocol defined by θ corresponding to the minimum cost J(θ). In the next section, we will explain how the iterative resolution of problems of the type ( 16) allows one to generate a sequence of sets {Ω k } N C k=1 such that the constraints violation on passing from Ω k+1 to k j=0 Ω j is smaller then η with a certain desired confidence probability 1δ.

Algorithm for RoA estimation

Given a target set Ω ⊆ R 5 + , our objective is to certify that the set Γ is such that there exists a control parametrization θ, for which at least 100 • (1η)% of the trajectories of [START_REF] Onofrio | A bi-parametric model for the tumour angiogenesis and antiangiogenesis therapy[END_REF], generated by the distributions of the initial states x 0 ∈ Γ and the uncertain parameters p, converge to Ω at time T , while satisfying constraints [START_REF] Onofrio | On the Dynamics of Tumor-Immune System Interactions and Combined Chemo-and Immunotherapy[END_REF], with a confidence higher than 1δ. Any solution of ( 16) defines a local control strategy that satisfies the constraints while minimizing the cost J(θ).

Γ generator

We suppose that we have a generator of sets Γ ⊆ R 5 + with a parametrized geometry providing a family of nested potential sets Γ, then we can compute the biggest one that is probabilistically certified through [START_REF] Iliadis | Optimizing drug regimens in cancer chemotherapy by an efficacy-toxicity mathematical model[END_REF]. In the case under study, we consider that the sets Γ have a polytopic form. Therefore, starting from Ω 0 which is in the certified region of attraction of many benign equilibriums without therapies, an iterative procedure can be designed to generate the sequence {Ω k } N C k=0 such that the trajectories starting in Ω k+1 end in k j=0

Ω j with the desired probability and without violating the constraints. In particular, we will consider sequences of sets such that Ω k ∩ Ω k+1 = ∅. Then, we keep doing this certification process until given Ω k-1 , the set Ω k is empty. Once the RoA probabilistic certification algorithm terminates, the candidate to be a probabilisitically certified RoA is the set

Ω C = N C i=1 Ω i . Note that, if x 0 ∈ Ω k , for k = 1, • • • , N C ,
this means that the trajectory of length

T will end in k-1 j=0
Ω j without violating the constraint with a certain probability, but no direct probabilistic guarantee is given regarding the convergence to the set Ω 0 .

It is not straightforward to derive a probabilistic bound on driving the states directly from the last set of the sequence Ω N C to Ω 0 . This is because the latter probability involves the accuracy and confidence parameters, η and δ. Another reason is that there is no guarantee that, given the initial state distribution X 0 (Ω k ), the distribution of the state at the end of the k-th therapeutic cycle is X 0 (Ω k-1 ), for which the probabilistic validation is performed. However, after deriving the sequence of certified sets, we can approximate the probability of driving the states from Ω N C to Ω 0 , with the corresponding certified control strategy, using Monte-Carlo simulations.

Algorithm 1 Sequence of probabilistically certified sets

Input: Ω 0 k ← 0 while Ω k ̸ = ∅ do Ω ← k j=0 Ω j repeat Generate non-empty Γ such that Ω ∩ Γ = ∅ until (16) is unfeasible for Γ if ∃ Γ such that (16) is feasible then Ω k ← Γ else Ω k ← ∅ end if k ← k + 1 end while N C ← k -1 Output: Ω C ← N C i=0 Ω i
Finally, by using Algorithm 1, we can obtain a sequence of certified sets, such that the output is the candidate to be a probabilistically certified RoA Ω C .

Note that since we focus on the estimation of the RoA for a specified control parametrization, the use of a cost function is not relevant although it could have been used in the case where some of the parameters defining the control are kept free.

Probabilistically certified RoA for a cancer model

As previously explained, considering N C treatment cycles, our objective consists in estimating the probabilistically certified RoA of model (1) that we denote Ω C . To this end, we certify a sequence of successive disjoint sets such that their union is the candidate to be a probabilistically certified RoA.

Moreover, the temporal control profiles that we consider correspond only to the hospitalization period (see Figure 1), meaning that the rest period is not included in the decision variable θ defined in Section 1, since we assume that this parameter can be esti-mated afterwards depending on the health conditions of the patient. Therefore, we propose a feedback control strategy that can be seen in an implicit way, such that at the end of each therapy period, we measure the states (patient health and tumor volume) and depending on the certified set Ω k where this measure lies, we can estimate the maximal possible recovery time (T c -T ) that the patient can take. At the end of the rest period, the certified therapy corresponding to this set is then applied, we keep doing this process until we reach the safe region Ω 0 .

The initial condition x 0 is assumed to be uniformly distributed in the set Γ while the following parameters vector is assumed to be uncertain:

p = (µ C , µ I , γ X , κ X , κ Y , δ Y , α Y , β Y , η Y ) (17)
and normally distributed in the following interval:

[0.9p nom , 1.1p nom ] , (18)
where p nom is the nominal value of each parameter, previously presented in Table 1.

The failure indicator function, which determines whether the constraints (3)-( 5) are satisfied or not, is defined on x(t|x 0 , p, θ) which is the state trajectory of (1) for a given control parametrization θ and a random sample of x 0 and p. We denote by x(T |x 0 , p, θ) the state trajectory evaluated at the end of the hospitalization period. Therefore, the failure indicator is defined as:

g c (θ, x 0 , p, Ω) := 0 if x 2 (t|x 0 , p, θ) ≥ c ∀t and x(T |x 0 , p, θ) ∈ Ω 1 otherwise ( 19
)
where Ω is a probabilistically certified target set which can be seen as the safe region to reach at the end of the cycle.

Using Algorithm 1, we can derive a sequence of probabilistically certified sets providing the probabilistically certified RoA. Firstly, we need to derive an initial target set Ω 0 , in order to initialize the certification algorithm.

Probabilistically certified initial target set Ω 0

Definition 2 Given p ∈ P (drawn according to the probability distribution P), x 0 following a uniform distribution on Ω start ⊆ R 5 + , and Ω end ⊆ R 5 + a certified target set in a neighborhood of benign equilibrium points of (1), generated by the realizations of p according to the probability distribution P. We define the certified safe reachability by:

Pr U (Ωstart)×P {x 2 (t|x 0 , p) ≥ c, ∀t > 0 and x(T |x 0 , p) ∈ Ω end } > 1 -η. ( 20
)
This means that the state trajectories having initial conditions in Ω start converge with a given probability to the set Ω end , after a time horizon T , in spite of all parametric uncertainties and without constraints violations. Note that Ω start can be equal to Ω end , in this case Ω start is slightly different than a probabilistically certified invariant set, since we don't require that the trajectories starting in Ω start stay in it, but rather to converge to it after some time T .

Firstly, we certify Ω eq such that:

Pr U (Ωeq)×P {x 2 (t|x 0 , p) ≥ c, ∀t > 0 and x(T |x 0 , p) ∈ Ω eq } > 1η.

Then, given p and x 0 following uniform distribution on Ω 0 , that we denote U(Ω 0 ), Ω 0 is determined such that:

Pr U (Ω 0 )×P {x 2 (t|x 0 , p) ≥ c, ∀t > 0 and x(T |x 0 , p) ∈ Ω eq } > 1 -η, (22)
for a given time T . Note that the set Ω eq is derived to be used as a target set for the determination of Ω 0 . In order to provide an estimation of Ω eq , we draw the distribution of the benign equilibriums of model (1) for many parameters vector samples (selected according to the probability distribution P). Then, we choose a geometry for Ω eq surrounding the benign equilibriums of the sample shown in Figure 2. Finally, we expand this set until ( 21) is not satisfied. After finding a proper geometry for the set Ω eq such that it satisfies [START_REF] Ledzewicz | On the role of the objective in the optimization of compartmental models for biomedical therapies[END_REF], we use Algorithm 1 in order to provide an estimation of the certified set Ω 0 . Note that in this case X 0 (Γ) corresponds to U(Ω 0 ) since we assume that x 0 is uniformly distributed on Ω 0 , and the target set for the states at time T denoted Ω in the definition of g c corresponds to Ω eq . Furthermore, since we deal with an uncontrolled system, (7) turns out to be a feasibility problem, where we need only to guarantee the probability condition in [START_REF] Merola | An insight into tumor dormancy equilibrium via the analysis of its domain of attraction[END_REF] by using the empirical mean over g c for N i.i.d. samples of (x 0 , p) mentioned in [START_REF] Iliadis | Optimizing drug regimens in cancer chemotherapy by an efficacy-toxicity mathematical model[END_REF], with θ = 0, n Θ = 1, and with the bound N given by [START_REF] Doban | A switching control law approach for cancer immunotherapy of an evolutionary tumor growth model[END_REF], for m = 1, η = 10 -2 and δ = 10 We assume that the set Ω 0 to be certified has the same geometry as the estimated nominal uncontrolled region of attraction Ωpnom 0 (derived in [START_REF] Moussa | Robust domain of attraction estimation for a tumor growth model[END_REF]) that we shrink until [START_REF] Merola | An insight into tumor dormancy equilibrium via the analysis of its domain of attraction[END_REF] is not satisfied given the confidence probability 1δ. There is clearly no guarantee that the set Ω 0 that we obtain is the biggest possible certified set, however, in this case, proving the existence of a set Ω 0 satisfying ( 22) is enough, since Ω 0 is only used as a target set for the Algorithm 1 allowing therefore to compute the sequence of certified sets.

Figure 2 shows the probabilistically certified RoA of the benign equilibriums Ω eq , the estimated uncontrolled nominal region of attraction Ωpnom 0 and the initial probabilistically certified target set Ω 0 for different T . Figure 3 shows the phase portrait of (1) with both the estimated nominal RoA Ωpnom 0 without control, and the certified initial target set Ω 0 for T = 60. We can see that the Ω 0 is smaller than Ωpnom 0 which shows the effects of parametric uncertainties consideration.

Validation of the estimation of Ω 0

In order to validate the estimation of the target set Ω 0 , we carry out 5000 Monte-Carlo simulations by randomly selecting the initial states as well as the model parameters according to their respective probability distributions. We can notice that in Figure 4 there are only 11 trajectories that converge to the malignant equilibrium, violating thereby the specified constraints. This corresponds to 99.78% of successful trajectories, validating therefore the imposed probabilistic bound on η.

Probabilistically certified region of attraction Ω C

We denote by Ω C the probabilistically certified region of attraction of system (1). We initialize Algorithm 1 with Ω 0 in order to derive the sequence of probabilistically certified sets providing the certified RoA for model [START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF].

We consider that the decision variable θ is defined by the following variables:

         σ I ∈ {0, 0.

Complexity analysis and computation time

The cardinality of Θ is n Θ = 300 giving the bound N ≥ 1863 according to [START_REF] Doban | A switching control law approach for cancer immunotherapy of an evolutionary tumor growth model[END_REF], for m = 1, η = 10 -2 and δ = 10 -3 . The number of simulations to be performed for each set certi-fication is N sim = N • n Θ = 558900. The required computational time to perform N sim simulations is less than 10 min using Matlab coder toolbox. Therefore, 1 simulation requires around 1.1 ms on a computer with an Intel(R) Core(TM) i5-10310U and a 2.21 GHz CPU. However, it is important to highlight two main points, the first one is the fact that the number of simulations in our case is due to the choice of solving the optimization problem by enumeration, since the cardinality of the input set Θ is relatively low. It is definitely possible to solve such problems using iterative algorithms (such as gradient descent) in order to reduce the computation time. The second point is that since the problem is completely scalable, it is possible to solve it using parallel computing which considerably reduces the computational time.

Figure 5 shows the 3 certified cycles for T = 5 days obtained using Algorithm 1, nominal and robust RoAs that have been estimated using the method presented in [START_REF] Moussa | Robust domain of attraction estimation for a tumor growth model[END_REF], where bang-bang control strategies were considered. We can see that, as the number of cycles increases, the certified RoA gets closer to the robust controlled one denoted Ω R . Even if the certified RoA remains smaller, it is important to recall that both the nominal and the robust RoAs in [START_REF] Moussa | Robust domain of attraction estimation for a tumor growth model[END_REF] have been estimated using continuous infusions of drugs, and do not consider the PK/PD effects. Note also that changing the type of the sets can change the size of the certified RoA. and the estimated robust RoA Ω R using continuous drugs infusions.

Validation of the estimation of Ω C

We approximated the probability of driving the states from Ω 3 to Ω 0 using 5000 Monte-Carlo simulations. We obtained that 99.4% of the trajectories of (1) having initial conditions in Ω 3 converge to Ω 0 using the probabilistic certified control strategies that we derived. Figure 6 shows the phase portrait of the 5000 Monte-Carlo trajectories. We can notice that only a small part of these trajectories violate the minimal constraint on immune cells density. The trajectories violating this constraints are presented in Figure 7. Figure 7: Monte-Carlo simulations, the trajectories that violate the minimal constraint on immune cells density, the blue polytope is the set Ω 3 where the initial states were selected.

Conclusion

In this paper, we presented a framework of probabilistic certification for regions of attraction which is based on the randomized methods, allowing to overcome the conservatism of worst-case robust approaches by proposing a tractable problem with probabilistic constraints. This framework has been used to derive a certified region of attraction for a cancer growth model. Furthermore, we provided a validation on the probability of driving the states to the certified safe target set with its corresponding control strategy.

The main advantages of this framework is that it is less conservative than the worstcase approach, since it is more tolerant to constraints violations in the presence of uncertainties. In addition to this, the methodology that we presented in this paper provides the control strategy corresponding to each certified initial states set, which allowed to validate the estimations using Monte-Carlo simulations. Furthermore, this approach is flexible in terms of computational complexity, since we can considerably reduce the computational time by solving the optimization problems iteratively instead of using enumeration or by using parallel computing since the problem is completely scalable.

The probabilistic certification of regions of attraction can be seen as a tool to tune the several parameters of the treatment protocols by properly choosing the model parameters and their distributions, the geometry of the regions of attraction to be certified and the control parametrization. All these choices impact the size of the certified region of attraction.

One of the future works in terms of RoA certification methodology is to enrich the set generator to consider different geometries, in order to approximate the biggest certifiable RoA. From a medical point a view, a future perspective would be to consider the synergy between the different injected drugs and to model their combined effects on the different compartments, in order to solve more challenging and seemingly realistic problems.

- 1 )u 1 2 Figure 1 :

 1121 Figure 1: Temporal open-loop control structure for each cycle, in black and yellow, respectively, the immunotherapy and the chemotherapy profiles.

Figure 2 :

 2 Figure 2: Probabilistically certified sets Ω 0 for different horizons T (days). The red circles are the benign equilibrium points of model (1) for 1000 parameters vector samples.

Figure 3 :

 3 Figure3: Phase portrait of (1), the colored trajectories represent the states (x 1 , x 2 ) for different initial conditions, the blue dots correspond to the 3 equilibrium points of the system, estimated nominal uncontrolled RoA Ωpnom 0 in dashed cyan and the estimated certified initial target set Ω 0 for T = 60 in blue.

Figure 4 :

 4 Figure 4: Monte-Carlo simulations to validate the certified target set Ω 0 .

Figure 5 :

 5 Figure 5: Probabilistically certified RoAs for 3 administration cycles, the estimated nominal RoA Ω pnom u

Figure 6 :

 6 Figure 6: Monte-Carlo simulations to validate the certified sequence of controls with their respective sets, the blue polytope is the set Ω 3 where the initial states were selected.

Table 1 :

 1 Definitions and nominal values of the parameters used in model[START_REF] Alamir | Robust feedback design for combined therapy of cancer[END_REF].

	Parameter	Definition	Numerical value	Unit	Reference
	µ C µ I	Tumor growth rate Tumor stimulated	1.0078 0.0029	10 7 • day -1 day -1	[23] [23]
		proliferation rate
	α Y	Rate of immune	0.0827	day -1	[23]
		cells influx
	β Y	Inverse threshold	0.0040	(-)	[23]
	γ X δ Y	Interaction rate Death rate	1 0.1873	10 7 • day -1 day -1	[11] [23]
	κ X	Chemotherapeutic killing parameter	1	10 7 •day -1	[11]
	κ Y	Immunotherapy administration parameter	1	10 7 • day -1	[11]
	x ∞ η Y	Fixed carrying capacity Chemo-induced loss	780 1	10 6 •cells mL/(µg•day)	[11] [23]
		on immune cells
	γ c	Patient response/ resistance	2.5	(-)	[5]
		to Etopside
	C 50c	Half-effect concentration of Etopside	10	µg • mL -1	[5]
	k 1	Chemotherapy clearance	1.6	day -1	[5], [16]
		rate from plasma
	k 2	Chemotherapy clearance	0.8	day -1	[5], [16]
		rate from tumor
	k 12	Link process between	0.4	day -1	[5], [16]
		plasma and the tumor
	V 1	Plasma volume	25	L	[5], [16]
	V 2	Effect site volume	15	L	[5], [16]
	γ i	Patient response/ resistance	2.5	(-)	[17]
		to Nivolumab
	C 50i	Half-effect concentration	32•10 -6	mg• mL -1	[17]
		of Nivolumab
	c i	Clearance rate of Nicolumab	11.6/24	day -1	[17]

Table 2 :

 2 Definitions of the different notations in Figure1.

	Parameter Definition
	σ I	Duration of immunotherapy infusion
	σ C	Duration of chemotherapy infusion
	d I	Concentration of immunotherapy infusion
	d C

Data Availability Statement

Data sharing is not applicable to this article as no new data were created in this study.


A Pivano

V O Dolocan

Chaotic dynamics of magnetic domain walls in nanowires

Keywords: numbers: 75.60.Ch, 75.10.Hk, 75.40.Mg

The nonlinear dynamics of a transverse domain wall (TDW) in Permalloy and Nickel nanostrips with two artificially patterned pinning centers is studied numerically up to rf frequencies. The phase diagram frequency -driving amplitude shows a rich variety of dynamical behaviors depending on the material parameters and the type and shape of pinning centers. We find that T-shaped traps (antinotches) create a classical double well Duffing potential that leads to a small chaotic region in the case of Nickel and a large one for Py. In contrast, the rectangular constrictions (notches) create an exponential potential that leads to larger chaotic regions interspersed with periodic windows for both Py and Ni. The influence of temperature manifests itself by enlarging the chaotic region and activating thermal jumps between the pinning sites while reducing the depinning field at low frequency in the notched strips.

I. INTRODUCTION

Nonlinear phenomena are studied in various fields ranging from atmospheric models to biophysics 1,2 . In magnetism, the nonlinear dynamics has developed with the study of spin-wave instabilities in magnetic spheres and the nonlinear motion of a single spin due to the nonlinearity of the Landau-Lifshitz equation [3][START_REF] Bertotti | Nonlinear Magnetization Dynamics in Nanosystems[END_REF][START_REF] Álvarez | [END_REF][6][7] . The emergence of chaotic motion in a single Bloch domain wall (DW) was also studied under applied magnetic field in magnetic bubble garnet thick films 8,9 . Recent results on non-linear DW dynamics discuss the current-driven motion of a confined DW in a constriction in relation to the DW spin-torque oscillators 10 , the dynamic resonant response of magnetostatically interacting DWs in parallel nanowires 11 and the stochastic resonance of a DW between two pinning centers in a nanowire 12 .

The magnetic domains separated by DWs in nanowires (or nanostrips) are highly studied nowadays for storage applications like the racetrack memory 13 or for logic devices 14 . The operation of the racetrack memory is based on the synchronous displacement by applied field or current of a series of DWs that are normally pinned at constrictions or other patterned traps. Several types of traps and nanowires (cylindrical or strip) were studied [15][16][17][18] along with the possible interaction between the DWs 19,20 . The traps create a pinning potential for the DW that depends on the DW type (vortex or transverse). In low width/diameter nanowires the transverse DW represents the stable state that can be pinned in the potential well created by the traps. Until now, the rf driven DW dynamics needed to operate at high velocity the possible devices was not studied in these systems even it is well known that the motion between potential wells can lead to chaos. The presence of chaos in devices is normally not desirable as affecting performance and different methods can be employed in electrical systems to render the systems response periodic. From a fundamental point a view, the study of the chaotic movement of a magnetic DW contributes to the understanding of chaos universal features and of the nonlinear magnetization dynamics.

In this article, we study numerically the nonlinear dynamics of a harmonic driven TDW between two pinning sites in a nanostrip up to rf frequencies. We investigate two pinning systems: one with two symmetric antinotches (T-traps) and one with two symmetric double notches (Fig. 1(a) and (b)). The pinning sites create a double potential well in both system types which leads to a complex nonlinear dynamics of the DW function of the rf applied field or current. We find that the antinotched strip leads to an well known Duffing potential, while for the notched strip a phenomenological potential was inferred 17 . The double well potential depends on the material parameters (Ni or Py), which results in different phase diagrams and chaotic regions. The temperature influences the dynamics of the DW in the nanostrip increasing globally the chaotic window and reducing the depinning field at low frequency.

This article is organized as follows. In Sec. II, we present the micromagnetic and the stochastic 1D model used to calculate the rf wall dynamics. In Sec. III, we compute and investigate the phase diagram of the DW dynamics in the different systems presented at T=0K and room temperature. Discussion and concluding remarks are presented in Sec. IV.

II. MODEL

We study numerically the dynamics of domain walls in nanostrips of Py or Ni with two types of pinning sites: notches and antinotches. The dimensions of notches and antinotches were chosen to have similar potential barriers between the pinning sites for the same material. The strip has a cross section of L y ×L z = 50×5nm 2 (with antinotches) or 60×5nm 2 (with notches). The antinotches have dimensions of l x × l y × l z = 70 × 150 × 5 nm 3 separated by 350nm. strip using the parameters of Nickel: saturation magnetization M s =477kA/m (µ 0 M s =0.6T), exchange stiffness A = 1.05 × 10 -11 J/m, damping parameter α=0.02. The dimensions of the notches are 20 × 10 × 5 nm 3 and are separated by 40nm for Ni nanostrip while for Py dimensions of 15 × 8 × 5 nm 3 were used with 30nm separation distance. Fig. 1(a) displays the equilibrium position of a head-to-head TDW in the notched strip using the parameters of Py: saturation magnetization M s =860kA/m (µ 0 M s =1.08T), exchange stiffness A = 1.3 × 10 -11 J/m, damping parameter α=0.01. An AC magnetic field or current were applied along the x-axis up to a frequency of 3GHz.

The DW dynamics was computed using 3D micromagnetic simulations with the Nmag package 21 and with the one-dimensional DW model 22,[START_REF] Thiaville | Spin Dynamics in Confined Magnetic Structures III[END_REF] . For the micromagnetic computations, the strips were discretized into a mesh with a cell size of 3nm, inferior to the exchange length (∼5nm). The average position of the DW center (x) is extracted for each applied field (in the axial x direction) along with the azimuthal angle (ψ) of magnetization in the yz plane. No magnetocrystalline anisotropy is considered and the temperature is set to T = 0K. The effect of the temperature is introduced using the stochastic 1D model. The 1D model of the DW supposes that the DW is rigid and gives a quasi-quantitative understanding of the motion of TDWs. The Langevin equations of motion of the DW are [START_REF] Boulle | [END_REF]25 :

(1 + α 2 ) Ẋ = - αγ∆ 2µ 0 M s S ∂E ∂X + γ∆ 2 H k sin 2ψ + γ 2µ 0 M s S ∂E ∂ψ + η ψ -αη X (1) (1 + α 2 ) ψ = - γ 2µ 0 M s S ∂E ∂X - γα 2 H k sin 2ψ - αγ 2∆µ 0 M s S ∂E ∂ψ + η X + αη ψ (2)
where X and ψ are the position and azimuthal angle (in the yz plane) of the DW, ∆ the DW width, S the section of the wire, γ the gyromagnetic ratio, M s the saturation magnetization, H k the DW demagnetizing field, α the damping parameter, η X and η ψ represent stochastic Gaussian noise with zero mean value and correlations

η i (t)η j (t ) = (2αk B T )/(µ 0 M s ∆Sγ 0 )δ ij δ(t -t ).
E is the potential energy of the DW that includes the internal energy, the Zeeman energy, the effects of current, the interaction energy with other DWs and the pinning energy. The pinning field associated with the pinning potential is given by H pin (x) = -1 2µ0MsS ∂Epin ∂X and the DW width variation is given by ∆(t) = ∆[Ψ(t)] = π 2A µ0M S 2 sin 2 Ψ+µ0M S H k . The potential pinning energy E pin and the equilibrium position in each system are determined from quasistatic micromagnetic simulations and are shown in Fig. 1(c) and (d) for the notched and the antinotched nanostrip. H k and ∆ are also estimated from the micromagnetic computations. The pinning potentials determined by fitting the micromagnetic results are of bistable Duffing-type for antinotched strips E pin = ax 2 + bx 4 (for Py an extend Duffing potential with a x 8 term was added). For notched strips an exponential phenomenological potential was used as before 17 :

E pin (x) = V 0 -V 1 exp - (x + x L ) 2 L 2 -exp - (x + x R ) 2 L 2
(3) with x L and x R the centers of each pinning site, L the length and V 1 the effective depth of pinning sites. For both types of pinning centers (notches and antinotches), the energy barrier between the pinning wells is controlled by the dimensions and distance between the centers.

To determine the apparition of chaos in the driven-DW dynamics as the control parameters are varied, we computed point-by-point phase diagrams for all systems with the 1D DW model for a large frequency range up to 3GHz. A similar micromagnetic computation will require an enormous execution time. The control parameters are the amplitude and frequency of the ac applied field (or current). The range of field amplitude (≤50 Oe) is chosen to have only viscous motion (low DW velocity, no precession) and/or inferior to the depinning field (for notched stripes). The frequency range (≤3GHz) is chosen to be on the same order of magnitude with access or reading/writing time in possible magnetic memories based on DW (∼nanosecond). To compare the results of 1D model with the results of micromagnetic simulations, we computed the bifurcation diagrams and the Poincaré sections showing strange attractors by both methods for all systems at fixed frequencies.

III. RESULTS

We start presenting our results with the dynamics of a DW under an applied magnetic field at zero temperature. Afterwards, we will investigate the dynamic behavior under applied current and the influence of the temperature.

A. Frequency -magnetic field phase diagram at T=0K

To identify and characterize chaos quantitatively we calculated the largest Lyapunov exponent which underlines the dynamics of chaotic behavior [START_REF] Hilborn | Chaos and Nonlinear Dynamics[END_REF]6 . Another type of quantifier which punctuates the geometric nature of trajectories in phase space, the fractal dimension (more precisely the correlation dimension), was also calculated in some special cases. The method of Lyapunov exponent is based on the exponential instability of nearby chaotic trajectories in phase space to variations in the initial conditions. Noting the difference in the nearby trajectories δx 0 (for each component of the phase space) at the initial time t 0 and δx n the difference at a later time t n the Lyapunov exponent is defined by (in the long-term limit):

λ = 1 n ln δx n δx 0 (4)
A positive average Lyapunov coefficient implies a chaotic behavior of the system. To characterize the strange attractors we also use the correlation dimension which gives the probability of finding two points in the same cell (of given radius). A noninteger value of the dimension implies a strange attractor.

We start our study with a DW pinned in a nanostrip with two symmetric antinotches and submitted to a harmonic applied magnetic field. The asymptotic movement of the DW is equivalent to a classical Duffing oscillator. The phase diagrams of largest Lyapunov exponent in the parameter space frequency -field amplitude are shown in Fig. 2(a) and (e) for Py and Ni respectively at T=0K.

The diagrams represent 1000×1000 point-by-point integration with a fourth order Runge-Kutta scheme. A clear difference is observed between the two diagrams. In the Ni strip, we determine only a small window of chaos in the parameter space. The chaos appears at low fields (≤ 12Oe) and frequencies inferior to 680MHz. Above these parameter values, the asymptotic motion of the DW remains periodic (period-1) between the two potential wells (below 1GHz) or in one potential well (above 1GHz). The periodicity of the asymptotic DW motion was confirmed micromagnetically for a number of cases at frequencies above 680MHz. A comparison between the micromagnetic and 1D model computed bifurcation diagram is shown in panel (f) for a frequency of 500MHz. The micromagnetic results display the same behavior as the 1D model. Chaotic windows appear between 4.5Oe and 10.2Oe with a large subharmonic period-3 motion in between. The transition to chaos starts with a period doubling bifurcation at 4.1Oe. The Poincaré section showing a strange (chaotic) attractor is shown in panel (g) for H=6Oe corresponding to the first chaotic window from (f) panel. The calculated correlation dimension (on 5000 points) of this attractor is 1.39.

In contrast, for Py antinotched strip, a large part of the parameter space is dominated by chaos. A periodic motion is determined for all frequencies at an applied field amplitude below 14Oe, and up to 50Oe for frequencies between 600MHz and 750MHz and above 2.5GHz. Up until 1.5GHz, the chaos starts at lower fields with increasing frequency. Two bifurcation diagrams are shown in panels (b) and (c) computed at 500MHz and 1GHz respectively. At 500MHz, the DW motion is periodic in one potential well until 22.5Oe, afterwards the motion stays periodic but takes place between the two potential wells. The onset of chaos starts with a period-doubling bifurcation at 35.7Oe and continues with a cascade of Feigenbaum-like period-doubling until 40.5Oe where it expands into a small window of chaos. At 41Oe, the 1D model predicts that the motion of the DW will revert to period-two motion with another cascade of perioddoubling onset of chaos. Micromagnetically, the fieldstep used of 1Oe does not allow us to precisely identify the periodic window between 41 and 43 Oe. In general, the 1D model gives quantitatively the same results as the micromagnetic simulation (computed with less points) at low fields and quasi-quantitatively at high fields or high frequency. At 1GHz, the DW motion is dominated by chaos with many narrow windows where chaos reverts to periodic motion (period-3 or period-5). The Poincaré section for a strange attractor is shown in panel (d) for an amplitude of the applied field of 50Oe and a frequency of 1GHz. The correlation dimension of this attractor is 1.70.

Typical properties of anharmonic oscillations like jump and hysteresis phenomena were also determined (panel (h)) for the DW motion in the antinotched nanostrip. In the frequency response spectrum above a certain critical frequency, which depends only on the nonlinear terms in the potential and the magnetic damping, bistability occurs. For Ni nanostrip, the jump frequency varies between 0.31GHz at 6Oe to 0.62GHz at 20Oe. For Py, the jump occurs at 0.94GHz for H=10Oe and at 1.5GHz for H=30Oe. A secondary resonance is observed for Py at lower frequencies and is less evident for Ni.

The motion of the DW in the notched nanostrip is detailed in Fig. 3. As two symmetric notches create a potential well that can pin a DW, the addition of another symmetric pair of notches at close distance to the first leads to a double well potential. The central barrier is highly dependent on this inter-notch distance and the potential has an exponential decay moving farther from the notch center. The external barriers are less important than for antinotches giving a depinning field around 50Oe (for antinotches H depin =300Oe). For Py (panel (a)), the Lyapunov phase diagram shows a harmonic motion of the DW for all fields up to a frequency of 1.19GHz. This is in total contrast with the antinotch case where the peri-odic motion was obtained at high frequencies. Therefore, the control of the periodicity of the movement for a large range of frequencies can be obtained by changing the type of the pinning centers. For Py notched strip, chaos appears in the DW motion at high applied field amplitude and moves to lower field with increasing frequency. A bifurcation diagram computed at 2.5GHz is shown in panel (b). At low fields, the DW motion is periodic with only one stable attractor (in one potential well). At 19Oe the asymptotic DW motion is still period-1, but now there are two stable attractors each corresponding to the two stable points (potential wells). The DW motion remains periodic in one potential well. The onset of chaos is sudden with an abrupt transition from periodic motion. One can observe several chaotic windows interspersed with periodic windows. The periodic windows are larger than in the case of antinotches with the largest being the first window (period-3) of 4.1Oe (data from 1D model). The DW asymptotic motion observed just before the period- 3 window shows intermittency (intervals of period-3 and chaotic motion). After the appearance of the period-3 window, the chaotic intervals become larger and more frequent. The third chaotic interval contains actually tiny periodic windows (mostly period-5). The Poincaré section for a strange attractor is shown in panel (c) at H=30Oe. The correlation dimension of this attractor is 1.84.

For Ni notched nanostrip (panel (d)), the DW motion is periodic at all fields at low (≤100MHz) or high frequencies (≥2.75GHz) and between 500MHz and 750MHz. The chaotic regions contain large pockets of periodic regions. A detailed bifurcation diagram (panel (e)) computed at 500MHz shows large periodic windows together with chaotic intervals. The chaotic windows contain several tiny periodic windows. The onset of chaos appears like an abrupt transition from a period-1 motion inside one potential well. The second periodic interval is also of period-1 with the motion taking place between the two potential wells. At higher applied field, the motion is period-3. The Poincaré section for a strange attractor is shown in panel (f) at H=6Oe. The correlation dimension of this attractor is 1.38.

B. Frequency -current phase diagram at T=0K

The effect of a spin-polarized rf-current (zero applied magnetic field) on the DW dynamics is shown in in Fig. 4 for the Py notched system at T=0K. The nonadiabatic spin-torque parameter β that allows the DW to be driven in the dissipative regime is taken here equal to 2α. There is still much debate over the actual value of β, with a consensus that is of the same order as the damping parameter α. The regime chosen here, β > α, corresponds to a similar behavior as the field-driven case for the velocity of the DW 25 showing a Walker-breakdown behavior. Therefore, we expect a similar phase diagram as in Fig. 3(a). In the two cases, field-driven and current-driven, we observed the same behavior, with a periodic motion until almost 1.5GHz, with chaotic motion appearing at higher frequencies. Same type of periodic windows inside the chaotic region are computed, the only difference being the slope of the boundary between the periodic and chaotic motion which depend on the value of the non-adiabatic spin-torque parameter.

The frequency-current phase diagrams for the antinotched strip (Py and Ni) are similar with the ones presented in Fig. 2 and are not shown.

C. Temperature dependence

The temperature influences the motion of the DW in two ways: on one hand the thermal noise disturbs the DW oscillation by modifying its shape and on the other hand the temperature activates jumps over the middle potential barrier at smaller applied field/current. These thermal activated jumps usually follow an Arrhenius-Néel type law [25][26][27] with the transition rate given by:

Γ = Γ 0 -∆V /k B T ( 5
)
where ∆V is the barrier height, and Γ 0 the attempt frequency at zero temperature. Γ 0 is estimated in the range 10 7 -10 12 Hz 26,27 and using an experimental time of 1ms 26 , the potential barrier a DW could surmount varies from 0.23 to 0.51eV at room temperature. This give a high estimate of the transition rate, but is more than enough to activate jumps between the two potential wells presented in section II. Using Eq. (1), we determined the stochastic motion of the DW by computing 100 realizations of the DW motion for a finite number of points (15×500) in the phase diagrams of Fig. 23. Globally, for antinotched and notched nanowires, the motion of DW at low fields (below 10-15Oe) is highly affected by thermal noise and the periodic motion can no longer be extracted. The motion is chaotic as the FFT of the DW motion shows noise for a large range of frequency. At higher fields, this thermal noise has lesser influence on the harmonic motion of the DW and periodic motion can be detected at 300K for some frequency range.

Typical temperature modification of the harmonic motion of the DW are shown in Fig. 5. The usual behavior at low fields is shown in panel (a) for the Py antinotched wire at 600MHz (H=15.7Oe). The DW motion takes place in one well (period-2) at 0K (continous line) with the FFT showing harmonics from 1 to 6 (image not shown). At 300K, the motion is still periodic in one well (period-2) with a jump probability to the second well of only 2%. We observe in this case an amplitude modulation of the DW motion (dotted line) which reflects in the FFT that shows noise in a large range of frequency (panel (b)). The noise peak around 1.2GHz (which here superposes with harmonic 2) appears in all the simulations for the Py antinotched wire at 300K. To quantify the noise in the system, we calculated the area under the large noise peak (weight intensity of the peak) for different magnetic fields and frequencies. We found that the weight of the noise peak decreases exponentially with increasing magnetic field and increasing frequency up to 1.5GHz. This corresponds to the noise being larger when the oscillations take place in one well and the middle potential barrier is higher. As the amplitude of the applied field is increased the height of the barrier diminishes and the oscillations take place between the two wells. This is the case in panel (c), where the DW motion is shown for an applied field of 22.6Oe (600MHz, Py antinotched wire). The motion is period-1 between two wells (with 100% probability) with a very small amplitude modulation at 300K. The FFT of the signal shows the odd harmonics with very small noise around them (panel (d)). Therefore, the phase diagram of the Py antinotched wire is modified as follows: for all frequecies at low applied fields (inferior of 15Oe) the DW motion becomes chaotic due to thermal noise with the motion occuring in one well. Above 15Oe and up to 1.5GHz, the phase diagram does not change much with the periodic motion showing small thermal noise and chaos appearing at smaller fields. Above 1.5GHz, the DW motion occurs in one potential well for applied fields up to 50Oe. The noise is important and the signal is chaotic. The temperature activates jumps between wells at lower fields (for example at 20Oe for 2.2GHz), with the motion still occuring in one well. At higher fields, the jumps appear more often with intermittent oscillations between the two wells.

The phase diagram for the Ni antinotched wire (Fig. 2(e)) changes more dramatically under the influence of temperature. Above 900MHz, for all applied fields the DW motion becomes chaotic (at 0K being periodic). Below 900MHz the DW motion is periodic for a large range of applied field, starting from 3Oe at 50MHz and from 7Oe at 200MHz. The periodic window diminishes with increasing frequency starting from 11Oe at 400MHz, 19Oe at 600MHz and 39.3Oe at 800MHz. An example of the temperature influence on the DW motion is shown in Fig. 5(e)-(f) for the Ni antinotched wire. For these control parameters (H=18Oe, f=400MHz), the motion is periodic between the two wells at 0K, showing a periodic attractor in the phase-space (thick red line in panel (f)).

To this periodic attractor it will correspond exactly one point in the Poincaré section in phase space. At 300K, the motion stays periodic with small noise (amplitude modulation) characterized by a periodic attractor which becomes enlarged consisting of several loops. Still the motion is periodic (or quasi-periodic) with the Poincaré section showing a cloud of points around the point corresponding to 0K. The FFT (not shown) is similar with the one in panel (d), showing harmonics 1, 3 and 5 with a small noise peak around 675MHz.

In the case of the notched wire, as the external potential barrier for depinning (escaping from the potential) is six times smaller than the antinotched wire, the depinning field is strongly affected by the temperature 19 . For Py, the thermal noise at low fields (inferior of 15Oe) is large as in the case of the antinotched wire and the motion is chaotic. Between 15 and 25Oe at low frequency (below 800MHz), the DW motion is periodic in one well with moderate noise and nonzero probability of jump to the second well. The jump probability increases with the applied field and is 29% at 15Oe and 93% at 25Oe (500MHz). Above 25Oe and up to 1.5GHz the DW depins. The depinning probabilities at 800MHz are 1% at 15Oe, 39% at 20Oe and 88% at 25Oe. Above 800MHz, the noise is too important and the motion is chotic. Therefore, the only periodic window is below 800MHz at fields between 15 and 25Oe with the window slowly reducing with increased frequency.

For the notched Ni strip, the periodic window is larger than for Py at 300K. As the middle potential barrier is lower, the DW motion takes place between the two wells even at 10Oe (at 100MHz). The DW motion is periodic between 10 and 30Oe at 100-200MHz and between 15 and 25Oe up to 1.2GHz. Above 1.2GHz, the motion is chaotic between two wells with some small periodic windows (in-termittency). A typical harmonic DW motion is shown in panel (g) of Fig. 5 for a field of 17.2Oe at 800MHz with small thermal noise at 300K (FFT in panel (h)).

IV. DISCUSSION AND CONCLUSION

The emergence of chaos in the DW motion depends largely on the damping parameter α and on the precise form of the potential and material parameters. For Ni, a value of α=0.02 was used which corresponds to the calculated value at room temperature. At low temperature, below 80K, the damping parameter increases for Ni up to a value of 0.1 28 . If an actual value of α = 0.1 is used, the phase diagram changes and for the antinotched strip only periodic motion of the DW is determined (no small chaotic region). For the notched strip, the chaotic region diminishes to a small region (inclined bubble-like) that stretches between 800MHz and 1.5GHz with the largest width of 10Oe (figure shown in Ref.29). This does not influence the behavior at room temperature, where the periodic motion is observed only up to 900MHz (antinotched wire) and up to 1.2GHz (notched wire). Similarly, for Py strip with notches, if the value of α is increased from 0.01 to a value of 0.1 the DW motion is almost completely periodic with a very small chaotic region around 2GHz and above 45Oe.

The influence of the saturation magnetization on the DW motion is more complex as the pinning energy (and demagnetization energy) and the DW width vary. For the nanostrip with notches, varying M s will change slightly the curvature of the potential and the central potential barrier height, which increases with M s 29 . When M s is decreased the DW movement will be more chaotic at high fields (more jumps between wells) and chaos will appear at slightly lower frequencies (above 40Oe). At higher M s , the chaotic motion moves to higher frequencies and larger periodic windows are determined inside the chaotic area. An example of this effect for Py notched strip is shown in Ref.29. However the dependence is not linear, and the changes in the phase diagram are small. In the case of the antinotched strip, the variation of the M s will change more drastically the form of the potential well (and equilibrium positions) and for values too far apart from the actual value (going to 1000kA/m from 860kA/m for Py) the potential changes too much and cannot be fitted anymore to a classic Duffing potential. Globally, the same behavior is expected with apparition of chaos at lower frequencies and fields.

Varying the exchange constant A will modify mostly the potential barrier height for the notched strip, while the curvature of the potential stays almost the same. The chaotic region is displaced to higher frequencies and fields when A is diminished, with changes inside the chaotic region as the pockets of periodic motion increase. Therefore, in principle, diminishing the exchange constant and increasing the saturation magnetization should displace the chaotic region to higher fields and higher frequencies.

In conclusion, we have shown that the dynamics of a TDW between two artificial pinning centers in a nanostrip can have a complex behavior. We have found that for certain parameter range, the DW motion is dominated by chaos or shows intermittency. At low temperature, there is a large region of periodic motion of the DW under applied field or current up to 3GHz and 50Oe. For Py, changing the type of pinning centers from antinotches to notches results in a high frequency to low frequency periodic DW motion, while for Ni the antinotches result in almost periodic DW motion. As typical experimental observations are realized at room temperature (for example by electrical detection of the DW motion), our stochastic results predict that periodic DW motion occurs only below 1.2-1.5GHz for all types of pinning centers (for Py and Ni). The smaller region of periodic DW motion is predicted for Py notched nanostrip, while the largest for Ni antinotched strip. In general, the Ni notched and antinotched wires seem better suited to be used at room temperature up to GHz frequencies for periodic displacement of DWs.

The chaotic dynamics is not less important in these systems with the possibility of synchronization of chaotic DW signals that could have potential applications in secure communications 30 . A future direction to be explored is chaos control in these systems which could have potential applications 31,32 . Our analysis could apply to other similar systems where chaos may appear as vortices in superconductors 33 or skyrmions in magnetic nanostructures 34 .

  Fig. 1(b) shows the equilibrium position of a head-to-head TDW in the antinotched

FIG. 1

 1 FIG. 1. (Color online) Simulated structures: planar nanowire with two symmetric double notches (a) and with two symmetric antinotches (b). The equilibrium position of a pinned DW is shown in each case. The arrows indicate the direction of magnetization. Normalized potential pinning energy for notches (c) and antinotches (d) is shown for two different materials (Ni and Py) as determined by micromagnetic simulations.

FIG. 2 .

 2 FIG. 2. (Color online) Normalized Lyapunov phase diagram for a DW at T=0K in a Py (a) and Ni (e) nanostrip with two symmetric antinotches. Dark color represents chaotic motion. The frequency f is given in units of γ0Ms, which corresponds to an absolute scale of 0-3GHz. The amplitude of the applied field h is given in units of Ms, that corresponds to 0-50Oe in an absolute scale. Bifurcation diagrams are shown in (b) and (c) computed at 500MHz and 1GHz respectively for Py. (d) Example of strange attractor obtained at H=50Oe from panel (c). (f) Bifurcation diagram for Ni at 500MHz. (g) Strange attractor for a field of H=6Oe from panel (f). The bifurcation diagrams and the strange attractor are computed with a 1D analytical model (upper panels) and by 3D micromagnetic simulations (lower panels). (h) Frequency response spectrum with jump phenomena for Py at H=10Oe and for Ni at H=20Oe.

FIG. 3

 3 FIG. 3. (Color online) Normalized Lyapunov phase diagram for a DW at T=0K in a Py (a) and Ni (d) nanostrip with two double notches. Dark color represents chaotic motion. The frequency f is given in units of γ0Ms, which corresponds to an absolute scale of 0-3GHz. The amplitude of the applied field h is given in units of Ms, that corresponds to 0-50Oe in an absolute scale for Py and to 0-30Oe for Ni. Bifurcation diagrams are shown in (b) and (e) computed at 2.5GHz for Py and 500MHz for Ni respectively. (c) Example of strange attractor obtained at H=30Oe from panel (b). (f) Strange attractor for a field of H=6Oe from panel (e). The bifurcation diagrams and the strange attractor are computed with a 1D analytical model (upper panels) and by 3D micromagnetic simulations (lower panels).

FIG. 4 .

 4 FIG. 4. (Color online). Frequency-current phase diagram for a DW at T=0K in a Py nanostrip with two double notches.

FIG. 5 .

 5 FIG. 5. (Color online). Comparaison between the DW motion at 0K (full line) and 300K (dotted line) for the Py antinotched strip (a) and (c), Ni antinotched strip (e) and Ni notched strip (g). The control parameters are: (a) 15.7Oe, 600MHz, (c) 22.6Oe, 600MHz, (e) 18Oe, 400MHz, (g) 17.2Oe, 800MHz. The power spectrum at 300K is shown in (b), (d) and (h) corresponding to (a), (c) and (g) respectively. (f) Phase-space plot of the DW motion from panel (e). The thick line corresponds to 0K.

* voicu.dolocan@im2np.fr 1 F. C. Moon, Chaotic and Fractal Dynamics: Introduction for Applied Scientists and Engineers (John Wiley & Sons, 2004).

ACKNOWLEDGMENTS

This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program "Investissements d'Avenir" supervised by the Agence Nationale pour la Recherche.

Chaotic dynamics of magnetic domain walls in nanowires -Supplemental Material

(b) M s = 900kA/m. The frequency f is given in units of γ 0 M s , which corresponds to an absolute scale of 0-2.79GHz for M s = 800kA/m and 0-3.14GHz for M s = 900kA/m. The amplitude of the applied field h is given in units of M s , that corresponds to 0-46.5Oe for M s = 800kA/m and to 0-52.3Oe for M s = 900kA/m. They are to be compared to the actual M s value (860kA/m) used in Fig. 3(a) of the main text. Same field and frequency step were used in all diagrams, so the number of points in the above diagrams is not the same. parameter α=0.1 which is closer to the expected low temperature value than α=0.02 which corresponds to the room temperature value.

II. MICROMAGNETIC SIMULATIONS -VIDEOS

Four videos are presented with their names self-explanatory:

• two videos that present a periodic (period-3) and chaotic movement of DW in a Py nanowire with notches at 2.5GHz. The type of movement is shown in the bifurcation diagram of Fig. 3(b) of the main text with the chaotic movement corresponding to the attractor shown in Fig. 3(c).

• two videos that present chaotic movement in a nanowire with antinotches: one in a Py nanowire at 1GHz and 25Oe (see Fig. All the videos are 20s long that corresponds to 40 periods of the applied magnetic field.