domain: spi info

text:
Yassir Idmessaoud
email: yassir.id-messaoud@laas.fr

Didier Dubois
email: dubois@irit.fr

Jérémie Guiochet
email: jeremie.guiochet@laas.fr

Uncertainty Elicitation and Propagation in GSN Models of Assurance Cases

Keywords: Uncertainty propagation, Belief elicitation, Goal structuring notation, Dempster-Shafer application, Safety cases

Goal structuring notation (GSN) is commonly proposed as a structuring tool for arguing about the high-level properties (e.g. safety) of a system. However, this approach does not include the representation of uncertainties that may affect arguments. Several works extend this framework using uncertainty propagation methods. The ones based on Dempster-Shafer Theory (DST) are of interest as DST can model incomplete information. However, few works relate this approach with a logical representation of relations between elements of GSN, which is actually required to justify the chosen uncertainty propagation schemes. In this paper, we improve previous proposals including a logical formalism added to GSN, and an elicitation procedure for obtaining uncertainty information from expert judgements. We briefly present an application to a case study to validate our uncertainty propagation model in GSN that takes into account both incomplete and conflicting information.

Introduction

Due to its expressiveness, the goal structuring notation (GSN) has became a defacto standard for graphical documentation of argument structures. It is notably used to argue about the safety of critical systems. However, even a well-designed GSN may include uncertainties that may question the final statement of the GSN. There is a lack of consensus about how to model these uncertainties in the argument structure. An interesting proposal [START_REF] Wang | Safety Case Confidence Propagation Based on Dempster-Shafer theory[END_REF] is to use Dempster-Shafer Theory (DST), since incomplete information can be explicitly modeled and calculated with. Several research works are investigating its use, but as presented in [START_REF] Graydon | An investigation of proposed techniques for quantifying confidence in assurance arguments[END_REF], the proposed uncertainty propagation schemes are often not clearly justified. This is mainly due to a lack of a clear definition of the logical relations between GSN elements. We investigate this issue in this paper, using DST and logical representations of arguments with new propagation models. We do not replace GSN informal notation, but build a formal model on top of it to propagate uncertainties. We also study how expert judgments can be elicited to feed our models.

The paper is structured as follows. Section 2 presents background and some related works. Sections 3 and 4 present the uncertainty propagation and elicitation methods respectively. Finally, Section 5 presents some experimental results gained by the proposed approach.

Background and related work

Goal structuring notation (GSN) is a graphical notation/language which represents argument structures (i.e., safety and assurance cases) in form of directed acyclic graphs (directed trees or arborescences). It breaks down a top claim, called "goal", into elementary sub-goals following a specific strategy and in accordance with a particular context. Each sub-goal is associated with pieces of evidence, called solutions, which support the conclusion. Figure 1 represents a typical hazard avoidance GSN pattern. To be considered as "acceptably safe" (G 1 ) all hazards (G 2 to G n ) of the system (X), listed in the context box (C 1 ), should be provably handled (Sn 1 , Sn 2 , ...) following the strategy (S 1 ). However, this symbol-based language does not specify the nature of the logical links between G 1 , G 2 , . . . G n , nor does it capture the uncertainty that may exist in the argument structure. Previous works [START_REF] Duan | Reasoning about confidence and uncertainty in assurance cases: A survey[END_REF][START_REF] Graydon | An investigation of proposed techniques for quantifying confidence in assurance arguments[END_REF] stated and discussed proposals that deal with the issue of uncertainty. An important part of these studies use probability theory to address it [START_REF] Denney | Towards measurement of confidence in safety cases[END_REF][START_REF] Guiochet | A model for safety case confidence assessment[END_REF]. For instance, some authors [START_REF] Nešić | A probabilistic model of belief in safety cases[END_REF] transform GSN into a Bayesian network (BBN) and propagate probabilities accordingly. Due to the limited expressiveness of the probabilistic framework when information is lacking, such approaches can properly deal with uncertainties due to aleatory phenomena, but they poorly represent epistemic uncertainties due to incomplete information. In addition, these methods are also very greedy in terms of data, which requires much time to collect and process.

As a generalization of probability theory, Dempster-Shafer theory [START_REF] Shafer | A mathematical theory of evidence[END_REF] (DST) offers tools to model and propagate both aleatory and epistemic uncertainty. A mass function, or basic belief assignment (BBA), is a probability distribution over the power set of the universe of possibilities (Ω), known as the frame of discernment. Formally, a mass function m : 2 Ω → [0, 1] is such that E⊆Ω m(E) = 1, and m(∅) = 0. Any subset E of Ω such that m(E) > 0 is called a focal set of m. m(E) quantifies the probability that we only know that the truth lies in E; in particular m(Ω) quantifies the amount of ignorance. A mass assignment induces a so-called belief function Bel : 2 Ω → [0, 1], defined by: Bel(A) = E⊆A m(E). It represents the sum of all the masses supporting a statement A. Belief in the negation ¬A of the statement A is represented by: Disb(A) = Bel(¬A); the value U ncer(A) = 1 -Bel(A) -Disb(A) quantifies the lack of information about A. In this paper, a conjunctive rule of combination is used for uncertainty propagation. This rule combines multiple pieces of evidence (represented by mass functions m i , with i = 1, 2, ..., n) coming from independent sources of information:

m ∩ = m 1 ⊗ m 2 such that m ∩ (A) = E1∩E2=A m 1 (E 1 ) • m 2 (E 2 ).
In DST, an additional step eliminates conflicts that may exist by means of a normalization factor (dividing m ∩ by Fig. 1. GSN example adapted from the Hazard Avoidance Pattern [START_REF] Kelly | Safety case construction and reuse using patterns[END_REF] 1 -m ∩ (∅)). This is Dempster rule of combination [START_REF] Shafer | A mathematical theory of evidence[END_REF]. This step is omitted here to indicate the presence of possibly conflicting pieces of information.

Our approach builds on some previous works (mainly [START_REF] Cyra | Support for argument structures review and assessment[END_REF][START_REF] Wang | Safety Case Confidence Propagation Based on Dempster-Shafer theory[END_REF]) that define a number of argument types and associate to each of them an uncertainty propagation formula in the setting of DST. However, in [START_REF] Cyra | Support for argument structures review and assessment[END_REF], no logical framework is provided, which prevents a formal justification of uncertainty propagation formulas. An implicit logical setting is offered in [START_REF] Wang | Safety Case Confidence Propagation Based on Dempster-Shafer theory[END_REF]. But it remains questionable since, for instance, rules that represent the relations between premises and the top-goal are modelled by equivalences. In our work we explicitly build propagation rules on a logical framework and we adopt a more flexible format using implications. A second issue is the elicitation process that collects information from experts and transforms it into belief and disbelief pairs in DST. For that, the method proposed in [START_REF] Cyra | Support for argument structures review and assessment[END_REF] and taken over in [START_REF] Wang | Safety Case Confidence Propagation Based on Dempster-Shafer theory[END_REF] is ad hoc. This transformation between expert information and (belief, disbelief) pairs is not one to one when the expert expresses no information. It yields some anomalous cases as discussed in [START_REF] Idmessaoud | Belief functions for safety arguments confidence estimation: A comparative study[END_REF]. Finally, no proposal was given in [START_REF] Cyra | Support for argument structures review and assessment[END_REF] to elicit belief on rules, while in [START_REF] Wang | Safety Case Confidence Propagation Based on Dempster-Shafer theory[END_REF] negative beliefs can be obtained, which is not acceptable. In this paper, we propose a new better-behaved elicitation approach based on the pignistic transform proposed in [START_REF] Smets | Decision making in the tbm: the necessity of the pignistic transformation[END_REF] that solves the two last issues.

From GSN to Dempster-Shafer Theory

As defined by [START_REF] Kelly | Arguing Safety -A Systematic Approach to Safety Case Management[END_REF][START_REF] Kelly | Safety case construction and reuse using patterns[END_REF], Goal Structuring Notation (GSN) is a non-formal representation that does not formally specify how premises support a conclusion. In order to model such a relation, we use logical expressions. Then we shall attach degrees of uncertainty to these logical expressions and explain how to propagate these degrees of uncertainty in the GSN, in agreement with classical logic.

(C)

The system X is acceptably safe (p) All test results are conclusive Fig. 2. A conclusion supported by one premise in GSN.

(C)

The battery is acceptably safe

(p 1 )
The risk of chemical leakage is treated

(p 2 )
The risk of explosion is treated Fig. 3. A conclusion supported by two premise in GSN.

Logical modeling of GSN

Figure 2 represents a conclusion (C) supported by a single premise (p). It describes the situation in which the conclusion (C) is true if the premise (p) supporting it is also true. This statement can be expressed using a logical implication connective: p ⇒ C standing for ¬p ∨ C, using negation ¬ and disjunction ∨. It is obvious that such an expression can only assert the validity of the conclusion (in case p holds), i.e., whether C is provably true, not whether it is provably false. Note that even if C can only be true or false, we may fail to know it. So we work in a three-state universe (belief, disbelief and ignorance). To establish disbelief in C, we need to add an implication of the form ¬p ⇒ ¬C. It describes the situation where the conclusion (C) would be false, when the premise (p) is false. We call such logical expressions "rules". Those that induce belief in C are called direct rules and those that induce disbelief are called reverse rules.

With complex systems, it is more likely to find claims supported by more than one piece of evidence. In these cases, it is necessary to consider the relationship between the premises that support the same claim. On the other hand, logical implications remain the only connective that links the evidence domain to the conclusion. Through the different GSN patterns encountered in the literature, we can identify three types:

-Conjunctive (C-Arg): It describes the case when all premises are needed to support the conclusion. The direct rule is obtained by translating this definition into a logical expression: (∧ n i p i ) ⇒ C. On the other hand, the reverse one is obtained by reversing the direct one: ¬(∧ n i p i ) ⇒ ¬C, which is equivalent to ∧ n i (¬p i ⇒ ¬C), a conjunction of simple rules. -Disjunctive (D-Arg): It describes the case when one premise is enough to support the whole conclusion. The corresponding rules are: and(∧ n i ¬p i ) ⇒ ¬C (reverse). -Hybrid (H-Arg): It describes the case where each premise supports the conclusion to some extent, but their conjunction does it to a larger extent. This rule type could be considered as a general type which includes the two previous ones. In fact, conjunctive and disjunctive types correspond to limit cases of the hybrid one. Figure 3 represents an example of the conjunctive type. To assert that the battery is acceptably safe, all risks of chemical leakage and explosion should be treated. It gives the expression: (p 1 ∧ p 2 ) ⇒ C. On the other hand, if one of the risks remains present we may assert that the battery is unsafe, which gives the expressions: ¬p 1 ⇒ ¬C and ¬p 2 ⇒ ¬C.

∧ n i (p i ⇒ C) (direct),
All rules defined above will be used to build our uncertainty propagation model. Since the conjunctive and disjunctive types represent a special case of the hybrid one, we will only present the last one. However, it is simple to deduce their expressions from the general formula.

Uncertainty propagation model

In order to build our uncertainty propagation model, we define two kinds of parameters:

-Uncertainty on premises: It is modeled as a mass function on each premise of the argument: < m 1 p , ..., m n p >. m i p assigns a mass to the premise p i , one on its negation (¬p i ) and one on the tautology (Ω, representing ignorance) summing to 1.

-Uncertainty on rules: It is used to evaluate the impact of premises on a conclusion. We associate a simple support function [START_REF] Shafer | A mathematical theory of evidence[END_REF] to each rule r of the argument type. Each simple support function consists in assigning a mass m r (r) = s to the rule and another one m r (Ω) = 1 -s to the tautology, these weights summing to 1. The set of mass functions is formally defined as: < m ⇒ , m i ⇒ , m ⇐ , m i ⇐ >, where: m ⇒ and m ⇐ represent, respectively, direct and reverse conjunctive mass functions that assign support to rules (∧ n i p i ) ⇒ C and (∧ n i ¬p i ) ⇒ ¬C, respectively. m i ⇒ , and m i ⇐ respectively, assign support to elementary rules p i ⇒ C and ¬p i ⇒ ¬C occurring in the disjunctive type.

Using the conjunctive rule of combination presented in section 2, to merge the masses on the rules (conjunctive and disjunctive ones) with the masses on premises (m

∩ = m ⇒ ⊗ m ⇐ ⊗ m i ⇒ ⊗ m i ⇐ ⊗ m i p )
, we quantify the uncertainty on the conclusion C [START_REF] Chatalic | An approach to approximate reasoning based on Dempster rule of combination[END_REF]. Since we work on a two-state frame of discernment for both premises Ω p = {p i , ¬p i } and conclusion Ω C = {C, ¬C}, masses and (dis-)belief degrees on premises, rules and the conclusion are equal. For instance,

m C (C) = Bel C (C) and m C (¬C) = Bel C (¬C) = Disb C (C).
We can prove the following results, by projecting m ∩ on the universe Ω C = {C, ¬C}:

Bel C (C) =Bel ⇒ ([∧ n i=1 p i ] ⇒ C) • n i=1 {Bel i p (p i ) • [1 -Bel i ⇒ (p i ⇒ C)]} + {1 - n i=1 [1 -Bel i p (p i ) • Bel i ⇒ (p i ⇒ C)]} (1)
Disb C (C) =Bel ⇐ ([∧ n i=1 ¬p i ] ⇒ ¬C) • n i=1 {Disb i p (p i ) • [1 -Bel i ⇐ (¬p i ⇒ ¬C)]} + {1 - n i=1 [1 -Disb i p (p i ) • Bel i ⇐ (¬p i ⇒ ¬C)]} (2)
Where:

-

Bel C (C) (resp. Disb C (C)): the degree of belief (resp. disbelief) in the con- clusion C obtained by projection of m ∩ on Ω C . -Bel i p (p i ) (resp. Disb i p (p i )): the degree of belief (resp. disbelief) in the i th premise. -Bel ⇒ ([∧ n i=1 p i ] ⇒ C) (resp. Bel i ⇐ (¬p i ⇒ ¬C))
: the degree of belief in the direct conjunctive rule (resp. i th reverse rule).

We can notice that each formula (1) and ( 2) is the result of the summation of two terms. The first part expresses a generalized conjunction (the product), and the second part reflects a generalized disjunction (the probabilistic sum 1 -(1 -a)(1 -b)). To extract propagation formulas for the pure conjunctive type (C-Arg), it is enough to set to zero the masses on the direct rules (Bel i ⇒ (p i ⇒ C)) and the mass on the conjunctive reverse rule (Bel ⇐ ([∧ n i ¬p i ] ⇒ ¬C)). Similarly, to derive the pure disjunctive formulas (D-Arg), we set to zero the mass on the conjunctive direct rule (Bel ⇒ ([∧ n i=1 p i ] ⇒ C)) and the masses on the reverse rules (Bel i ⇐ (¬p i ⇒ ¬C)). We obtain:

C-Arg :

Bel C (C) = Bel ⇒ ([∧ n i=1 p i ] ⇒ C) • n i=1 Bel p (p i ) Disb C (C) = 1 - n i=1 [1 -Disb i p (p i ) • Bel i ⇐ (¬p i ⇒ ¬C)] D-Arg : Bel C (C) = 1 - n i=1 [1 -Bel i p (p i ) • Bel i ⇒ (p i ⇒ C)] Disb C (C) = Bel ⇐ ([∧ n i=1 ¬p i ] ⇒ ¬C) • n i=1 Disb i p (p i )
Note that the belief (resp. disbelief) degree of the conclusion (Bel C (C)) only depends on the belief (resp. disbelief) degree of premises (Bel i p (p i )) and of the corresponding direct (reverse) rules (Bel ⇒ and Bel i ⇒ ). However, we observe in some cases that the sum of belief and disbelief of the conclusion, as calculated above, is greater than 1 which is not coherent. This is when the mass m ∩ (∅) > 0. It is then counted in both sums defining the degrees of belief and disbelief. It indicates the presence of conflict between premises and rules. The coherence property Bel C (C) + Disb C (C) ≤ 1 always hold if m ∩ (∅) = 0. If it is not null, the conflict mass (3) should be subtracted from both belief and disbelief values, in order to get genuine contradiction-free degrees of belief and disbelief that respect the coherence property.

In [START_REF] Idmessaoud | Quantifying confidence of safety cases with belief functions[END_REF], we provided a recursive equation to compute m n ∩ (∅) for n premises when we know m n-1 ∩ (∅):

m n ∩ (∅) = Bel n-1 C (C) • m n (¬p n ∧ ¬C) + Disb n-1 C (C) • m n (p n ∧ C) + m n-1 ∩ (∅) (3)
where:

-Bel n-1 C (C) = {1 - n-1 i=1 [1 -Bel i p (p i ) • Bel i ⇒ (p i ⇒ C)]} -m n-1 ∩ (∅) -Disb n-1 C (C) = {1 - n-1 i=1 [1 -Disb i p (p i ) • Bel i ⇐ (¬p i ⇒ ¬C)]} -m n-1 ∩ (∅) -m i (p i ∧ C) = Bel i p (p i ) • Bel i ⇒ (p i ⇒ C) -m i (¬p i ∧ ¬C) = Disb i p (p i ) • Bel i ⇐ (¬p i ⇒ ¬C
) Remark: D-Arg and C-Arg are conflict-free. Assuming that rule masses are maximal (= 1), for n = 2 we get:

Bel C (C) = Bel 1 p (p 1 ) • Bel 2 p (p 2 ) (for C-Arg) and Bel C (C) = Bel 1 p (p 1 ) + Bel 2 p (p 2 ) -Bel 1 p (p 1 ) • Bel 2 p (p 2 ) (for D-Arg).

Belief and Disbelief elicitation

The model of uncertainty propagation presented above requires two types of inputs in order to compute belief and disbelief degrees of a conclusion: Belief/Disbelief on the rules and on the premises. These two information items will be directly collected from experts. To give their assessment about a premise or a rule, experts are asked to fill in an evaluation matrix, presented in Figure 4. Each point of this matrix corresponds to a strength of decision, denoted by Dec(A), and a degree of confidence in this decision, denoted by Conf (A) attached to a proposition A. In a scale of 5 equidistant items, decision describes which side the expert leans towards: From the rejection (Dec(A) = 0) of a claim A, to its acceptance (Dec(A) = 1). It is formally the same as a degree of probability. On the other hand, confidence reflects the amount of information an expert possesses that can justify a decision. There are 6 equidistant levels of the confidence scale, from "Lack of confidence" Conf (A) = 0 to "For sure" Conf (A) = 1.

In Figure 5, we present four extreme expert assessments (see the black dot). The upper matrices represent the case of total confidence. The assessor rejects (resp. accepts) the claim in Figure 5.a (resp. 5.b). It corresponds to a maximal disbelief (resp. belief) degree. In contrast, the lower matrices represent resp. the cases of total conflict (Figure 5.c) and ignorance (Figure 5.d). In both cases, the expert cannot make a clear decision either because he has as a lot of information both to support and reject the claim (Conf (A) = 1), or because he has no information (Conf (A) = 0). In contrast to other works [START_REF] Cyra | Support for argument structures review and assessment[END_REF][START_REF] Wang | Confidence Assessment Framework for Safety Arguments[END_REF], we allow the assessor to use a midpoint value (Dec(A) = 1/2) to show full hesitancy. Uncertainty on premises: To be used in equations ( 1) and ( 2), the pair (decision, confidence) is translated into a triple (belief, disbelief and uncertainty).

To do so, we use the formula proposed in [START_REF] Cyra | Support for argument structures review and assessment[END_REF], which defines confidence as the sum of belief and disbelief degrees (equation (4), left). On the other hand, we consider decision as the pignistic transform [START_REF] Smets | Decision making in the tbm: the necessity of the pignistic transformation[END_REF] that turns a mass into a probability (equation (4), right). So, we solve the following system for Bel(p) and Disb(p):

Conf (p) = Bel(p) + Disb(p); Dec(p) = 1 + Bel(p) -Disb(p) 2 (4)
However, as indicated in [START_REF] Idmessaoud | Belief functions for safety arguments confidence estimation: A comparative study[END_REF], the pignistic transform can generate negative belief and disbelief values when the pair (Dec, Conf ) given by an expert lies outside the triangle shown in Figure 4. Known as "Josang Triangle" [START_REF] Jøsang | Subjective logic[END_REF], it represents a constraint that brackets decision Dec(p) between two values:

1 -Conf (p) 2 ≤ Dec(p) ≤ 1 + Conf (p) 2 (5)
It guarantees that all clear-cut decisions (rejection or acceptance) are made only when the confidence level is maximal. To avoid negative belief and disbelief values, we must adjust the decision value to respect constraint [START_REF] Duan | Reasoning about confidence and uncertainty in assurance cases: A survey[END_REF]. Therefore, when Dec(p) < 1-Conf (p) 2 (rejection: black dots in Fig. 4), we set

Dec(p) = 1-Conf (p) 2
. On the other hand, when Dec(p) > 1+Conf (p) 2 (acceptance: grey dots in Fig. 4), we set Dec(p) = 1+Conf (p) 2 .

Example 1. Suppose we get the following assessments on two goals (p 1 ) and (p 2 ):

p 1 : Opposable with high confidence (Dec(p 1 ) = 0.25, Conf (p 1 ) = 0.6).

p 2 : Acceptable with very high confidence (Dec(p 2 ) = 1, Conf (p 2 ) = 0.8).

To calculate Bel(p i ) and Disb(p i ), we write them in terms of Dec(p i ) and Conf (p i ), from (4):

Bel(p) = Conf (p)-1 2 +Dec(p), Disb(p) = Conf (p)+1 2 -Dec(p).
We can notice that the assessment for p 1 is inside the triangle in the matrix (Figure 4). Hence, there is no need to adjust the values: On the other hand, the assessment for p 2 is situated outside the triangle. In this case, we can be sure that the decision degree must be adjusted in accordance with the confidence value to get correct inputs. Before adjustment, we find a negative value of disbelief, which does not make sense: Bel(p 2 ) = 0.8-1 2 + 1 = 0.9 and Disb(p 2 ) = 0.8+1 2 -1 = -0.1. Following the description above, we set Dec(p 2 ) = 1+Conf (p2) 2 = 1+0.8 2 = 0.9. Then we find that Bel(p 2 ) = 0.8,

Bel(p 1 ) = 0.
Disb(p 2 ) = 0 and U ncer(p 2 ) = 1 -Bel(p 2 ) -Disb(p 2 ) = 0.2.
Uncertainty on rules: Assuming clear-cut knowledge about some (or all) premises (Bel i p (p i ), Disb i p (p i ) ∈ {0, 1}) and total ignorance about the others (U ncer i p (p i ) = 1), we notice that Bel C (C) and Disb C (C) in ( 1) and ( 2) are equal to rule masses. For instance, in the case of a conclusion C supported by two premises p 1 and p 2 , assuming total acceptance of these two premises with maximal confidence, we get:

Bel C (C) = Bel ⇒ ([p 1 ∧ p 2 ] ⇒ C).
While assuming total rejection with maximal confidence of p 1 , and total ignorance about p 2 , we get:

Disb C (C) = Bel ⇐ (¬p 1 ⇒ ¬C).
In order to collect masses on rules, under the assumption mentioned above (sure truth, sure falsity or ignorance on premises) we use the same approach as for eliciting uncertainty on premises. First, using the evaluation matrix (Figure 4), we take the expert opinions about the conclusion (which corresponds to the rules masses under those assumptions). Then, we change them to belief values using transformation formulas [START_REF] Denney | Towards measurement of confidence in safety cases[END_REF].

Moreover, we assume that a rule is either accepted or discarded, but not negated. In fact, for any rule R : p ⇒ C we do not consider a positive disbelief because this would imply a belief in ¬(p ⇒ C) = p ∧ ¬C, i.e., ¬R which is not a rule. So we only assign mass to a rule or to the tautology; the latter is the extent to which a rule is discarded. This constraint impacts the allowed pairs (Dec, Conf) for the expert. The latter is constrained to choose only a decision on the positive side (from "No decision" to "acceptable") for direct rules. On the contrary, (s)he can only choose a negative decision (from "rejectable" to "No decision") for the reverse rules. Formulas in (4) are used to derive the degrees of belief on rules. We can notice in this example that both cases respect the Josang constraint (5). Hence, there is no need to adjust the decision value. Using (4) for the direct rule R 1 :

Bel ⇒ (R 1 ) = Bel C (C) = (0.6)-1 2
+ (0.75) = 0.55 and we set Bel ⇒ (¬R 1 ) = 0. In the same way, for the reverse rule R 2 :

Bel ⇐ (R 2 ) = Disb C (C) = (0.8)+1
2 -(0.25) = 0.65 and we set Bel ⇐ (¬R 2 ) = 0.

Uncertainty assessment procedure

In this section, we present our approach to uncertainty propagation from premises to the top goal of a GSN. As illustrated on Figure 6, this procedure is structured in two phases.

The first one, called modeling phase, collects expert opinions on rules, expressed with qualitative scores (Dec, Conf ), and translates them into numerical mass assignments to rules. It will be conducted by asking (2n + 2) questions to the assessor using the evaluation matrices, n being the number of premises. The first (2n) questions concern masses on elementary rules (direct and reverse). For instance, to get, respectively, the values of Bel i ⇐ (¬p i ⇒ ¬C) and Bel i ⇒ (p i ⇒ C) the expert will be asked the following questions (in case n = 2):

1. Supposing no knowledge about premise p 1 (resp. We assume that once these masses on rules are evaluated, they can be used for the considered system using the second phase explained below.

The second phase, called application phase, concerns the collection of expert data on premises. One question per premise is then formulated to the experts: considering the knowledge on the pieces of evidence (also called solutions in GSN), what is your "Decision" and "Confidence" concerning premise p i ? Grouped in a questionnaire, these (3n + 2) questions will be asked in form of matrices to be filled in by the assessor (for rules, some matrices may be pre-filled, see Figure 6). Then, these values (on rules and premises) are used to calculate the belief/disbelief in the conclusion (eqs. ( 1) and ( 2)). Finally, we may transform the resulting triple (Belief, Disbelief, Uncertainty) concerning the conclusion, to a pair (Decision, Confidence) using formulas (4) and approximate them by qualitative values. In this section, we use a portion of GSN proposed in [START_REF] Damour | Towards certification of a reduced footprint acas-xu system: A hybrid ml-based solution[END_REF] to test and validate our uncertainty propagation approach. That study proposed a hybrid architecture of a collision avoidance system for drones, Urban Air Mobility and Air Taxis with horizontal automatic resolution. It is named ACAS-X (Next-Generation Airborne Collision Avoidance System). It replaces a set of lookup tables (LUTs) (that provide anti-collision maneuvering guidance according to the speed of the two aircrafts, their relative positions, and the time until the loss of vertical separation occurs) by a neural network (NN) of much smaller size. In addition to the NN-based controller, this architecture includes a safety net which contains a portion of LUTs (already established as safe) for unsafe areas (where the NN may give results different from those of the LUTs), and a check module which controls the switch between these two sub-systems (NN and LUTs). The GSN section (figure 7) in which we are interested, argues that "G 1 : Real world situations where MLM 3 is not robust are identified and mitigated ". To demonstrate this statement, the top goal (G 1 ) is broken down into two sub-goals (G 2 ) and (G 3 ). (G 2 ) ensures that the property was correctly defined to identify all unsafe situations (G 4 ) and formally checked (G 6 ) in each of the areas (called p-boxes) into which the input space was correctly decomposed (G 5 ). (G 3 ) ensures that unsafe situations were properly mitigated via an appropriate architecture (G 7 ).

Modeling phase

Real-world situations where the MLM is not robust are identified and mitigated

(G1)

Ensure that all unsafe situations are correctly identified and mitigated

(ST)

All unsafe situations are identified

(G2) (G3)

All unsafe situations identified are mitigated

(G7)

Architecture mitigation (switch to the LUT when appropriate)

The LUT property is correctly defined

(G4)

The LUT property is checked in each p-box

(G6)

The input space (ODD) is correctly decomposed to p-boxes

(G5)

Certified development process DO178

(S4)

Formal verification results

(S3)

The validation is trivial (it consists of mathematical decomposition on the whole 3D input space)

(S2)

The validation of this property is trivial (S1) Fig. 7. Assurance Case -ML subsystem robustness [START_REF] Damour | Towards certification of a reduced footprint acas-xu system: A hybrid ml-based solution[END_REF] Table 1 groups the degrees of belief on the rules involved in this case. Following the assessment procedure above, these values are the result of a questionnaire 4 answered by a safety expert about this system. We can notice that all direct conjunctive rules receive maximal weights and the elementary rule weights for (G 1 ) and (G 2 ) are null. Thus, we deduce that this GSN represent a conjunctive type where all sub-goals are needed to support (G 1 ). As seen in [START_REF] Idmessaoud | Quantifying confidence of safety cases with belief functions[END_REF], C-Arg tends to propagate the premises that support the conclusion with the least weight, increasing along with it the uncertainty level. Thus, we can explain why we go from acceptable premises with very high confidence (G 6 , G 7 ), high confidence (G 5 ) and for sure (G 4 ) to a tolerable top goal (G 1 ) with low confidence (Dec = 0.692, Conf = 0.384). Graphs in figures 8 and 9 present, respectively, the sensitivity of decision and confidence degrees of the conclusion (G 1 ) to the sub-goals (G 4 ), (G 5 ), (G 6 ) and (G 7 ). To determine the latter, we vary the value of a premise from its minimal to its maximal value, while we keep the values of the other premises to their base values. We can notice that all values, are indeed included in the interval [0,1]. We can also notice that the pair (decision, confidence) on the goal (G 1 ) varies from "Rejectable for sure" (Dec = 0, Conf = 1) to "Tolerable with high confidence" (Dec = 0.82, Conf = 0.64). The sub-goal (G 4 ) has the lowest influence on decision and the highest influence on confidence; the opposite applies for sub-goal (G 5 ).

Conclusion

In this paper, we propose an extensive approach to the elicitation and propagation of uncertainty in a logical GSN model and report on a preliminary case study for testing our approach. However, some issues still need to be addressed. First of all, our propagation model does not consider all GSN components (such as Justification, Assumption, etc.). In addition, our elicitation model seems to encourage experts to give extreme values of (decision, confidence) so that we often end up with a conjunctive or disjunctive type. But these two types are not the only types that exist in literature. Finally, the transformation of expert opinion from quantitative to qualitative values is also a source of uncertainty. In a future work, we plan to develop a purely qualitative approach to information fusion based on [START_REF] Dubois | A possibilistic counterpart to shafer evidence theory[END_REF], and compare it to the quantitative one.

Acknowledgement

A special thanks to the authors of [START_REF] Damour | Towards certification of a reduced footprint acas-xu system: A hybrid ml-based solution[END_REF], especially to Christophe GABREAU for answering the questionnaire concerning the assessment of the GSN presented in our case study.

Fig. 4 .

 4 Fig. 4. Evaluation matrix

Fig. 5 .

 5 Fig. 5. Extreme assessments (black dot)

6 -1 2 + 2 -

 622 0.25 = 0.05, Disb(p 1 ) = 0.6+1 0.25 = 0.55 and U ncer(p 1 ) = 1 -Bel(p 1 ) -Disb(p 1 ) = 0.4 for the amount of ignorance.

Example 2 .

 2 Consider the case of Figure 2: -Direct rule (R 1 : p ⇒ C): Assuming Dec(p) = 1, expert assigns "Tolerable with high confidence" to C: Dec(C) = 0.75, Conf (C) = 0.6 -Reverse rule (R 2 : ¬p ⇒ ¬C): Assuming Dec(p) = 0, expert assigns "Opposable with very high confidence" to C: Dec(C) = 0.25, Conf (C) = 0.8

  p 2 ) : (Dec = 0.5, Conf = 0) and minimal Dec value (rejectable for sure) in premise p 2 (resp. p 1 ): (Dec = 0, Conf = 1), what is your Decision/Confidence in the conclusion? 2. Supposing no knowledge about premise p 1 (resp. p 2 ): (Dec = 0.5, Conf = 0) and a maximal Dec value (acceptable for sure) concerning premise p 2 ( resp. p 1 ) : (Dec = 1, Conf = 1), what is your Decision/Confidence in the conclusion? The additional two questions concern the conjunctive rules (resp. reverse and direct): 3. Supposing minimal Dec value (rejectable for sure) concerning both premises p 1 , p 2 : (Dec = 0, Conf = 1), what is your Decision/Confidence in the conclusion? 4. Supposing maximal Dec value (acceptable for sure) concerning both premises p 1 , p 2 : (Dec = 1, Conf = 1), what is your Decision/Confidence in the conclusion?

Fig. 6 .

 6 Fig. 6. Schema of the assessment framework for safety argument 5 Case study

G1 (i = 2 , n = 3 ) 1 OpposableFig. 8 .Fig. 9 .

 23189 Fig. 8. Decision sensitivity on the top goal G1

Table 1 .

 1 Elicited belief degrees on rules

	Goal (Gi)	Belief degree on rules

Machine learning Model.


Christof Dues

Marius J Müller

Sangam Chatterjee

Claudio Attaccalite

Simone Sanna
email: simone.sanna@theo.physik.uni-giessen.de

Nonlinear optical response of ferroelectric oxides: first-principles calculations within the time-domain and the frequency-domain

Keywords: Ferroelectrics, DFT, SHG, THG, spectroscpy, LiNbO3, LiTaO3, KNbO3

come

Nonlinear optical response of ferroelectric oxides: First-principles calculations within the time and frequency domains

I. INTRODUCTION

Nonlinear optical phenomena are at the basis of many technological applications [1][START_REF] Brown | The Optics Encyclopedia: Basic Foundations and Practical Applications[END_REF][START_REF] Shen | [END_REF][START_REF] Butcher | The elements of nonlinear optics[END_REF].

They range from, e.g., well established frequency converters [START_REF] Tucker | [END_REF] to exciting novel approaches for white light generation [START_REF] Alfano | The supercontinuum laser source: fundamentals with updated references[END_REF][START_REF] Rosemann | [END_REF][8][9][10][11][12], which is currently understood as a combination of several nonlinear optical effects of different order [13,14]. The demand for a theoretical description of nonlinear effects and the accurate calculation of the related susceptibilities has grown parallel to the employment of nonlinear optical materials in different devices.

Unfortunately, the ab initio modelling of the nonlinear optical response of a medium is one of the most challenging tasks in theoretical physics. While the linear-response optical properties are successfully calculated from first principles within Green's function theory [15], this is not the case for nonlinear optical susceptibilities. The state-of-the-art approach for the computation of the linear optical response combines band structures calculated in many-body perturbation theory [16] (e.g., by G 0 W 0 ) to include quasiparticle effects, with the solution of the Bethe-Salpeter equation to account for the electron-hole attraction [17].

Within this approach, it is difficult to elaborate expressions for the nonlinear optical susceptibilities which include the many-body effects [18]. Accurate second-order susceptibilities can be calculated, e.g., on the basis of the electronic wavefunctions from a combination of density functional and a perturbation theory [19,20]. However, the complexity of the expressions derived within perturbation theory quickly grows with the perturbation order [21], and makes nonlinearities of higher orders de facto inaccessible [22]. The applications are generally limited to smaller systems such as periodic crystals [18]. Up to our knowledge, a single attempt to solve the Bethe-Salpeter equation for second harmonic generation can be found in the literature [21].

The calculation of the nonlinear optical susceptibilities can be also performed in the time domain from the dynamical polarization [23][24][25][26]. The response of a medium to a timedependent electrical field may be expanded into a power series:

P α (ω) = β χ (1) αβ (-ω; ω)E β (ω) + βγ χ (2) αβγ [-ω = -(ω + ω ); ω , ω ]E β (ω )E γ (ω ) + βγδ χ (3) αβγδ [-ω = -(ω + ω + ω ); ω , ω , ω ]E β (ω )E γ (ω )E δ (ω ) + • • • .
(1) Thus, if the dynamical polarization is calculated, e.g., by numerical integration of the equations of motion in presence of a laser field, optical susceptibilities of virtually any order can be obtained. Moreover, several nonlinear phenomena such as sum-and differencefrequency generation or four-wave mixing can be calculated simultaneously, as they are described by the same equations of motion. Thereby, the dynamical polarization must be calculated as a geometric Berry phase as described in the modern theory of polarization [27],

if periodic boundary conditions are applied; this is generally required for the description of crystalline solids. In the past, one of the authors has presented a practical implementation of this approach [28], in which the equations of motion are derived following the scheme introduced by Souza et al., [29] based on the generalization of the Berry phase to the dynamical polarization [30,31].

In the present work, we calculate and the nonlinear optical susceptibilities of different ferroelectric oxides both within the frequency-domain and within the time-domain approach and corroborate the results with cooresponding experimental data. In a first step, we calculate the second harmonic generation (SHG) spectrum of LiNbO 3 up to energy of the incoming photons of 6 eV, based both on the momentum matrix approach and on the time evolution of the polarization. The excellent agreement of the spectra calculated by the different methods validates both approaches against each other, furthermore confirming earlier theoretical and experimental results [32]. In a second step, we calculate the SHG spectrum as well as the third harmonic generation (THG) spectra of LiNbO 3 , LiTaO 3 as well as KNbO 3 , and, if possible, compare them with the experimental results and calculations from the literature. Furthermore, we investigate the role of quasiparticle effects on the nonlinear optical response by making use of the computationally efficient independent particle approximation and incorporating quasiparticle effects by a previously determined scissors-shift. In addition, we show that the major resonances in the nonlinear spectra can be explained by multi-photon processes within the fundamental bandgap.

II. METHODOLOGY

A. Second harmonic generation in the frequency-domain: the momentum matrix approach In this work, we employ the perturbative approach explained in detail in Refs. [20,32].

The method is based on the calculation of the momentum matrix elements

p α nm ( k) = n k pα m k (2)
where m k and n k are two Bloch states.

Considering the transition energies ω n,m = ε m k -ε n k between the states m and n at the reciprocal space point k, and the notation (derived from the anticommutator of the momentum operators in the cartesian directions β and γ)

p β ml ( k)p γ ln ( k) = 1 2 p β ml ( k)p γ ln ( k) + p γ ml ( k)p β ln ( k) , (3)
the SHG susceptibility is calculated for a complex frequency ω = ω + iη, in which the small positive imaginary part iη adiabatically switches the electromagnetic field on. In our work, we set η to 0.2 eV for all the calculations.

The expression for the SHG susceptibility reads

χ (2) αβγ (-2ω; ω, ω) = - ie 3 ω3 2 m 3 V k nml 1 ω mn ( k) -2ω ×   f nl ( k)p α nm ( k) p β ml ( k)p γ ln ( k) ω ln ( k) - ω + f ml ( k)p α nm ( k) p γ ml ( k)p β ln ( k) ω ml ( k) - ω   . (4)
In this expression, the eigenvalues ε n k can be evaluated either within the independent Partice Approximation (IQA, i.e., DFT with (semi)local xc-potentials) or within the independent quasiparticle approximation (IQA, e.g., from G 0 W 0 calculations).

Equation 4 can be divided into a two-band term

χ (2),two αβγ (-2ω; ω, ω) = - ie 3 2 m 3 V k nm    16f nm ( k)p α nm ( k) ∆ β mn ( k)p γ mn ( k) ω mn ( k) 4 ω mn ( k) -2ω - f nm ( k)p α nm ( k) ∆ β mn ( k)p γ mn ( k) ω mn ( k) 4 ω mn ( k) -ω    (5)
and a three-band term

χ (2),three αβγ (-2ω; ω, ω) = - ie 3 2 m 3 V k nml n =m =l p α nm ( k) p β ml ( k)p γ ln ( k) ω ln ( k) -ω ml ( k) ×    16f nm ( k) ω mn ( k) 3 ω mn ( k) -2ω + f ml ( k) ω ml ( k) 3 ω ml ( k) -ω + f ln ( k) ω ln ( k) 3 ω ln ( k) -ω    . (6)
In equation 5, the matrix elements of the intraband transitions ∆ β mn ( k) = p β mm ( k)p β nn ( k) are calculated as

p β mm ( k) = m e ∂ k β ε n ( k). (7)
B. Hyperpolarazibilites in the time domain: Time evolution of the polarization

For the calculation of the nonlinear optical response in time-domain from the time evolution of the polarization, we employ the procedure described in Ref. [33], which we briefly outline in the following. The starting point are the (zero-field) Kohn-Sham equations of the form, Ĥ0,IPA = -

2 2m i ∇ 2 i + VeI + VH [ρ 0 ] + Vxc [ρ 0 ] (8)
where V eI is the electron-ion interaction and V H and Vxc are the Hartree and the exchangecorrelation potentials, respectively. Improvements beyond the IPA might be introduced at this step to consider quasiparticle shifts or effects originating from the response of the effective potential to density fluctuations in a time-dependent screened Hartree-Fock manner.

The corresponding Hamiltonian is denoted in this case as Ĥ0 instead of Ĥ0,IPA .

The laser excitation of frequency ω L of the form E(t) = E 0 sin(ω L t) is used to define the non hermitian field coupling operator

ŵ k = ief 4π m 3 α=1 ( a α • E)N kα σ=± σ ṽ k σ α ,m v k,m . (9)
In this expression,

ṽ k ± α ,n = m (S( k, k ± α )) mn v k ± α ,m (10)
and S mn are matrix elements defined by equation 12.

Replacing the coupling operator ŵ k with the hermitian form ŵ k + ŵ † k , the equations of motion are obtained:

i d dt v k,m = Ĥ0 + ŵ k ( E) + ŵ † k ( E) v k,m . (11)
To account for quasiparticle effects, a scissor operator can be included in equation 11, which modifies the eigenvalues of the Hamiltoniana Ĥ0 without modifing the corresponding eigenvectors. The equations of motion are then solved to obtain the lattice-periodic part v k,n of the Bloch states n k . From these, the overlap integrals

S mn ( k, k + q α ) = v k,m v k+ qα,n (12)
are calculated, which are the elements of the matrix S, from which the time dependent polarization is calculated as

P α = - ef 2πv a α N k ⊥ α k ⊥ α   N kα -1 i=1 Tr ln S( k i , k i + q α )   . (13)
In a final step, the nonlinear susceptibilities are obtained by postprocessing the polarization in a signal analysis procedure as described in Ref. [33].

C. Computational parameters and electronic groundstate

The calculation of the nonlinear optical response in frequency-domain from the momentum matrix elements is performed within the density functional theory (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP, Version 5.4.4 [34,35]). The exchange-correlation potential in the formulation of Perdew, Burke, and Ernzerhof [36] and projector-augmented waves potentials [37] are employed, that include the 2s 1 , valence electrons in case of lithium, 2s 2 2p 4 in the case of oxygen, and 4s 2 4p 6 4d 4 5s 1 for niobium, respectively. The basis set for the expansion of the wavefunctions contains plane waves with kinetic energy up to 400 eV. The unit cell for the simulation of LiNbO 3 (and LiTaO 3 ) is rhombohedral and consists of two formula units. In this case, the integration in the reciprocal space is performed on a 6 × 6 × 6 Monckhorst-Pack grid [38], which reflects the symmetry of the unit cell and corresponds to 38 k points in the irreducible Brillouin zone.

With these computational parameters, the atomic structure of lithium niobate closely matches that of earlier calculations [32,39,40] and reproduces the experimental values within 1 % [START_REF] Räuber | Chemistry and physics of lithium niobate[END_REF][START_REF] Volk | Lithium Niobate, Defects, Photorefraction and Ferroelectric Switching[END_REF]. The corresponding structure is shown in Fig. 1 together with partial charge densities associated to the valence and (at Γ degenerate) conduction band edges.

The electronic band structure of LiNbO 3 is shown in Fig. 2. It features a very flat dispersion of both valence and conduction states, which is characteristic of the material [32,43,44]. The indirect (direct) fundamental electronic bandgap amounts to 3.52 eV (3.42 eV) which also is in agreement with earlier calculations [32,39].

The rhombohedral unit cell of LiNbO 3 contains 64 electrons, which occupy the lowest 32 Kohn-Sham-states. In total, 256 bands are employed for the calculation of the linear optical properties because of the slow convergence of the real-part of the dielectric function with respect to the number of conduction bands. The resulting dielectric function (absolute value, along with the real and imaginary part) is shown exemplarily for the ε zz component in Fig. 3 (lhs). The structured peak of the imaginary part at about 4 eV represents the most important spectral feature in agreement with previous results [32]. The onset of the optical absorption thus corresponds to the bandgap energy.

The calculation of the nonlinear optical response in the time-domain from the timeevolution of the dynamical polarization is performed on the basis of the electronic groundstate calculated within the DFT as implemented in the QUANTUM ESPRESSO code [45,46] (Version 6.5) using also in this case the Perdew-Burke-Ernzerhof [36] (PBE) functional. The set of SG15 ONCV pseudopotentials [47,48] is used to describe the electron-core interaction of the involved atoms. The wavefunctions are expanded in a plane wave basis

│χ (2) zzz (-2ω;ω,ω)│ [pm / V] ℏω [eV] │χ (2) │ Im[ χ (2) ]
Re[ χ (2) ] FIG. 3. Left-hand side: zz component of dielectric tensor calculated for ferroelectric LiNbO 3 within DFT-PBE in the IPA. Right-hand side: zzz component of the SHG tensor of LiNbO 3 calculated on the same footing. For better comparison, the real and imaginary part of χ (2) are also multiplied by (-1) (dashed lines).

up to a cut-off energy of E = k The YAMBO [49] code is employed to calculate the optical properties. To this end, the The substitution of the alkali metal Li (group 1, period 2) with the heavier, isovalent K (group 1, period 4) leads to a fundamentally different crystal, KNbO 3 . KNbO 3 crystallizes, with decreasing temperature, in a cubic, tetragonal, orthorhombic, and monoclinic phase [53]. In this work, we focus on the tetragonal phase (space group P 4mm, c/a ratio of 1.023), as it is computationally convenient and because recent studies have shown that it features the largest SHG coefficients in comparison to the other phases [54]. The DFT-PBE calculated atomic structure reproduces within1 % the experimentally determined structural parameters [50,55].

Similarly to LiNbO 3 and LiTaO 3 , the valence band top originates from the anionic sublattice and has O character, while the cations (Nb) determine the conduction band bottom, as shown in Fig. 4.

For the calculation of the optical response, a very dense 16×16×16 k-point mesh (corresponding to a total of 720 k-points) is necessary. For the time evolution of the wavefunctions of the 40 electrons contained in the unit cell, a Kohn-Sham basis made of 15 occupied states and 27 unoccupied states is considered. The nonlinear susceptibilities are calculated sampling the energy axis with points at a distance of 0.08 eV.

D. Experimental setup

All investigated samples are commercially available X-cut wafers (Surface Net GmbH, Germany). We use a 50 fs pulse duration 5 kHz repetition rate Ti:Sa amplifier to drive an optical parametric amplifier including mixing stages to generate sub-100 fs pulses from 800 nm to 1580 nm. The beam diameter is reduced by two concave mirrors and split into two beams using a symmetric Michelson interferometer. SHG is measured by placing the crystal in the unfocused beam under normal incidence. Lock-in technique on Si photodiodes are used for detection. The fundamental laser is suppressed by the appropriate dielectric shortpass filters. The data are corrected against a simultaneously acquired reference on z-cut quartz for which the literatrue absolute value of d 11 = 0.3 pm/V at 1064 nm is assuemd [56].

We use Miller's constant-delta condition [57] with a Miller delta of δ 11 = 1.328 × 10 2 m/C to account for the dispersion of the second-order nonlinear coefficient of the quartz reference, as was already performed, for example, in Ref. [54].

III. RESULTS AND DISCUSSION

We start our discussion with the comparison of the SHG spectra obtained with the momentum matrix approach and with the time-evolution of the dynamical polarization. As the focus is partially set on the comparison of the two methods, results within the IPA are presented. Many-body effects are considered in a second step, based on the quasiparticle shifts reported in Refs. [32,54].

After expansion of the wavefunctions from the irreducible part to the full Brillouin zone, and the calculation of their derivatives, the second order optical susceptibility

χ (2) 2 = χ (2) 2 + χ (2) 2 (14)
can be estimated in the frequency domain as a sum over momentum matrix elements.

We find that the two-band contributions are negligible and the three-band contributions play the crucial role. Fig. 3 (lhs) shows the result of this procedure exemplarily for the χ

zzz (-2ω, ω, ω) tensor component. The most prominent feature is the first peak at 1.9 eV with an intensity of 183 pm/V confirming the strong optical non-linearity of LiNbO 3 . In the static limit, χ

zzz (0) is as high as 42 pm/V. The calculated spectrum is in overall good agreement with the results reported in Ref. [32].

The position of the resonances in the SHG spectrum can be understood on the basis of the optical absorption. To this end, we show in Fig. 5 (lhs) the absolute values of the SHG coefficients χ

(2) zzz and the imaginary part of the dielectric tensor ε zz (ω) as well as ε zz (ω/2). The latter is characterized by spectral signatures roughly at energies for which also the SHG spectrum shows the main peaks. This means that the SHG spectral features might be understood as two-photon processes.

In a second step, the SHG spectrum of LiNbO 3 is calculated from the time evolution of the polarization. After calculation of P (t) and expansion in a power row up to the 6 th order, a Fourier analysis yields the spectrum shown in Fig. 5 (rhs). For energies within the fundamental bandgap, two peaks at 1.9 eV and 2.5 eV represent the most prominent signatures. Both peaks match the spectral features calculated in the frequency domain, both in peak position and intensity.

The direct comparison of the LiNbO 3 SHG spectra calculated with the two approaches is shown in Fig. 6. In the investigated frequency range, the agreement is excellent concerning both the position and intensity of the spectral features. The main difference between the two spectra is represented by the magnitude of the deeps as calculated in the time-domain, which might be due, however, to a different choice of the damping parameter which determines the resonance widths and to the use of different pseudo-potentials. The agreement between the results obtained with two fundamentally different approaches is by no means a matter of course, and represents a major mutual validation of the two methods and their implementations.

The major advantage of the time-domain approach over the frequency-domain approach is the simultaneous calculation of the higher harmonics. Fig. 7 shows exemplarily the zzzz component of the THG tensor of LiNbO 3 . It is characterized by a dominant peak at 1.4 eV with a very strong intensity of about 1.7•10 5 pm 2 /V 2 , a minor feature at 2.6 eV and a small shoulder at 3.9 eV. A comparison with the corresponding component of the dielectric function shows clearly that the first peak is due to three-photon processes, the second feature to two-photon processes, and the minor shoulder with roughly the bandgap energy is due to

zzz component of the second-order optical susceptibility calculated for the investigated crystals, while Fig. 9 shows a comparison of calculated and measured SHG spectra for LiNbO 3 and LiTaO 3 . Notably, the calculations correctly predict the spectral depence of the measured values in the experimentally accessible range from 0.78 eV to 1.55 eV and also the magnitude agrees exeptionally well, i.e., within a factor of two.

The SHG coefficients of LiNbO 3 and LiTaO 3 are similar concerning both peak position and intensity, although the LiTaO 3 spectrum is blueshifted due to the higher bandgap. We observe, however, that the energy shift between the two spectra (the distance between the first peak of each spectrum is of about 0.45 eV) is higher than the difference between the two bandgap energies (about 0.15 eV). This suggests that not only the band edges but also additional states close to the gap contribute to this spectral feature. As the bandstructure of the two ferroelectrics are very similar, the energy shift might be explained by the different lattice parameters of the two compounds. In particular, the two spectra differ in the low energy region before the first peak. While LiNbO 3 is characterized by a continuous increase, the SHG coefficient of LiTaO 3 remains below 15 pm/V till 1.8 eV, followed by a steep gradient which brings the optical susceptibility to 101 pm/V at 2.4 eV. This behavior, which also holds true when quasiparticle effects are included in the calculation, might explain why lower SHG coefficients are measured for LT in the visible range.

The SHG spectrum of KNbO 3 is fundamentally different, from those of LiNbO 3 or LiTaO 3 .

In the long-wavelength limit, the zzz component of the SHG coefficient (ca. 40 pm/V) is much larger than in the case of LiNbO 3 or LiTaO 3 . After a steep gradient beginning at 1.5 eV, the SHG coefficient reaches a global maximum of 330 pm/V at 2.6 eV. Thus, the zzz component of the SHG coefficient as calculated within the time-domain is in qualitative agreement with the calculations presented in Ref. [54] and based on the momentum-matrix approach.

Figure 8 (rhs) shows how the χ the investigated crystals is affected by quasiparticle effects. Besides an expected blueshift of half the value of the shifts (indicated in the figure), the quasiparticle effects lead to an overall reduction of the signal intensity. A similar effect has been predicted by Riefer et al. [32] for LiNbO 3 . Nevertheless, no major redistribution of the spectral weights is observed, and the relative differences between the spectra closely mirror those described for the IPA spectra.

The χ

(3) zzzz component of the THG spectra (reported in Fig. 10, lhs) as calculated in the IPA for the investigated oxides has a rather similar behavior for all crystals and does not feature the differences observed for the corresponding SHG tensor component. All spectra are dominated by main THG peak between 1 eV and 2 eV, whose onset mirrors the order of the bandgap energies of KNbO 3 , LiNbO 3 and LiTaO 3 . Yet, the third-order energy differences amount to one third of the bandgap energy, so that the three main peaks are rather close to each other. The height of the maxima grows from KNbO 3 (115000 pm 2 /V 2 ) to LiNbO 3 (169000 pm 2 /V 2 ) which is the reverse order of the SHG maxima. On contrary, the order of the THG activity in the long-wavelength limit mirrors that of the corresponding SHG component, with the highest value for KNbO 3 (15500 pm 2 /V 2 ) and the lowest for LiTaO 3 (11500 pm 2 /V 2 ).

When quasiparticle effects are considered (Fig. 10, rhs) the intensity of the THG signatures is drastically reduced by about 50% of the IPA value. Yet, LiNbO 3 and LiTaO 3 do not undergo large redistributions of the relative spectral weights. However, due to the quasiparticle-induced shifts the THG spectra of LiNbO 3 and LiTaO 3 become almost indistinguishable in the visible range.

Differently than in the case of LiNbO 3 and LiTaO 3 , the quasiparticle effects strongly modify the KNbO 3 THG spectrum. The intensity of the main peak is somewhat less affected than in the case of the other ferroelectrics, and, more important, the intensity of the peak at 3.4 eV is strongly increased to 47000 pm 2 /V 2 . This is an unusual feature, as the general decrease of the optical susceptibility upon widening of the fundamental gap is a known effect based on the sum rules for the harmonic susceptibilities [58,59].

We remark that neither the crystal local fields nor the electron-hole interaction have been considered in our calculations. While local fields are expected to reduce the optical nonlinearities, excitonic effects typically enhance them [60,61], so that some error cancelation might be expected. Yet exciton effects are in some cases not strong enough enought to compensate the quasi-particle effects in THG [62].

IV. CONCLUSIONS

The nonlinear optical response of the ferroelectric oxides LiTaO Due to the efficiency of the time-domain approach and the possibility to calculate higher harmonics, this method seems particularly suitable to investigate complex nonlinear processes from first principles. This might be particularly important, e.g., to understand the mechanisms underlying the generation of a high directional light supercontinuum upon infrared irradiation recently observed in organotetrel molecular crystals [63].

FIG. 1 .

 1 FIG. 1. Rhombohedral LiNbO 3 unit cell with Li, Nb and O in grey, white and red, respectively. Partial charge densities of the electronic groundstate are shown as isosurfaces (0.011 eV/ Å3 ). The topmost occupied valence state (left part) and the lowermost unoccupied conduction states (central and right part) strongly resemble the O 2p and Nb 4d orbitals, respectively.

FIG. 2 .

 2 FIG. 2. Electronic band structure of LiNbO 3 . The energies are given relative to the valence band maximum. The direct (indirect) fundamental electronic bandgap amounts to 3.52 eV (3.42 eV).

2 /

 2 (2m) = 80 Ry ≈ 1090 eV. Integration of the reciprocal space has been performed with a Γ-centered 10 × 10 × 10 k-Mesh, which consists of 172 individual k points.

  DFT (single particle) wave functions calculated with QUANTUM ESPRESSO are imported to build a Kohn-Sham basis set consisting of 22 topmost valence bands and the 43 lowest conduction bands. These values have been carefully chosen on the basis of the convergence tests shown in Appendix 1. Test calculations of the linear optics show an indirect (direct) fundamental electronic bandgap of 3.45 eV (3.52 eV). The band gap energy values, the band structure and the dielectric function are in overall very good agreement with the VASP calculations. Furthermore, additional occupied or unoccupied states only lead to minor changes in the imaginary part of the dielectric constant. Six additional states modify the dielectric function by less than 1 %.The optical susceptibilities are explicitly calculated for ω representing 64 energy values between 0.2 eV and 6 eV, as well as for an incident electric field with a field strength of 10 14 W/m 2 and a damping of 0.2 eV. The Crank-Nicholson algorithm is used to calculate P (t) for 5364 time steps 0.01 fs apart. Taking dephasing of eigenmodes introduced by the sudden switch-on of the electric field into consideration, only the last 20.68 fs of the polarization are eligible for frequency analysis, which are used to extract P (2ω) and P (3ω). Next, this result is converted into χ(2) (-2ω; ω; ω) and χ(3) (-3ω; ω, ω, ω) for the outbound radiation corresponding to 2 ω and 3 ω, respectively, i.e. the harmonic generation of inbound radiation with energy ω of second and third order, respectively.In order to explore the effect of chemical variations in the cationic sublattice on the nonlinear optical response of LiNbO 3 , the ferroelectric oxides LiTaO 3 and KNbO 3 were modeled. The substitution of the transition metal Nb (group 5, period 5) with the heavier, isovalent Ta (group 5, period 6) leads to the formation of an isomorph crystal, ferroelectric LiTaO 3 , which likewise crystallizes into the R3c space group and can be modeled by a rhombohedral unit cell containing two formula units. The material as well as its properties and applications are in general very similar to LiNbO 3 . However, the Curie[50] and melting (1923 K vs 1526 K)[START_REF]Korth kristalle gmbh[END_REF] temperatures, as well as the coercive fields (17 kV/cm vs 40 kV/mm)[START_REF] Kim | [END_REF] are lower, and the optical nonlinearities in the visible range less pronounced. In our simulations, LiNbO 3 and LiTaO 3 are modeled with identical numerical parameters. Besides a different value of the fundamental electronic bandgap (3.56 eV), LiTaO 3 has a very similar, likewise flat electronic bandstructure to LiNbO 3 , whereby the valence band maximum and the conduction band minimum have O 2p and Ta 5d character, respectively. The dielectric function of LiTaO 3 strongly resembles that of LiNbO 3 (see Fig.3, rhs) although all spectral features are slightly shifted to higher energies due to the larger bandgap energy.

FIG. 4 .

 4 FIG. 4. Tetragonal structure of KNbO 3 , with K, Nb and O atoms depicted in yellow, white and red, respectively. Partial charge densities corresponding to the highest valence band (left-hand side) and the lowest conduction band (right-hand side) are shown. Isosurfaces correspond to a charge density of 0.017 eV/ Å3 .

FIG. 5 .

 5 FIG. 5. Absolute values of the zzz component of the SHG tensor of LiNbO 3 calculated by the momentum matrix elements (left-hand side) and by the time evolution of the polarization (righthand side). Each of them is compared with the imaginary part of the dielectric function ε zz (ω) and its energy-scaled counterpart ε zz (ω/2).

FIG. 6 .

 6 FIG. 6. Absolute values of the SHG susceptibility of LiNbO 3 (the zzz component is representatively shown) as calculated in the time-domain (dashed blue line) and in the frequency domain (solid red line).

Figure 8 (

 8 Figure 8 (lhs) shows a comparison of the χ

FIG. 7 .

 7 FIG. 7. Absolute values of the zzzz component of the THG tensor of LiNbO 3 as calculated from the time evolution of the polarization. The imaginary part of the corresponding component of the dielectric tensor zz (ω) and its energy-scaled counterparts zz (ω/2) and zz (ω/3) are reported for comparison.

( 2 )FIG. 8 .

 28 FIG. 8. Absolute values of the zzz component of the SHG tensor in the IPA (left-hand side) and in the independent quasiparticle approximation (right-hand side) calculated for different ferroelectric oxides.

FIG. 9 .

 9 FIG. 9. Absolute values of the zzz component of the SHG tensor of LiNbO 3 and LiTaO 3 as calculated by DFT in the IPA (left-hand side) and measured (right-hand side).

3 ,

 3 LiNbO 3 , and KNbO 3 is investigated from first principles and compared with specifically performed measurements of the SHG spectra of LiTaO 3 and LiNbO 3 . To this end, two different approaches, based on the calculation of the momentum-matrix elements and of the dynamical polarization in its Berry phase formulation, are employed, respectively. While the second order optical susceptibility is calculated in the frequency domain, second-order and third-order nonlinearities are calculated within the time-domain. The excellent agreement between the two approaches concerning the calculated second order optical susceptibility validates both methods against each other, and allows to reproduce the experimental results as well as earlier results where available. KNbO 3 shows the highest SHG coefficients, while all three materials are comparable concerning THG. The main spectral features of all investigated ferroelectrics can be explained by parametric processes and multi-photon resonances. The consideration of quasiparticle effects blueshifts the spectra without affecting the relative differences between them, and reduces the optical nonlinearities. This reduction is particularly severe concerning THG, which is reduced by about 50% with respect to the IPA calculations. The only exception is an anomalous peak at 3.4 eV in the THG of KNbO 3 , which gains in intensity upon consideration of quasiparticle effects. Our calculations furthermore reveal that LiNbO 3 and LiTaO 3 have distinct SHG spectra but almost indistinguishable THG spectra. It must be therefore possible by means of LiNb x Ta 1-x O 3 solid solutions to tune the crystal's second order optical susceptibility without affecting the third harmonic generation.

ACKNOWLEDGEMENT

  This work is supported by the German Science Foundation (DFG) through the research groups FOR2824 and FOR5044. Calculations for this research were conducted on the Lichtenberg high performance computer of the TU Darmstadt and at the Höchstleistungrechenzentrum Stuttgart (HLRS). The authors furthermore acknowledge the computational resources provided by the HPC Core Facility and the HRZ of the Justus-Liebig-Universität Gießen.V. APPENDIX: CONVERGENCE OF NUMERICAL PARAMETERS[1] R. W. Boyd, Nonlinear optics (Academic press, 2020).

FIG. 11 .

 11 FIG. 11. Lhs: Band structure of LiNbO 3 , in which the energy relative to the valence band maximum is given on the left axis and the corresponding band indices (cf. rhs) are given on the right axis. Rhs: Resulting SHG spectra of LiNbO 3 calculated with different Kohn-Sham basis sets. In the bottom (top) half of the figure, the lower (higher) boundary of the band range of the Kohn-Sham basis is varied, while the higher (lower) boundary is fixed. In each case, the corresponding band indices are given in the key. These calculations are all performed with a 4 × 4 × 4 grid of k-points.VB and CB stay for valence and conduction bands, respectively.

12 FIG. 12 .

 1212 FIG. 12. Resulting SHG spectra of LiNbO 3 calculated with different number of subdivisions regarding the discretization of the Brillouin zone. All calculations are performed with a Kohn-Sham basis consisting of bands 11 -75.