Elyes Cherfa
email: elyes.cherfa@univ-ubs.fr

Soraya Mesli-Kesraoui

Chouki Tibermacine
email: chouki.tibermacine@lirmm.fr

Régis Fleurquin
email: regis.fleurquin@univ-ubs.fr

Salah Sadou
email: salah.sadou@univ-ubs.fr

Identifying Metamodel Inaccurate Structures During Metamodel/Constraint Co-Evolution

Keywords: Model-driven Engineering, Metamodelling, OCL, Co-Evolution

Metamodels are subject to evolution over their lifetime. UML metamodel for instance evolved through different versions, ranging from 0.8 to 2.5 minors. These metamodels are sometimes accompanied with constraints defined using OCL (Object Constraint Language). Many works in the literature developed methods for managing and assisting the co-evolution of metamodels and their constraints. These methods enable a developer to update, in an automated (or semi-automated) way, the constraints associated to a metamodel starting from the deltas identified between versions of this metamodel. In this work we complement this assistance by notifying the developer with potential inaccurate structures in the metamodel that may be introduced during evolution. We introduce in this paper an original evolution assistance method which focuses rather on the problem (notifying metamodel inaccurate structures) than on the solution (generating OCL constraints using patterns of them). The ultimate goal of this assistance is not only to enable the developer to complete existing/updated constraints with new ones, but also to accompany her/him to further check existing constraints and to test whether they still hold. A case study is presented to show the relevance of the method.

I. INTRODUCTION: CONTEXT & PROBLEM STATEMENT

Metamodels define the syntax and part of the semantics of modeling languages. In the Model-Driven Engineering field, these are defined using languages like MOF [START_REF][END_REF]. Their specification includes sometimes constraints defined using the OMG's Object Constraint Language (OCL [START_REF] Jos | The object constraint language: getting your models ready for MDA[END_REF]). These constraints enable to make the semantics of the metamodel more precise in order to be compliant with the business domain.

For many reasons, metamodels are subject to evolution. Sometimes, this evolution is motivated by a syntax change or the introduction of new modeling elements (eg. the evolution of UML metamodel from version 1.5 to 2.0). It can also be motivated by refactoring the metamodel in order to improve and simplify its structure (eg. from UML 2.4 to UML 2.5). However, to keep the metamodel consistent, it is important to also evolve the constraints associated with it accordingly. This problem is known under the name Metamodel/Constraint coevolution. Many works in the literature proposed methods for coevolving metamodels and their constraints [START_REF] Herrmannsdoerfer | Cope: A language for the coupled evolution of metamodels and models[END_REF], [START_REF] Marković | Refactoring ocl annotated uml class diagrams[END_REF], [START_REF] Meyers | A generic in-place transformation-based approach to structured model co-evolution[END_REF], [START_REF] Wachsmuth | Metamodel adaptation and model co-adaptation[END_REF]. These works provided efficient tools for automating as much as possible the identification of impacted constraints when a change occurs in the metamodel. Sometimes they even offer a solution for rewriting the constraints impacted by the change [START_REF] Hassam | Assistance system for ocl constraints adaptation during metamodel evolution[END_REF], [START_REF] Djamel | Metamodel and constraints co-evolution: A semi automatic maintenance of ocl constraints[END_REF]. Thus, they offer a precious assistance to developers when changing their metamodels, by freeing them from managing the changes on OCL constraints.

However, these works concern only the constraints that have been already defined. It is true that a modification on the metamodel can certainly induce a modification on the existing constraints, but sometimes it can generate the need to add new constraints. Indeed, if a portion of a metamodel is enough precise without constraints, its modification can be less precise and thus requires constraints. Existing works do not propose any means to identify potential new OCL constraints when the metamodel evolves.

In this work, we tackle the problem of new constraints needed for a consistent Metamodel/Constraint coevolution. Thus, we propose a new method for assisting this coevolution which complements existing solutions. This method considers a new way of providing assistance to developers, by focusing on the problem domain (what are the structures in the metamodel that potentially cause problems and which need OCL constraints?) instead of exploring the solution domain (generating the missing constraints). Concretely we reuse two existing works to provide such assistance: one concerns the needed operators to evolve a metamodel [START_REF] Kruse | Co-Evolution of Metamodels and Model Transformations: An operator-based, stepwise approach for the impact resolution of metamodel evolution on model transformations[END_REF] and the other concerns the definition of Metamodel Inaccurate Structures (MIS) [START_REF] Cherfa | On investigating metamodel inaccurate structures[END_REF]. We believe that this way of assistance is informative, because it highlights the potential bad structures in the new version of a metamodel. This implies either writing new OCL constraints or reconsidering some existing ones.

In the rest of the paper, we present in Section II, an example illustrating the problems tackled by our work. We detail in Section III the approach we have proposed to solve this problem. Section IV presents the application of our approach on the case study: StateMachine metamodel and its evolution from version 1.5 to 2.0. The results of this evaluation are presented and discussed in Section V. Section VI discusses the different threats to the validity of the results of the of this case study. We end this paper with a conclusion (Section VIII) after presenting the related work (Section VII).

II. PROBLEM ILLUSTRATION

To illustrate in a simple way the problem addressed in our work, we take a small example from the UML StateMachine metamodel and its evolution from version 1.5 to version 2.0 (Fig. 1). This evolution was made by applying the following evolution operations: During this evolution, the reference submachineState initially typed by SubmachineState class has been Pulled up to the CompositeState class and then typed by State class. According to the adaptation operations existing in the literature, the application of this evolution operation would require a coevolution of the constraints relating to the submachineState reference. In fact, all constraints of type:

OCLExpression. submachineState [Rest], such as OCLExpression is of type StateMachine, will be coevolved after applying the first PullUp operation in:

OCLExpression. submachineState ->select(v | v. oclIsTypeOf(SubmachineState)) -> forAll(s | s[Rest])
and after applying the second PullUp operation in:

OCLExpression. submachineState -> select(v | v.oclIsTypeOf(CompositeState)) -> select(v | v.oclIsTypeOf(SubmachineState)) -> forAll(s | s[Rest])
But, the PullUp operation also generated a new MIS which may require a new OCL constraint. Indeed, all the subclasses of the State class, which now contains the reference submachine, inherit this reference. It may be possible that these sub-classes do not need to inherit the reference and thus require a new constraint to reduce the multiplicity of this inherited reference. Indeed, in the StateMachine 2.0 version, the following constraint has been added. This constraint reduces the multiplicity of the inherited reference submachine in the FinalState sub-class.

context FinalState inv inv: self.submachine->isEmpty()

In the same vein, we studied in this paper the different metamodel inaccurate structures that may require the addition of new constraints or the removal of the existing one(s), after the application of an evolution operation.

III. BACKGROUND AND PROPOSED APPROACH

Our approach relies on OCL constraint coevolution operators and Metamodel Inaccurate Structures. Thus, for a better understanding of our approach, we need to explain the main elements of these two works from the literature.

A. Metamodel coevolution operators

During the coevolution phase, a metamodel undergoes several changes. Examples of needs leading to changes on the metamodel are adding/removing concepts or refactoring it in order to simplify its comprehension and make it easy to maintain. In the literature, the coevolution problem concerning metamodels has been addressed by many studies aiming at identifying the changes that are usually applied. As a result, a set of 16 evolution operators was introduced by Wachsmuth [START_REF] Wachsmuth | Metamodel adaptation and model co-adaptation[END_REF] for coevolving metamodels and models.

Besides this, Marković and Barr [START_REF] Marković | Refactoring ocl annotated uml class diagrams[END_REF] proposed 15 rules to coevolve class diagrams and OCL constraints. These have been completed by [START_REF] Hassam | Assistance system for ocl constraints adaptation during metamodel evolution[END_REF] with 7 operations for coevolving metamodels and OCL constraints. Table I taken from [START_REF] Kruse | Co-Evolution of Metamodels and Model Transformations: An operator-based, stepwise approach for the impact resolution of metamodel evolution on model transformations[END_REF] summarises all identified EMOF metamodel evolution operators which help in solving coevolution problems. The evolution operators are grouped into four categories: 1) Those for adding or removing new modeling elements like: class, association, attribute, operation, generalization, package, etc; 2) Those that help in the manipulation of the properties such as: Move property, PullUp property, etc; 3) The third category groups the operators that manipulate hierarchies (Inheritance relations); 4) The last group represents all complex operators that were proposed to refactor structures based on patterns. OCL constraints are the main artifacts that are affected by the metamodel changes. For example, removing an attribute needs to be followed by the removal of the OCL constraints that target this attribute. Following this rule, previous research works have identified a set of OCL coevolution operators that are triggered by one or more metamodel coevolution operators. These approaches are efficient when it comes to modification and removal operators. However, when the changes concern a part of the metamodel which was not associated with constraints, these operators give no indication. This is normal since their goal is to identify the constraints that need to be changed after a change in the metamodel. But sometimes, a change in a metamodel may raise the need for new constraints to be more precise, as shown in the example of the previous section.

B. Metamodel Inaccurate Structures

In [START_REF] Cherfa | On investigating metamodel inaccurate structures[END_REF], the authors have analyzed a dataset containing 10 metamodels and more than 800 OCL constraints. The analysis was done by studying all the metamodel structures that are targeted by at least one constraint. These structures have been grouped into sets, each of which can be defined by a generic structure. The latter is named Metamodel Inaccurate Structure (MIS). Thus, the presence of a MIS in a metamodel suggests the presence of weaknesses and therefore the need for OCL constraints to make the relevant part of the metamodel more precise. One of these weaknesses is inheritance. Inheritance (Generalization) allows to generalize common properties from many classes and factorize them into one class which will be inherited from all the other classes which need one or all of these properties. The problem that may occur is that in order to generalize a property that can be found in many classes, the multiplicity and the set of possible values of this property need to be relaxed to contain all the possible multiplicities or values that the property can take in the sub-classes. As a consequence, some classes will inherit a property with a wider multiplicity or a wider range of values, which needs to be precised (restricted) by adding a new OCL constraint. The authors have pointed out 10 MIS as illustrated in Table II. 2, taken from the StateMachine metamodel, illustrates the Enumeration Type Relation MIS. This MIS occurs when a class contains an Enumeration typed attribute, and one or more incoming associations. Often, constraints are specified to precise the attribute literals value that the association must have (i.e. the type of the association is restricted by the enumeration). Here, the Pseudostate (Fig. 2) instances that are linked to the ConnectionPointReference instances through the association entry must be restricted to the literal entryPoint. This is also the case for the association exit with the literal exitPoint. A MIS suggests that a structure is possibly inaccurate. So, the metamodel designer (domain expert) needs to analyze each MIS occurrence in order to determine whether the structure is really not precise and hence needs to be completed with OCL constraints.

C. Proposed Approach

Initially, as illustrated in figure 3, we have identified 30 operators (A) that evolve metamodel from (B) version to a new version (C). To guarantee that the existing OCL constraints (E) are conforming to the new metamodel version, a set of constraints coevolution operators (D) is applied to these constraints to coevolve them to be conforming to the new metamodel version (F). Since the set of constraints (F) is incomplete and hence does not cover all the concepts of the new metamodel version, a set of MIS that can arise from the evolution operators (H) is automatically searched (G). This set of MIS instances is then analysed by metamodel designer (I) to remove all the instances that do not require OCL constraints. After this sorting, the resulted MIS instances set (J) encompasses only instances that require OCL constraints. To this end, we have developped a tool that proposes OCL (K) constraints for each MIS instance. The proposed constraints are instances of OCL patterns that are frequently used to precise MIS. After that, the domain expert will validate the set of constraints (L) relying on its domain knowledge. Consequently, the resulted set of OCL constraints (M) includes all the new OCL constraints that are applied to precise new metamodel elements resulting from the metamodel evolution. This set of constraints (M) union the set of constraints resulted from coevolving the OCL constraints (F) from the old metamodel using coevolution operators give place to a complete set of OCL constraints (N) that are conforming to the new metamodel version.

To take benefit from the knowledge acquired from the MIS and to complete the OCL coevolution, we performed an analysis over each metamodel coevolution operator. This analysis aims at identifying all the MIS that may arise following the application of a metamodel evolution operator. To do so, for each evolution operator, we have taken a set of structures and applied it, then searched for MIS instances to know which MIS may result from which metamodel evolution operator.

As a result, we obtained a set of MIS for each evolution operator. This results show which MIS can be generated after the application of an evolution operator. Table III summarizes our findings.

To be more precise, we have divided some operators such as pullUp property into pullUp association, and pullUp attribute to precise the resulted MIS. The table contains three types of metamodel evolution operators (I). Each type encompasses a set of operators (II) that we have taken from the literature. For each operator, we state if its application may cause a change that impacts OCL constraints, which require the coevolution of at least one constraint (III). For example, adding a class does not affect any OCL constraint, but removing an attribute requires removing/refactoring all the constraints that were referenced. After that, for each operator, we state if it may engender one or many MIS or not, and if it is the case, we list all the MIS that may be triggered from an evolution operator (IV). Finally, we give some conditions predominantly related to the metamodel structure, making the rise of a MIS possible (V). For example, the application of the operator "add association" may trigger the MIS "Cycle" only if the association is reflexive.

Figure 4 illustrates the application of Add Association operator in a metamodel structure. From this evolution results the Inherited Association Multiplicity Restriction MIS. This MIS occurs when a super-class contains an association with a multiplicity. Often, an OCL constraint is added to one of the sub-classes to restrict the association multiplicity. The operator Add Association does not affect the correctness of existing constraints. However, the added association may need constraints. This shows that the existing OCL coevolution operators are not sufficient to make the evolved metamodel precise. Consequently, relying on coevolution operators to coevolve existing constraints and search MIS instances raised from the metamodel evolution is an efficient combo to cover all the metamodel structures with constraints. The strength of our approach lies in its capacity to identify the need for new constraints that did not appear in the old metamodel version.

IV. CASE STUDY

To evaluate the performance of our method in OCL coevolution assistance, we conducted a case study. The purpose of this study is to provide answers to the following research questions.

A. Research Questions

• RQ1. What is the performance of our method in the coevolution of OCL constraints? This question assesses the ability of our approach to correctly coevolve the OCL constraints during the evolution of their meta-model. • RQ2. What is the performance of our method in the notification of potential new constraints related to MIS? This question aims to study the ability of our method to notify structures which will potentially require new constraints or the removal of already existing constraints and thus complete the coevolution of OCL constraints during the evolution of a metamodel. The results of this question will highlight the originality of our method compared to methods from the state of the art.

B. Selected case study

We chose to evaluate our method on the State machine metamodel evolved from version 1.5 (Fig. 5) to version 2.0 (Fig. 6). This choice is guided by: 1) accessibility of the metamodel and its OCL constraints with a complete documentation; 2) the number of significant evolution operations (73 operations that were required to evolve the metamodel from Fig. 3. Proposed approach to coevolve OCL constraints Fig. 4. Inherited Association Multiplicity Rest. resulted from Add Association operator 1.5 to 2.0), this covers many evolution operators [START_REF] Djamel | Metamodel and constraints co-evolution: A semi automatic maintenance of ocl constraints[END_REF], which allows us to validate the co-evolution of constraints such as existing approaches from the literature. Also, it enables the identification of new constraints that are needed to refine the new metamodel version; 3) the metamodel is widely used in practice which proves its good construction; 4) maneuverability: the metamodel is concise compared to the whole UML metamodel. In the following, we present the StateMachine metamodel based on its 2.0 version (Fig. 6).

State machines offer a range of concepts to model the discrete behavior of a system as a graph of vertices and transitions. The vertices are of different types: states, pseudostates, and connection point references (Fig. 6). The states describe a situation in which a condition is verified [START_REF]Object Management Group[END_REF]. States can be simple, composite, sub-machines, or final states. Unlike a simple state, a composite state can contain other states and transitions often clustered in regions. A composite state, containing at least one region, may contain several regions. In this case, the state is called orthogonal because the regions can run in parallel. A final state is the last state to go through when running a state machine. Pseudostates are transient states (Execution never stops in these states). The state machine diagram offers 10 pseudostates (entryPoint, exitPoint, initial, deepHistory, shallowHistory, join, fork, junction, terminate, or choice). On the other hand, the connection point characterizes the entry and the exit points of a sub-machine.

A transition is an arc from a source vertex to a target vertex. A transition is fired at the reception of a given event and/or if its related guard is true. On the version 2.0 [START_REF]Object Management Group[END_REF], three types of transitions were defined: internal, local, and external. An internal transition is a transition that connects a state to itself and it does not cause a state change during its execution. A local transition is a transition that does not come out of a composite state, it is always performed in the composite state. On the other hand, the execution of an external transition leads to the exit of the composite state from which it is outgoing.

C. Data

We studied the evolution of state machine diagrams from version 1.5 to 2.0. For both versions, we used the official state machine specification provided in [START_REF]Object Management Group[END_REF]. Then, for both versions, we built an ECORE metamodel and a ".ocl" file containing the OCL constraints associated with this metamodel.

The metamodel 1.5 (Table IV) is composed of 47 modeling elements: 14 classes, 1 enumeration, 5 attributes, 2 operations, and 25 associations. On this metamodel, 30 OCL constraints were specified [START_REF]Object Management Group[END_REF]. Of these 30 constraints, we were able to use only 25 constraints. The other 5 constraints have been removed either because they do not parse (they contain errors) or they relate to modeling elements contained in other packages of the UML metamodel that were not taken into account in this case study.

On the other hand, the metamodel 2.0 (Table IV) contains 69 modeling elements and 49 OCL constraints [START_REF]Object Management Group[END_REF]. Only 38 out of 49 constraints were used in this study. 11 OCL constraints were excluded for the same reasons mentioned above.

Once the OCL metamodels and constraints were collected, we identified the set of evolution operations that make Statemachine metamodel 1.5 evolving to 2.0. We manually compared the two versions of the metamodel and extracted all the evolution operations. We have identified 73 evolution

D. Data processing & method

To automatically execute these 73 operations on our tool, we have developed a Java program. This program loads the metamodel 1.5 with its 25 constraints and applies the evolution operations one by one. The application of each operation returns three outcomes: modification of the metamodel by integrating the changes generated by the operation, adapting the OCL constraints after the application of the operation, and finally, generating a set of notifications related to MIS occurrences. After the application of all 73 operations on metamodel 1.5, we obtained: 1) an evolved metamodel corresponding to version 2.0; 2) a set of adapted constraints; 3) a set of notifications for new potential constraints (corresponding to the identified MIS).

To answer the first research question, we calculated the number of OCL constraints, automatically coevolved (COEV OL) by our tool, and compared them to the constraints present in version 2.0 (coevolution made by UML designers). As a result of this comparison, we obtained 4 sets: constraints coevolved by our method and by the UML designers (COEV OL T P ); constraints coevolved by our tool and not coevolved by UML designers (COEV OL F P ), constraints not coevolved by our tool and not coevolved by UML designers (COEV OL T N ), and constraints not coevolved by our tool but coevolved by UML designers (COEV OL F N ). With these four sets, we calculated the precision and the recall of our approach, as follows:

P recision COEV OL = COEV OL T P COEV OL T P + COEV OL F P (1)
Recall COEV OL =

COEV OL T P COEV OL T P + COEV OL F N

(2)

To answer the second search question, we calculated the number of the notified MIS occurrences generated by our tool. Then, for each notification, we looked for an OCL constraint that was added or removed in version 2.0, and that avoids the MIS occurrence indicated by our tool. Indeed, notifications generated by our tool emit a suspicion that a constraint must be added or removed after the application of an evolution operation. This comparison resulted in three sets: 1) a set of MIS for which OCL constraints were added by UML designers in version 2.0 (M IS T P ); 2) a set of MIS for which constraints have not been added by the UML designers (M IS F P ); 3) a set of constraints added by UML designers that do not correspond to the identified (M IS F N ). With these three sets, we calculated the precision and the recall of our method.

P recision M IS =

M IS T P M IS T P + M IS F P

(3)

Recall M IS = M IS T P M IS T P + M IS F N (4)
V. RESULTS

A. RQ1. What are the performances of our approach in the coevolution of OCL constraints to the evolution of their meta-model?

The application of the 73 evolution operations on the StateMachine metamodel 1.5 has enabled 15 OCL constraints to be adapted from the 25 existing constraints. We obtained a precision of 1 and a recall of 1. The precision shows that the constraints coevolved by our tool have been well coevolved manually by the UML designers. This result is very interesting because the coevolution of these 15 OCL constraints took only a few seconds, which represents a considerable time saving compared to manual coevolution.

During this experiment, we noticed that the UML designers adopted other coevolutions than the ones proposed by our tool. For example, when deleting a metamodeling element, our tool proposes the removal of constraints related to this element. However, UML designers in some cases do not remove the constraints but rewrite them. Taking as an example the following constraint.

context Transition inv inv:(self.stateMachine-¿notEmpty() and not oclIsKindOf(self.stateMachine, ActivityGraph)) implies self.target.oclIsKindOf(Pseudostate) implies ((self.target.oclAsType(Pseudostate).kind = #join) implies (self.source.oclIsKindOf(State)))

Our tool proposed to remove this constraint because it refers to the stateMachine association between the Transition class and the StateMachine class, which was removed in version 2.0 of the StateMachine metamodel. However, UML designers preferred to modify it rather than to delete it as follows.

context Transition inv inv: ((self.target.oclAsType(Pseudostate).kind = #join) and (self.source.oclIsKindOf(State)))

On the other hand, automatic coevolution was able to detect constraints that should be coevolved and which have been omitted by the UML designers. Indeed, during the application of the evolution operation Remove association "top" from "StateMachine" class, our tool removed the following constraint. It took 10 years for UML designers to correct this constraint when our tool immediately detected it. These results show that the manual coevolution of OCL constraints is a long and error-prone process. The assistance of designers during the evolution of metamodels by automatically coevolving the OCL constraints is then a very useful and helpful process.

B. RQ2. What are the performances in the notification of potential new constraints related to MIS?

Regarding MIS, our tool has generated 103 notifications of potential MIS. Out of these 103 notifications (M IS T P + M IS F P ), 20 (M IS T P ) were actually about constraints that have been added in version 2.0 by UML designers. With these performances, we obtained a precision of 0.20 and a recall of 0.62 (M IS F N = 12). In the light of these results, it seems important for us to improve the performance of our approach especially concerning the number of false positives generated (reduction of the number of notifications). Indeed, we found that some notifications overlapped and were therefore counted twice. However, when looking at the results in more detail, other factors also had an impact on the results.

False negatives (constraints that do not match the MIS) are not really false negatives, because our approach looks only for MIS but not for all the constraints added in version 2.0. Indeed, the 12 constraints added by the experts and which do not correspond to MIS, are specific constraints and are highly dependent on the business domain and not MIS-related. If we consider false negatives as all constrained MIS that our tool could not notify, then we get a recall of 1. Our tool can identify all MIS related to an evolution operation.

We also found that some false positives (notification of MIS but no written constraints for this MIS) may be omitted by UML designers. Indeed, we found that 3 MIS identified by our tool have been completed in version 2.5 by 3 OCL constraints. These MIS were about the addition of the Kind enumeration in the Transition class. We strongly believe that other MIS require OCL constraints but that experts have forgotten them. For example, the operation Add association "connectionPoint" from State to "Pseudostate" generated one potential MIS that points out the need to add an OCL constraint in order to restrict the multiplicity of the connectionPoint association for the FinalState sub-class. In other words, maybe FinalState should not have connection Points. While this seems logical, no OCL constraint was found in the specification that meets this MIS. But, we are convinced that some of this MIS does require a constraint.

In the end, we found that the UML designers added 23 new OCL constraints between version 1.5 and version 2.0. Out of these 23, our tool was able to notify 11 constraints, which corresponds to 48% of new constraints added. Knowing that, we are convinced that some constraints have been omitted in the specification.

VI. THREATS TO VALIDITY

In this section, we discuss below the strategies taken to avoid some threats to the validity of the results and ensure a certain quality of our study. According to [START_REF] Runeson | Guidelines for conducting and reporting case study research in software engineering[END_REF], the threats to the validity of a case study are: construct validity, internal validity, external validity, and reliability.

A. Construct validity

Construct validity ensures that the carried measures allow to measure what it was intended to be measured. At the level of coevolution of OCL constraints, there are no potential threats of this type. Because, we have compared OCL constraints. As the OCL language is formal, there is no risk of confusion or ambiguity. However, regarding the notification of MIS, a risk may emerge because we have compared the existing constraints, written in OCL, with the notifications (textual form). To avoid this risk, we have defined a constraint template for each MIS that allows us to avoid it. We then compared these templates with the constraints written by the experts.

B. Internal validity

Internal validity concerns the consistency of the conclusions drawn on cause and effect between results. Our study is less affected by this type of validity.

C. External validity

External validity refers to the generalizability of the conclusions to other populations, domains, etc. We have studied this validity along two axes: OCL coevolution and MIS notification. On the OCL coevolution aspect, we can say that the results obtained can be generalized to other metamodels. But, on the MIS notification one, the results cannot be generalized, because they are strongly dependent on the type of the applied evolution operations. Indeed, in [START_REF] Cherfa | On investigating metamodel inaccurate structures[END_REF], studies have shown that the number of MIS in a metamodel is strongly dependent on its metamodeling elements and therefore on the type of the applied evolution operation. To avoid this threat, we chose a case study with a large and varied set of evolution operations. In fact, the case study allowed us to apply 18 types of evolution operation among the 31 previously mentioned.

D. Reliability

The reliability of the study consists in the ability of the results to be reproduced and replicated by other researchers. To increase the reliability of our study, we chose to use a metamodel known by the community and all the documentation of this case study is freely accessible via the Internet. We have also provided in this article, the set of evolutionary operations and their impact on MIS.

VII. RELATED WORK

Related works can be divided into three categories. The first one concerns the approaches that evolve metamodels. The second one encompasses the OCL constraint coevolution approaches. Lastly, the third category concerns the approaches that study the conjunctive use of metamodels, class diagrams, and the Object Constraints Language.

A. Metamodel evolution

The evolution of metamodels is a topic that was studied by many research works. These works can be classified into three categories according to [START_REF] Louis M Rose | An analysis of approaches to model migration[END_REF]. The first one concerns the approaches where the modelers encode manually the migration strategy using programming languages such as Java, or model transformation languages like QVT. In [START_REF] Garcés | Managing model adaptation by precise detection of metamodel changes[END_REF], Garces et al. compared two metamodel versions to capture the differences between them. The transformation rules are then used to adapt models automatically. [START_REF] Gruschko | Towards synchronizing models with evolving metamodels[END_REF] proposes an approach to coevolve models. It starts by detecting change either by comparing metamodels or by analyzing the change sequence applied to the metamodel old version. Then, the identified changes are classified depending on their impact on the model instances. Finally, an appropriate migration algorithm for model migration is determined. The above mentioned approaches rely on the copy rules, which can be very complex and timeconsuming to do manually. Consequently, the second category of metamodel evolution approaches appeared. This category concerns the approaches that use matching techniques to infer migration strategies by comparing the old and new metamodel versions. In these approaches [START_REF] Herrmannsdoerfer | Cope: A language for the coupled evolution of metamodels and models[END_REF], [START_REF] Mantz | Co-transformation of type and instance graphs supporting merging of types and retyping[END_REF], [START_REF] Mantz | Coevolving meta-models and their instance models: A formal approach based on graph transformation[END_REF], [START_REF] Meyers | A generic in-place transformation-based approach to structured model co-evolution[END_REF], [START_REF] Taentzer | Customizable model migration schemes for meta-model evolutions with multiplicity changes[END_REF], the coevolution rules are specified as transformation rules using a unified metamodel that represents both metamodel versions. Then, an analysis is performed to remove modeling elements that are not present in the metamodel's new version. Hence, the models can be updated to be conform to the new metamodel. The third category encompasses the approaches that capture the metamodel changes as evolution operators [START_REF] Herrmannsdoerfer | An extensive catalog of operators for the coupled evolution of metamodels and models[END_REF].

The performance of these approaches depends on the completeness of the set of evolution operators. Wachsmuth [START_REF] Wachsmuth | Metamodel adaptation and model co-adaptation[END_REF] recorded a set of operators to evolve metamodels and coevolve models. Marković and Barr [START_REF] Marković | Refactoring ocl annotated uml class diagrams[END_REF] proposed rules to coevolve class diagrams and OCL constraints. Hassam and al. [START_REF] Hassam | Assistance system for ocl constraints adaptation during metamodel evolution[END_REF] proposed operations for coevolving metamodels and OCL constraints. Kruse et al. [START_REF] Kruse | Co-Evolution of Metamodels and Model Transformations: An operator-based, stepwise approach for the impact resolution of metamodel evolution on model transformations[END_REF] gathered all the operators present in the literature and summarised them. Also, inspired by Fowler [START_REF] Fowler | Refactoring: improving the design of existing code[END_REF] who proposed code refactoring operators, Kruse et al. have proposed some complex metamodel evolution operators such as Introduce Composite Pattern. Each of the above mentioned papers have contributed in completing the metamodel evolution library by providing additional operators which improved considerably the operator-based metamodel evolution approach.

B. OCL constraints coevolution

The coevolution of OCL constraints that follows the evolution of metamodels was the center of attention of many research works. We may distinguish semi-automatic coevolution methods. For example, Hassam et al. [START_REF] Hassam | Adapting ocl constraints after a refactoring of their model using an mde process[END_REF], [START_REF] Hassam | Assistance system for ocl constraints adaptation during metamodel evolution[END_REF] proposed to highlight obsolete constraints after metamodel evolution. The semi-automatic approach shows the updates that need to be applied to the OCL constraints using the Query View Transformation (QVT) language [12]. Khelladi et al. [START_REF] Djamel | Detecting complex changes and refactorings during (meta) model evolution[END_REF], [START_REF] Djamel | Metamodel and constraints co-evolution: A semi automatic maintenance of ocl constraints[END_REF] propose a semi-automatic approach for recording metamodel changes in chronological order, and thus detecting changes and applying resolution strategies to adapt the constraints. Kusel et al. [START_REF] Kusel | Systematic co-evolution of ocl expressions[END_REF] propose semi-automatic resolution actions for the coevolution of OCL expressions in model transformations as a response to metamodel evolution. The common point between these approaches is the fact that they all rely on existing OCL constraints sets. This strategy gives positive results when the set of OCL constraints that undergoes the coevolution is complete. Also, these approaches are straightforward when the metamodel does not undergo too many evolutions that change its structure. The weakness of these approaches lies in the inability to add new constraints that complete new modeling elements.

On the other hand, other works propose fully automated methods. For example, Cabot et al. [START_REF] Cabot | Automatic integrity constraint evolution due to model subtract operations[END_REF] focus mainly on removal operations. Hence, this automatic approach has as a goal to delete constraints or parts of them. It targets the elements that are not present in the metamodel anymore. This approach focuses on removing constraints that cause inconsistencies with the conceptual schemes after subtracting operations. The approach is very powerful to avoid inconsistencies, but needs to be completed with other methods to take into consideration the other aspects of the OCL coevolution. Demuth et al. [START_REF] Demuth | Automatically generating and adapting model constraints to support coevolution of design models[END_REF], [START_REF] Demuth | Supporting the co-evolution of metamodels and constraints through incremental constraint management[END_REF] propose a template-based approach that defines a generic structure for OCL constraints that are then instantiated to update the constraints after metamodel evolution. Markovic et al. [START_REF] Marković | Refactoring ocl annotated uml class diagrams[END_REF] proposed an approach in which they formalize the most important refactoring rules for class diagrams and classify them with respect to their impact on annotated OCL constraints. This approach focuses mainly on classifying the refactoring rules depending on their impact on the OCL constraints. Batot et al. [START_REF] Batot | Heuristic-based recommendation for metamodel -ocl coevolution[END_REF] propose an automatic two-steps process to automatically coevolve metamodels and OCL constraints using a genetic algorithm. For each OCL constraint, the genetic algorithm explores the possible changes in the modeling space. The authors have defined a set of objective functions that helps to chose the changes that respect these objectives. The strength of this approach is the use of randomization using genetic operators, which can give rise to multiple OCL constraints. Just like the previous approaches, it does not fully exploit the new metamodel structure to extract all the new information from it.

While all the work that coevolve and refactor OCL constraints rely on the metamodel and the existing constraint set to perform updates, our approach relies only on metamodel structure. In addition to this, our approach can be used to add new constraints that did not appear on the old version of the metamodel.

C. Metamodel structures analysis

In our approach, we relied on the results of the analysis made in [START_REF] Cherfa | On investigating metamodel inaccurate structures[END_REF] on existing metamodels and their OCL constraints. This study had as a goal to identify metamodel structures that are often constrained with the OCL language. Other studies have been performed on metamodels and class diagrams to determine the weaknesses of both languages. Wahler et al. [START_REF] Wahler | Efficient analysis of pattern-based constraint specifications[END_REF]- [START_REF] Michael | Using patterns to develop consistent design constraints[END_REF] have captured the class diagram limits in the form of Anti-Patterns, and proposed after that a set of OCL constraints patterns that can be used to complete these imprecise structures. The main reason behind our use of MIS instead of Anti-patterns is due to the large empirical evaluation performed in [START_REF] Cherfa | On investigating metamodel inaccurate structures[END_REF]. From another perspective, Cadavid et al. have performed an analysis over metamodels and their OCL constraints in order to investigate the conjunctive use of both MOF and OCL languages. The study led to the identification of a set of OCL patterns that are repeatedly found in the metamodels. These works are powerful when it comes to writing new OCL constraints, especially during the metamodel creation phase. However, the metamodel structure analysis needs to be complemented with artifacts when it is used for refactoring metamodels.

VIII. CONCLUSION & PERSPECTIVES

In this paper, we have proposed a novel approach to assist the coevolution of metamodels and their related OCL constraints by notifying raised MIS occurrences. This approach aims at completing existing techniques that intend to coevolve metamodels with their existing OCL constraints. Our approach relies on the metamodel structure to identify the changes that may need OCL constraints.

To evaluate our approach, we have performed a case study with the UML's State Machine metamodel by comparing the 1.5 version with the 2.0 version. The results showed that out of 23 new OCL constraints, our approach was able to notify the need to define 48% of them, which is a promising result. Our approach performances are still under improvement, especially to reduce the number of false positives. Indeed, we did an exhaustive MIS search on all the metamodel, without using any other technique to reduce the false positive without losing precision.

We plan to improve our help system by using machine learning, which exploits all the data already collected, first, by studying metamodels to propose different metamodel categorizations. This may guide MIS search according to the metamodel family. We also think about using intelligent pattern recognition algorithms to search MIS instances provoked by the metamodel evolution only around the new concepts. These ideas are under study and will be presented in our future works.

Fig. 1 .

 1 Fig. 1. Evolution of the StateMachine (excerpt)

Fig. 2 .

 2 Fig. 2. Enumeration Type Relation from State Machine

Fig. 5 .Fig. 6 .

 56 Fig. 5. StateMachine metamodel version 1.5 (from [9])

TABLE V LIST

 V OF EVOLUTION OPERATIONS.

	Add class	4	Remove association	4
	Add association	23	Restrict property	2
	Add attribute	4	Generalize property	1
	Add operation	4	Extract class	1
	Add enumeration	1	Move property	2
	Add literal	6	PushDown Association 1
	Add generalization	4	PullUp attribute	1
	Remove class	4	PullUp association	2
	Remove attribute	1	Rename Element	7
	Remove generalization 1	Total	73

  context Transition inv inv: self.source.oclIsKindOf (Pseudostate) implies (self.source.oclAsType(Pseudostate).kind = #initial) implies (self.source.container = self.stateMachine.top) implies ((self.trigger-> isEmpty()) or (self.trigger.stereotype.name =' create') However, this constraint has not been removed in StateMachine version 2.0 and it is one of the 11 unparsed constraints. It was also surprising to note that this constraint dragged into all other versions of StateMachine (from 2.0, 2.1.2, 2.2, 2.3, and 2.4) without correction. It was corrected by UML designers in the latest version of UML, namely version 2.5.


C Attaccalite

M Grüning

A Marini

A real-time approach to the optical properties of solids and nano-structures: the time-dependent Bethe-Salpeter equation

Keywords: numbers: 78.20.Bh Theory, models, and numerical simulation, 78.47.je Time resolved light scattering spectroscopy, 73.22.-f Electronic structure of nanoscale materials and related systems

Many-body effects are known to play a crucial role in the electronic and optical properties of solids and nano-structures. Nevertheless the majority of theoretical and numerical approaches able to capture the influence of Coulomb correlations are restricted to the linear response regime. In this work we introduce a novel approach based on a real-time solution of the electronic dynamics. The proposed approach reduces to the well-known Bethe-Salpeter equation in the linear limit regime and it makes possible, at the same time, to investigate correlation effects in nonlinear phenomena. We show the flexibility and numerical stability of the proposed approach by calculating the dielectric constants and the effect of a strong pulse excitation in bulk h-BN.

I. INTRODUCTION

Real-time methods have proven their utility in calculating optical properties of finite systems mainly within time-dependent density functional theory (TDDFT). 1,2 On the other hand extended systems have been mostly studied by using many-body perturbation theory (MBPT) within the linear response regime. 3 The different treatment of correlation and nonlinear effects mark the range of applicability of the two approaches. The real-time TDDFT makes possible to investigate nonlinear effects like second harmonic generation 4 or hyperpolarizabilities of molecular systems. 2 However the standard approaches used to approximate the exchange-correlation functional of TDDFT treat correlation effects only on a mean-field level. As a consequence, while finite systemssuch as molecules-are well described, in the case of extended systems-such as periodic crystals and nanostructures-the real-time TDDFT does not capture the essential features of the optical absorption 5 even qualitatively.

On the contrary MBPT allows to include correlation effects using controllable and systematic approximations for the self-energy Σ, that is a one-particle operator nonlocal in space and time. Σ can be evaluated within different approximations, among which one of the most successful is the so-called GW approximation. [START_REF] Aulbur | Solid State Physics[END_REF] Since its first application to semiconductors 7 the GW self-energy has been shown to correctly reproduce quasi-particle energies and lifetimes for a wide range of materials. [START_REF] Aulbur | Solid State Physics[END_REF] Furthermore, by using the static limit of the GW self-energy as scattering potential of the Bethe-Salpeter equation (BSE), 3 it is possible to calculate response functions including electron-hole interaction effects.

In recent years, the MBPT approach has been merged with density-functional theory (DFT) by using the Kohn-Sham Hamiltonian as zeroth-order term in the perturbative expansion of the interacting Green's functions. This approach is parameter free and completely ab-initio, 5 and in this paper will be addressed as ab-initio-MBPT (Ai-MBPT) to mark the difference with the conventional MBPT. However the Ai-MBPT is a very cumbersome technique that, based on a perturbative concept, increases its level of complexity with the order of the expansion. As an example, this makes the extension of this approach beyond the linear response regime quite complex, though there have been recently some applications of the Ai-MBPT in nonlinear optics. [8][9][10] Another stringent restriction of the Ai-MBPT is that it cannot be applied when non-equilibrium phenomena take place: for example it cannot be applied to study the light emission after an ultra-fast laser pulse excitation. A generalization of MBPT to non-equilibrium situations has been proposed by Kadanoff and Baym. [START_REF] Leo | Quantum Statistical Mechanics[END_REF] In their seminal works the authors derived a set of equations for the real-time Green's functions, the Kadanoff-Baym equations (KBE's), that provide the basic tools of the nonequilibrium Green's Function theory and allow essential advances in non equilibrium statistical mechanics. [START_REF] Leo | Quantum Statistical Mechanics[END_REF] Both the standard MBPT and non-equilibrium Green's Function theory are based on the Green's function concept. This function describes the time propagation of a single particle excitation under the action of an external perturbation. In the equilibrium MBPT, due to the time translation invariance, the relevant variable used to calculate the Green's functions is the frequency ω. Instead, out of equilibrium, in all non steady-state situations, the time variables acquire a special role and much more attention is devoted to the their propagation properties. The time propagation avoids the explosive dependence, beyond the linear response, of the MBPT on high order Green's functions. Moreover the KBE are non-perturbative in the external field therefore weak and strong fields can be treated on the same footing.

One of the first attempts to apply the KBE's for investigating optical properties of semiconductors was presented in the seminal paper of Schmitt-Rink and coworkers. [START_REF] Schmitt-Rink | [END_REF] Later the KBE's were applied to study quantum wells, 13 laser excited semiconductors, 14 and luminescence 15 . However, only recently it was possible to simulate the Kadanoff-Baym dynamics in real-time. [16][17][18][19] In this work we combine a simplified version of the KBE's with DFT in such a way to obtain a parameterfree theory that is able to reproduce and predict ultrafast and nonlinear phenomena (Sec. II). This approach, that we will address as time-dependent BSE, reduces to the standard BSE for weak perturbations (Sec. II C) but, at the same time, naturally describes optical excitations beyond the linear regime. After discussing some relevant aspects of the practical implementation of our approach (Sec. III), we exemplify how it works in practice by calculating the optical absorption spectra of h-BN and the time dependent change in its electronic population due to the perturbation by means of an ultra-fast and ultrastrong laser pulse (Sec. IV).

II. THE TIME-DEPENDENT BETHE-SALPETER EQUATION

We derive here a novel approach to solve the time evolution of an electronic system with Hamiltonian coupled with an external field, Ĥ = ĥ + Ĥmb + Û ,

where U represents the electron-light interactions (see Sec. III A for its specific form). As usually done in MBPT, Ĥ is partitioned in an (effective) one-particle Hamiltonian ĥ and a part containing the many-particle effects Ĥmb .

In our derivation, we take as starting point the KBE's that we briefly introduce in Sec. II A (see e.g. Refs. 20 for a systematic treatment). Then, in Sec. II B we proceed in analogy with the equilibrium Ai-MBPT: first, we define ĥ as the Hamiltonian of the Kohn-Sham system, second we introduce the same approximations for the self-energy operator. As a result we obtain the analogous of the successful GW +BSE approach for the non-equilibrium case. Indeed in Sec. II C we show that our approach, the time-dependent BSE, reduce to the GW +BSE in the linear regime.

A. The Kadanoff-Baym equations

Within the KBE's, the time evolution of an electronic system coupled with an external field is described by the equation of motion for the non-equilibrium Green's functions [START_REF] Leo | Quantum Statistical Mechanics[END_REF][START_REF] Kremp | Quantum Statistics of Nonideal Plasmas[END_REF][START_REF] Schafer | Semiconductor Optics and Transport Phenomena: From Fundamentals to Current Topics[END_REF] , G (r, t; r ′ t ′ ). To keep the formulation as simple as possible and, being interested only in long wavelength perturbations, we expand the generic G in the eigenstates {ϕ n,k } of the ĥ Hamiltonian for a fixed momentum point k:

[G k (t 1 , t 2 )] n1n2 ≡ G n1n2,k (t 1 , t 2 ) = ϕ * n1k (r 1 )G (r 1 , t 1 ; r 2 , t 2 ) ϕ n2k (r 2 )d 3 r 1 d 3 r 2 . ( 2
)
As the external field does not break the spatial invariance of the system k is conserved. Within a second-quantization formulation of the manybody problem, the equation of motion for the Green's function described by Eq. ( 2) are obtained from those for the creation and destruction operators. However the resulting equations of motion for G k are not closed: they depend on the equations of the two-particle Green's function that in turns depends on the three-particle Green's function and so on. In order to truncate this hierarchy of equations, one introduces the self-energy operator Σ k (t 1 , t 2 ), a non-local and frequency dependent oneparticle operator that holds information of all higher order Green's functions. A further complication arises with respect to the equilibrium case because of the lack of time-traslation invariance in non-equilibrium phenomena that implies that Σ k (t 1 , t 2 ) and G k (t 1 , t 2 ) depend explicitly on both t 1 , t 2 . Then, one can define an advanced

Σ a k (G a k ), a retarded Σ r k (G r k ), a greater and a lesser Σ > k , Σ < k (G > k , G < k )
self-energy operators (Green's functions) depending on the ordering of t 1 , t 2 on the time axis. Finally, the following equation for the G < k is obtained (see e.g. Ch. 2 of Ref. 20 for more details):

i ∂ ∂t 1 G < n1n2k (t 1 , t 2 ) = δ(t 1 -t 2 )δ n1n2 + h n1n1k (t 1 )G < n1n2k (t 1 , t 2 ) + n3 U n1n3k (t 1 )G < n3n2k (t 1 , t 2 ) + n3 dt 3 Σ r n1n3k (t 1 , t 3 )G < n3n2k (t 3 , t 2 ) + Σ < n1n3k (t 1 , t 3 )G a n3n2k (t 3 , t 2 ) . (3)
This equation, together with the adjoint one for i ∂ ∂t2 G < , describes the evolution of the lesser Green's function G < k that gives access to the electron distribution (G < k (t, t)) and to the average of any one-particle operator such as for example the electron density [Eq. ( 10)], the polarization [Eq. (32)] and the current. However, in general Σ r , Σ < and the G a k depend on G > k , so that in addition to Eq. (3) the corresponding equation for the G > k has to be solved. Then, in principle, to determine the non-equilibrium Green's function in presence of an external perturbation one needs to solve the system of coupled equations for G > k , G < k , known as KBE's. Indeed, this system has been implemented within several approximations for the self-energy in model systems, 16,17 in the homogeneous electron gas, 18 and in atoms 19 . The possibility of a direct propagation in time of the KBE's provided, in these systems, valuable insights on the real-time dynamics of the electronic excitations, as their lifetime and transient effects. [16][17][18][19] Nevertheless, the enormous computational load connected to the large number of degrees of freedom de facto prevented the application of this method to crystalline solids, large molecules and nano-structures. In the next subsection we show a simplified approachgrounded on the DFT-that while capturing most of the physical effects we are interest in, makes calculation of "real-world" systems feasible.

B. The Kohn-Sham Hamiltonian and an approximation for the self-energy

In analogy to Ai-MBPT for the equilibrium case, we choose as ĥ in Eq. ( 1) the Kohn-Sham Hamiltonian, 22 ĥ = -

2 2m i ∇ 2 i + VeI + V H [ρ] + V xc [ρ], ( 4
)
where VeI is the electron-ion interaction, V H is the Hartree potential and V xc the exchange-correlation potential. Within DFT, the Kohn-Sham Hamiltonian corresponds to the independent particle system that reproduces the ground-state electronic density ρ of the full interacting system ( ĥ + Ĥmb ), that is

ρ = nk f nk |ϕ(r)| 2 , (5)
where f nk is the Kohn-Sham Fermi distribution. Equation ( 3) can be greatly simplified by choosing a static retarded approximation for the self-energy,

Σ r (t 1 , t 2 ) = Σ cohsex (t 1 ) -V xc δ(t 1 -t 2 ) (6a) Σ < (t 1 , t 2 ) = 0 (6b)
where the usual choice is Σ cohsex , the so-called Coulombhole plus screened-exchange self-energy (COHSEX). In Eq. 6a we subtracted the correlation effects already accounted by Kohn-Sham Hamiltonian ĥ. 23 The COHSEX is composed of two parts:

Σ sex (r, r ′ , t) = iW (r, r ′ ; G < )G < (r, r ′ , t), (7)
Σ coh (r, r ′ , t) = -W (r, r ′ ; G < ) 1 2 δ(r -r ′ ), (8)
where W (r, r ′ ; G < ) is the Coulomb interaction in the random-phase approximation (RPA). These two terms are obtained as a static limit of the GW self-energy (see Ch.4 of Ref. 20 and Refs. 23 and 24). With the approximation in Eqs. (6a)-(6b), Eq. ( 3) does not depend anymore on G > and it is diagonal in time:

i ∂ ∂t G < n1,n2,k (t) = h k + U k (t) + V H k [ρ] -V H k [ρ] +(Σ cohsex k (t) -V xc k [ρ]), G < k (t) n1,n2 . ( 9
)
where ρ is the density obtained from the G < as

ρ(r, t) = i n1n2k ϕ n1k (r)ϕ * n2k (r)G < n2n1k (t). (10)
Equation ( 9) is conserving [START_REF] Leo | Quantum Statistical Mechanics[END_REF] and satisfies the sum rules for the response functions because both the one-and twoparticle self-energies are obtained from the same Σ cohsex and the system of equations is solved self-consistently. 25 However, despite the full real-time COHSEX dynamics [Eq. 9] is an appealing option considerably simplifying the dynamics with respect to the KBE's, it neglects the dynamical dependence of the self-energy operator. This, in practice, induces a consistent renormalization of the quasiparticle charge 7 in addition to an opposite enhancement of the optical properties 26 . In the COHSEX approximation both effects are neglected. At the level of response properties for most of the extended systems dynamical effects are either negligible or very small (while recently it has been shown their importance for finite systems, see Refs. 27 and 28) and, for practical purposes, it has been shown that they partially cancel with the quasiparticle renormalization factors. 26 Therefore we modify Eq. 9 in order to include only the effect of the dynamical self-energy on the renormalization of the quasi-particle energies, that is the most important effect. Also in this case, our idea is to proceed in strict analogy with Ai-MBPT and to derive a real-time equation that reproduces the fruitful combination of the G 0 W 0 approximation-for the one-particle Green's function-and of the BSE with a static selfenergy-for the two-particle Green's function. Indeed the G 0 W 0 +BSE is the state-of-the-art approach to study optical properties within the Ai-MBPT. 5 To this purpose Eq. ( 8) is modified as:

i ∂ ∂t G < n1n2k (t) = h k + ∆h k + U k + V H k [ρ] -V H k [ρ] +Σ cohsex k [G < ] -Σ cohsex k [ G< ], G < k (t) n1n2 . (11)
∆h is a scissor operator 5 that applies the G 0 W 0 correction to the Kohn-Sham eigenvalues, e KS n1k ,

[∆h k ] n1,n2 = e G0W0 n1k -e KS n1k δ n1,n2 , (12)
and G< nn ′ is the solution of Eq. ( 11) for the unperturbed system (U = 0)

G< nn ′ k = i f nk δ nn ′ , (13)
where we assume that the Kohn-Sham Fermi distribution is not changed by the scissor operator. Note further that in Eq. 11,

V xc [ρ] cancels out because it is independent of G < (t).
Equation (11) is the key result of this work. It is equivalent to assume that the quasi-particle corrections modify only the single particle eigenvalues leaving unchanged the Kohn-Sham wave functions. Within AiMBPT this approximation is very successful for a wide range of materials characterized by weak correlations (see e.g. Refs. 5 and 6).

C. The linear response limit

When an external perturbation U (t) is switched on in Eq. (11), it induces a variation of the Green's function,

∆G < k (t) = G < k (t) -G< k .
In turns, this variation induces a change in the self-energy and in the Hartree potential. In the case of a strong applied laser field these changes depend on all possible orders in the external field. However for weak fields the linear term is dominant. In this regime it is possible to show analytically that Eq. ( 11) reduces to the G 0 W 0 +BSE approach 3,[START_REF] Aulbur | Solid State Physics[END_REF] . Proceeding similarly to Ref. 27 we consider the retarded density-density correlation function:

χ r (r, t; r ′ , t ′ = -i [ ρ(r, t)ρ(r ′ , t ′ ) -ρ(r, t) ρ(r ′ , t ′ ) ] θ (t -t ′ ) . ( 14
)
χ r describes the linear response of the system to a weak perturbation, represented in Eq. ( 1) by U ,

χ r (r, t; r ′ , t ′ ) = δρ(rt) δU (r ′ t ′ ) U=0 . ( 15
)
We start by expanding χ(r) in terms of the Kohn-Sham orbitals:

χ r (r, t; r ′ , t ′ ; q) = i,j,k l,m,k ′ χ r i,j,k l,m,k ′ (t, t ′ ; q) × ϕ i,k (r)ϕ * j,k+q (r)ϕ * l,k ′ (r ′ )ϕ m,k ′ +q (r ′ ), ( 16
)
where q is the momentum, and we define the matrix elements of χ r as,

χ r ij,k lm,k ′ (t, t ′ ; q) = d 3 rd 3 r ′ ϕ * i,k (r)ϕ * m,k ′ +q (r ′ )ϕ j,k+q (r)ϕ l,k ′ (r ′ ). ( 17
)
Since we are interested only in the optical response, in what follows we restrict ourselves to the case q = 0 and drop the q dependence of χ r (for the extension to finite momentum transfer see Ref. 18). Inserting the expansion for χ [Eq. ( 16)], ρ [Eq. (10)] and U (U m,nk ≡ mk|U |nk ) in Eq. ( 15) we obtain the following relation linking the matrix elements of χ r to the matrix elements of G < :

χ r ij,k lm,p (t, t ′ ) = δ iG < ji,k (t) δU lm,p (t ′ ) U=0 . (18)
Then, we can obtain the equation of motion for the matrix elements of χ r by taking the functional derivative of Eq. ( 11) with respect to U l,m,k (t),

-i ∂ ∂t χ r ij,k lm,p (t, t ′ ) = δ δU lm,p (t ′ ) [h k + ∆h k + U k (t) + V H k [ρ(t)] -V H k [ρ] + Σ k [G < (t)] -Σ k [ G< ], G < k (t)] ji . ( 19
)
Making use of the definitions in Eqs. ( 12) and ( 13), together with Eq. ( 18), it can be verified that the functional derivative of the one-electron Hamiltonian and of the external field give the contribution

δ δU l,m,p (t ′ ) [h k + ∆h k + U k , G < k (t)] ji U=0 = (e G0W0 jk -e G0W0 ik )χ r ji,k lm,p (t -t ′ ) + i(f ik -f jk )δ jl δ im δ kp . (20)
Note that, since the perturbation is weak, the Hamiltonian of the system is invariant with respect to time translation and thus χ r depends only on t -t ′ . The term in Eq. ( 19) containing the Hartree potential, that is not directly depending on the external perturbation, is expanded with respect to U l,m,k (t) by using the functional derivative chain rule and the definition in Eq. ( 18) as:

δV H ij,k [ρ (t)] = n,n ′ ,p l,m,k ′ dt ′ dt ′′ δV H ij,k [ρ (t)] δG < n ′ n,p (t ′ ) × χ r n,n ′ ,p lm,k ′ (t ′ , t ′′ )δU lm,k ′ (t ′′ ) , (21)
A similar equation can be obtained for Σ cohsex ij,k [G < (t)]. Equation (21) for Hartree potential and its analogous for the self-energy can be explicited by using

V H mn,k (t) = -2i ij G < ji,k (t) v q=0 mn,k ij,k , (22)
Σ cohsex mn,k (t) =i ij,q G < ji,(k-q) (t) W mk,i(k-q) nk,j(k-q) , (23)
where the matrix elements of v q=0 and W are labeled accordingly to Eq. ( 17). In Eq. ( 22) v q=0 is the long range part of the bare Coulomb potential, responsible for the local field effects in the BSE. Then by inserting Eq. ( 22) in Eq. ( 21) the functional derivative for the Hartree term is

δ lm,p (t ′ ) V H k [ρ(t)] -V H k [ρ], G < k (t) ji U=0 = 2i 2 (f ik -f jk ) st v q=0 ji,k st,k χ r st,k lm,p (t -t ′ ). (24)
Similarly, an analogous equation is obtained for the self-energy (see also Appendix A),

δ δU lm,p (t ′ ) Σ k [G < (t)] -Σ k [G < (t)], G < k (t) ji U=0 = -i 2 (f ik -f jk ) st,q W jk,s(k-q) ik,t(k-q) χ r st,(k-q) lm,p (t -t ′ ), ( 25
)
where we neglected the part containing the functional derivative of the screened interaction with respect to the external perturbation. This is a basic assumption of the standard BSE that is introduced in order to neglect high order vertex corrections. 3 Finally, we insert Eqs. ( 20), ( 24) and ( 25) in Eq. ( 19), and by Fourier transforming with respect to (t -t ′ ) we obtain

ω -ǫ G0W0 jk -ǫ G0W0 ik χ r ij,k lmp (ω) = i (f ik -f jk ) [δ jl δ im δ k,p + +i st,q {W jk,s(k-q) ik,t(k-q) -2v q=0 ji,k st,k }χ r st,k-q lm,p (ω) . ( 26
)
formally equivalent to the standard BSE.

III. OPTICAL PROPERTIES FROM A TIME-DEPENDENT APPROACH A. Practical solution of the time-dependent BSE

To solve Eq. ( 11) for a given electronic system [Eq. ( 1)], we start from ĥ, with its eigenvalues and eigenstates determined from a previous DFT calculations, and from the corrections ∆h k , determined e.g. from a previous G 0 W 0 calculation. Then, we switch on the external perturbation U and integrate the equations of motion using the same scheme as in Ref. 16 for the diagonal part of the G < , that is equivalent to a second order Runge-Kutta. Specifically, in Eq. ( 1) we choose to treat the interaction with the external electric field E within the direct coupling-or length gauge,

Û = -er • E(t). ( 27
)
Other choices are possible and indeed in the literature the electron-light interaction is often described within the minimal coupling-or velocity gauge (p • A), with A the vector potential. As it has been pointed out in Ref. 29, the length and velocity gauges lead to the same results only if a gauge transformation is correctly applied. However, in this respect the velocity gauge presents two main drawbacks. First, the wave functions and the boundary conditions have to be transformed by a time-dependent gauge factor T (r, t) = exp{iA(t) • r} and accordingly, in the Green's function formalism also the self-energy and the dephasing term have to be transformed. Second, within perturbation theory the velocity gauge induces divergent terms in the response function that in principle cancel each other, but that in practice lead to artificial divergences in the optical response 1,30 due to numerical precision and incomplete basis sets.

The interaction Hamiltonian U is evaluated in terms of unperturbed Kohn-Sham eigenfunctions as mk|U |nk = -E(t) mk|r|nk = -E(t) r mn,k , (28) where the dipole matrix elements r mn,k , for m = n are calculated by using the commutation relation i[H, r] = p + i[V nl , r] where V nl is the non-local part of the Hamiltonian operator. 31 Since we are interested in calculating the dielectric properties at zero momentum, we choose to work with an homogeneous electric field E(t), with no space dependence except its direction, 1 generated by a vector potential A(t) constant in space,

E(t) = - 1 c dA(t) dt . ( 29
)
Also in this case other choices consistent with the periodic boundary conditions would be possible, as for example an external potential with the cell periodicity, 32 or an electric field with a finite momentum 18 q such that q = k -k ′ . Instead, the particular form of E(t) as function of time is not specified a priori, but given as input parameter of the simulation. Indeed, the possibility of providing the form of the external field as an input is one of the key strengths of the real-time approach, potentially allowing to use the same implementation to simulate a broad range of phenomena and of experimental techniques. For example, as described in Sec. IV in order to calculate the linear optical susceptibility spectrum χ(ω) we will use a delta function E(t) = E 0 δ(t -t 0 ) (obtained from Eq.( 29) with A(t) = A 0 Θ(t -t 0 ), where t 0 is the time at which the external field is switched on). This electric field probes the system at all frequencies with the same intensity. Also, in the other example described in Sec. IV, we can use a quasi-monochromatic source E(t) = E 0 sin ω 0 t exp (-δ 2 (t -t 0 ) 2 /2) to selectively excite the system at a given frequency ω 0 . Furthermore, two or more electric fields can be used to simulate e.g. pump-probe, sum-of-frequency or wave-mixing experiments.

The macroscopic quantity that is calculated at the end of the real-time simulation is the induced polarization P(t), related to the electric displacement D(r, t) and the electric field E(r, t) by the so called material equation:

D(r, t) = ǫ 0 E(r, t) + P(r, t), (30)
that stems directly from the Maxwell equations. P(t) is obtained from G < [Eq. (11)] by,

P(t) = - 1 V n,m,k r mn,k G < nm,k (t), (31)
and from this quantity we can obtain the optical properties of the system under study. For instance, within linear response, the electric displacement D(r, t) is directly proportional, in frequency space, to the electric field as D(ω) = ǫ(ω)ǫ 0 E(ω). Therefore the polarization can be expressed as:

P(ω) = ǫ 0 (ǫ(ω) -Î)E(ω) (32)
and accordingly the optical susceptibility that describes the linear response of the system to a perturbation is χ(ω) = ǫ(ω)-Î. Then, the optical susceptibility χ(ω) can be calculated by Fourier transforming the macroscopic polarization P(ω) (or alternatively the current density j(ω)), by means of Eq. ( 32) as:

χ(ω) = P(ω) ǫ 0 E(ω) . (33)
Note that by choosing a delta-like E(t), the Fourier transform of P(t) provides directly the full spectrum of the optical susceptibility χ(ω). Beyond the linear regime, higher order response functions, χ(2) , χ(3) , . . . can be obtained (to calculate e.g. the second-or third-harmonic generation) by using a (quasi)monochromatic field as in e.g. Ref. 4; non-perturbative phenomena, such as highharmonic generation, can be analyzed instead from the power spectrum (|P (ω)| 2 ).

To summarize, the schematic flow of a time-dependent BSE simulation is shown in Fig. 1

B. Dissipative effects

In an excited electronic system dissipative effects are present due to inelastic electron scattering and (quasi)elastic scattering processes with other degrees of freedom, such as defects or phonons. Both effects contribute to the relaxation and decay of excited electronic population as well as of the decay of phase coherence, that is to a finite dephasing rate. Our approach, Eq. ( 11), does not account for dissipative effects: on the one hand the COH-SEX self-energy is real, so that the excitations lifetimes are infinite, on the other hand the electronic systems is perfectly isolated [Eq. ( 1)], so that there is no dephasing due to interaction with other degrees of freedom.

In practical calculations then we introduce a phenomenological damping to simulate dissipative effects. We implemented two different approaches. An a posteriori treatment, where at the end of the simulation (in the post-processing block of Fig. 1) the polarization (and the electric field) are multiplied by a decaying exponential function, e -t/τ , where τ is an empirical parameter. This parameter, that is compatible with the simulation length, effectively simulates the dephasing and introduces a Lorentzian broadening in the resulting absorption spectrum. This is in the same spirit of the Lorentzian broadening introduced in the linear response treatment to simulate the experimental optical spectra, and has the advantage of producing spectra with different broadening from the same real-time simulation. Nevertheless this approach is limited to the linear response case.

In order to treat dissipative effects beyond the linear regime, an imaginary term is added to the self-energy in the form of an additional term Γ k G < k (t) -G< k appearing on the r.h.s. of Eq. ( 11), with:

i [Γ k ] n1,n2 = Γ ph n1k + Γ ph n2k + Γ pop n1k δ n1,n2 , (34)
where Γ pop n1k and Γ ph n1k are respectively the lifetime of the perturbed electronic population and the dephasing rate, and are given as input parameters of the simulation.

IV. EXAMPLES

To illustrate and validate the time-dependent BSE approach and our numerical implementation, we present two examples on h-BN. This is a wide gap insulator whose optical properties are strongly renormalized by excitonic effects and for which all the parameters necessary in DFT, G 0 W 0 and response calculations, [START_REF]The DFT calculations for the h-BN have been performed using the unitary cell as in Ref. 35, a 8 × 8 × 4 k-point sampling[END_REF] are known from previous studies. [START_REF] Wirtz | [END_REF] In these examples we used Eq. ( 11), with and without including the self-energy term. We refer to the former approximation as TD-BSE, and to the latter as TD-HARTREE. Within equilibrium MBPT these two approximations would correspond to the BSE and RPA, and in fact they reduce to BSE and RPA within the linear response limit (Sec. II C).

In the first example (Fig. 2), we simulated h-BN interacting with a weak delta-like laser field. 36 As explained in Sec. III A a delta-like laser field probes all frequencies of the system and the Fourier transform of the macroscopic polarizability provides directly the susceptibility [Eq. ( 33)], and thus the dielectric constant [Eq. (32)]. Since we use a weak field, we expect negligible nonlinear effects. Then accordingly with Sec. II C, the results from TD-BSE and TD-HARTREE can be directly compared with the BSE and RPA within the standard Ai-MBPT approach. Indeed, in Figs. In the second example (Figs. 34) we exploit the potentiality of the TD-BSE approach by going beyond the linear regime and using a strong quasi-monochromatic laser field (see Sec. III A). This field excites the system selectively at one given frequency, moreover it is strong enough to induce changes in the electronic population of the system. To track these changes, during the dynamics we followed the evolution of N c (%), that is the percentage of valence electrons that are pumped by the electric field in the conduction bands (in our simulation we have 16 valence electrons in the h-BN unit cell, since core electrons are accounted using pseudopotentials). The total number of valence electrons in the system is given by the trace of G < , while N c (t) = -i ck G < cck (t) where c labels the empty states in the unperturbed system. We performed the simulations 37 for different intensities of the field (from 10 6 kW/cm to 10 9 kW/cm) and for two values of the field frequency, 5.65 eV and 8.1 eV, that depending on the level of the theory, are either at resonance or off-resonance with the system characteristic frequencies. More precisely, within TD-BSE 5.65 eV corresponds to the strong excitonic feature in the absorption spectrum, while at 8.1 eV the absorption is negligible; conversely within RPA at 5.65 eV the absorption is negligible, while 8.1 eV corresponds to the strongest feature in the spectrum (see inset of Figs. 34). The results of the various simulations are summarized in Fig. 3 that shows N c (%) as a function of the fluence-the pulse energy per unit area. For a comparison the ablation threshold of h-BN has been determined as 78 mJ/cm 2 in the femtosecond laser operational mode. 38 Finally, Figs. 4(a) and 4(b) show the evolution of N c (%) during the simulation for a field intensity of 10 9 kW/cm: one can clearly observe the enhancement in the electronic-population change due to resonance effects. The very different picture that is obtained within the two different approximations emphasizes the importance of accounting for excitonic effects (also) in the strong field regime.

V. SUMMARY

We presented a novel approach to the ab-initio calculation of optical properties in bulk materials and nano-structures that uses a time-dependent extension of the BSE. The proposed approach combines the flexibility of a real-time approach with the strength of MBPT in capturing electron-correlation. It allows to perform computationally feasable simulations beyond the linear regime of e.g. second-and third-harmonic generation, fourwave mixing, Fourier spectroscopy or pump-probe experiments. Furthermore, being the approach based on the non-equilibrium Green's Function theory, it is possible to include effects such as lifetimes, electron-electron scattering 27 and electron-phonon coupling 39 in a systematic way. Finally, we have applied the TD-BSE to the case of h-BN. First, we have calculated the optical absorption and compared it with the results from equilibrium Ai-MBPT validating our approach and numerical implementation. Then, we have shown the potentialities of the TD-BSE approach beyond the linear-regime by calculating the change in the electronic population due to the interaction with a strong quasi-monochromatic laser field.

FIG. 1 .

 1 FIG. 1. Schematic flow of a time-Dependent BSE simulation. See Sec. III A for details.

FIG. 2 .

 2 FIG. 2. h-BN:Comparison between the real-time approach and the standard RPA and BSE approaches based on the equilibrium MBPT. (a),(c): polarization P(t) generated by an electric field E(t) = Eoδ(t) within the TD-HARTREE [(a)] and TD-BSE [(c)] approximations.(b),(d): the corresponding absorption spectra (red circles) are compared with the RPA [(b)] and with the BSE [(d)] results (black line). The experimental absorption spectrum (grey shadow) is also shown as reference.

FIG. 3 .

 3 FIG. 3. h-BN:Percentage of valence electrons pumped to the conduction bands (Nc) by a quasi-monochromatic pulse as a function of the fluence. The pulse is centered either at 5.65 eV (blue boxes) or 8.1 eV (green circles) calculated within the TD-BSE (black line) and the td-HARTREE (red dashed line) approximations. In either case, each point corresponds to a separate simulation and the lines are drawn to help guide the eye. The inset shows the absorption spectra within the TD-BSE (black line) and TD-HARTREE (red dashed line) with the arrows pointing at the pump frequencies.

FIG. 4 .

 4 FIG. 4. h-BN: Percentage as function of time of valence electrons pumped to the conduction bands (Nc) by a quasi-monochromatic pulse with intensity 10 9 kW/cm 2 centered either at 5.65 eV (blue boxes) or 8.1 eV (green circles) and calculated either within the TD-BSE [(a)] or the TD-HARTREE [(b)] approximations. The brown shadow represents the fluence as function of time. The inset shows the absorption spectra within the TD-BSE (black line) and TD-HARTREE (red dashed line) with the arrows pointing at the pump frequencies.

  2(b), 2(d) the imaginary part of the dielectric constant (optical absorption) obtained by Fourier transform of the polarization in Figs. 2(a), 2(c) is indistinguishable from that obtained within equilibrium Ai-MBPT, validating our numerical implementation.

VI. ACKNOWLEDGMENTS

A. C. acknowledges useful discussions on this work with Ilya Tokatly and Lorenzo Stella.

The authors acknowledge funding by the European Community through e-I3 ETSF project (Contract Number 211956). A. M. acknowledges support from the HPC-Europa2 transnational programme (application N. 819). M. G. acknowledges support from the Fundação para a Ciência e a Tecnologia (FCT) through the Ciência 2008 programme.

This work was performed using HPC resources of the GENCI-IDRIS project No. 100063, of the CIMENT platform in Grenoble, of the CASPUR HPC in Rome, and of the Laboratory of Advanced Computation of the University of Coimbra.

Appendix A: An efficient method to update the COHSEX self-energy during the time evolution

In this appendix we show how we store and update the Σ cohsex self-energy in a efficient manner. First of all we neglect the variation of the screened interaction W (r, r ′ ; G < (t)) with respect to the G < (r, r ′ , t) by setting to zero the functional derivative ∂W/∂G (see Sec. II C). Within this approximation the Σ coh does not contribute to the time evolution, therefore only Σ sex needs to be updated:

(A1) The KBE involves the matrix elements m, k|Σ sex |m ′ , k :

where

In order to rapidly update Σ sex after a variation of G < (r, r ′ , t), we store the matrix elements:

) in such a way that Σ sex m,m ′ can be rewritten as

The M matrix can be very large, but its size can be reduced by noticing that: (i) the matrix M is Hermitian respect to the (m, m ′ ) indexes; (ii) the number of k and q points is reduced by applying the operation symmetries that are left unaltered by the applied external field; (iii) for converging optical properties only the bands close to the gap are needed (see section IV). As an additional numerical simplification we neglected all terms such that

, where M c is a cutoff that, if chosen to be M c ≃ 5 • 10 -3 does not appreciably affect the final results. In principle by using an auxiliary localized basis set 40 one can obtain a further reduction of the matrix dimensions, but in the present work we did not explore this strategy.