Qiang Tang
email: q.tang@utwente.nl

Arjan Jeckmans
email: a.j.p.jeckmans@utwente.nl

Towards a Security Model for Computational Puzzle Schemes

Keywords: computational puzzle, deterministic computation, parallelization, security model, repeated squaring computational puzzle, determinable difficulty, parallelization

In the literature, computational puzzle schemes have been considered as a useful tool for a number of applications, such as constructing timed cryptography, fighting junk emails, and protecting critical infrastructure from denial-of-service (DoS) attacks. However, there is a lack of a general security model for studying these schemes. In this paper, we propose such a security model and formally define two properties, namely the determinable difficulty property and the parallel computation resistance property. Furthermore, we prove that a variant of the RSW scheme, proposed by Rivest, Shamir, and Wagner, achieves both properties.

Introduction

A computational puzzle scheme [START_REF] Dwork | Pricing via processing or combatting junk mail[END_REF]19,[START_REF] Juels | Client puzzles: A cryptographic countermeasure against connection depletion attacks[END_REF][START_REF] Rivest | Time-lock puzzles and timed-release crypto[END_REF] enables a prover to prove to a verifier that a certain amount of resources has been dedicated to solve a puzzle. To motivate our discussions, consider the following two toy examples.

1. Let H be a one-way hash function. The verifier selects a random number x and set the puzzle to be H(x), while the prover searches for a number x , where H(x) = H(x ), as the puzzle solution. The verifier

H(r i ) = h i (1 ≤ i ≤ d).
In the first scheme, the verifier only needs to compute one hash in order to verify the solution, regardless of the amount of computation that the prover needs to perform. In contrast, in the second scheme, the verifier needs to perform the same amount of computation as that of the prover in order to verify the solution. Intuitively, the first scheme seems to be the preferred one, but it still has some drawbacks.

• First, the amount of computation the prover needs to perform is probabilistic. Suppose x is chosen from the domain [1, N]. In the best case the prover only needs to compute one hash, while in the worst case it needs to compute N hashes until find the right answer. As a consequence, the verifier never knows what is the exact amount of computation required to solve a puzzle, although the average is N 2 hashes.

• Second, if a prover has access to more than one computer, then it can speed up the puzzle solving process by having all of them work in parallel. For example, client A, which has access to x computers, could be x -1 times faster in finding a solution than client B, which has access to one computer. In practice, it is very difficult for a verifier to determine the amount of computing resources a prover can access, especially in the presence of malicious provers which control a large number of Zombie computers. In some cases such as timed cryptography [START_REF] Boneh | Timed commitments[END_REF][START_REF] Garay | Timed release of standard digital signatures[END_REF][START_REF] Rivest | Time-lock puzzles and timed-release crypto[END_REF], this will create an unfair situation for different provers. Arguably, in many other cases such as fighting against DoS attacks [START_REF] Juels | Client puzzles: A cryptographic countermeasure against connection depletion attacks[END_REF], the parallelism property of a computational puzzle scheme may also be undesirable.

Contribution

Most of existing computational puzzle schemes, as surveyed in Section 2, possess similar drawbacks. In fact, there is a lack of a general security model for formally studying relevant properties. In this paper, we propose a security model for computational puzzle schemes in general. In particular, in correspondence to the above drawbacks, we formally define two properties, namely determinable difficulty and parallel computation resistance.

• The determinable difficulty property implies that the verifier can precisely determine the required resource required from the prover in solving a puzzle.

• The parallel computation resistance property implies that the prover cannot accelerate the puzzle solving process by letting more than one computer work in parallel.

These two properties describe the fundamental characteristics of computational puzzles. This is the first model of this kind, though there are security models for some specific categories of computational puzzle schemes (e.g. [START_REF] Chen | Security notions and generic constructions for client puzzles[END_REF]). In addition, we also provide discussions on some other properties such as computation disparity between the verifier and a prover, puzzle hardness granularity, and puzzle statefulness. We prove that a variant of the RSW computational puzzle scheme, proposed by Rivest, Shamir, and Wagner, achieves the determinable difficulty and parallel computation resistance properties. This not only shows that the defined security properties are achievable but also solves the open issue about analyzing the security of the well-known RSW scheme in a rigorouslydefined security model. To our knowledge, most existing computational puzzles do not achieve the parallel computation resistance property.

Organization

The rest of the paper is organized as follows. In Section 2, we present a brief literature review on computational puzzle schemes and present a formal definition. In Section 3 we propose a general security model and define the properties. In Section 4, we present a variant of the RSW computational puzzle scheme and prove its security. In Section 5, we conclude the paper.

Literature Review

Merkle [START_REF] Merkle | Secure communications over insecure channels[END_REF] was the first to introduce the notion of puzzle which led to the invention of public key cryptography. In this context, the puzzle is required to be unsolvable by any polynomial-time entity. Dwork and Naor [START_REF] Dwork | Pricing via processing or combatting junk mail[END_REF] proposed the concept of pricing function to combat junk emails. Rivest, Shamir, and Wagner [START_REF] Rivest | Time-lock puzzles and timed-release crypto[END_REF] proposed the concept of timed-lock puzzle, which serves as a tool to realize the concept of timed-release crypto. Juels and Brainard [START_REF] Juels | Client puzzles: A cryptographic countermeasure against connection depletion attacks[END_REF] proposed the concept client puzzle and suggested to use it to prevent Denial of Service (DoS) attacks. Regardless of the different notations, pricing function, timed-lock puzzle, proof of work, and client puzzle share the same characteristic: they can be regarded as another type of puzzle (different from that of Merkle [START_REF] Merkle | Secure communications over insecure channels[END_REF]) which is moderately hard in the sense that a polynomial-time entity can successfully find a solution by spending a certain amount of resources. Without loss of generality, we use the term computational puzzle to refer to the second type of puzzle.

Roughly speaking, there are two types of computational puzzle schemes. One type is CPU-bound, where the computation is measured by the amount of CPU cycles needed to solve a puzzle. Some examples are those in [START_REF] Back | Hashcash -amortizable publicly auditable cost functions[END_REF][START_REF] Chen | Security notions and generic constructions for client puzzles[END_REF][START_REF] Dwork | Pricing via processing or combatting junk mail[END_REF]19,[START_REF] Juels | Client puzzles: A cryptographic countermeasure against connection depletion attacks[END_REF][START_REF] Rivest | Time-lock puzzles and timed-release crypto[END_REF][START_REF] Tritilanunt | Toward non-parallelizable client puzzles[END_REF][START_REF] Waters | New client puzzle outsourcing techniques for DoS resistance[END_REF], which form the majority of the existing computational puzzle schemes. Abadi et al. [START_REF] Abadi | Moderately hard, memory-bound functions[END_REF] first noticed the fact that CPU power varies a lot for different computers (such as PC, PDA, and Workstation), and introduced memory-bound computational puzzle schemes, where the computation is measured by the amount of memory look-ups needed to solve a puzzle. The schemes in [START_REF] Doshi | Efficient memory bound puzzles using pattern databases[END_REF][START_REF] Dwork | On memory-bound functions for fighting spam[END_REF] fall into this category. Regardless of the nature of different forms of computations involved in solving puzzles, the essence is the same, namely the prover needs to spend a certain amount of resources (either CPU cycles or memory look-ups) in finding a solution.

Computational puzzle schemes have been explored to realize a number of secure functionalities, such as timed encryption [START_REF] Rivest | Time-lock puzzles and timed-release crypto[END_REF], timed bit commitment [START_REF] Boneh | Timed commitments[END_REF], timed release of digital signatures [START_REF] Garay | Timed release of standard digital signatures[END_REF], uncheatable benchmarks [START_REF] Cai | Towards uncheatable benchmarks[END_REF], trustworthy website usage metering [START_REF] Franklin | Auditable metering with lightweight security[END_REF], and generating delays in lottery applications [START_REF] Goldschlag | Publicly verifiable lotteries: Applications of delaying functions[END_REF]. The rationale is that the process of spending resources in solving a puzzle automatically results in a time delay, which enables the verifier to control when the prover is able to access the functionalities. In other words, the functionalities are masked by puzzles, which need to be solved first before the prover can access the functionalities. Besides realizing various security functionalities, researchers have also applied computational puzzle schemes to mitigate a wide range of DoS attacks, such as fighting junk emails [START_REF] Dwork | Pricing via processing or combatting junk mail[END_REF][START_REF] Dwork | On memory-bound functions for fighting spam[END_REF], protecting authentication protocols [START_REF] Aura | DoS-resistant authentication with client puzzles[END_REF][START_REF] Lee | A public-key based authentication and key establishment protocol coupled with a client puzzle[END_REF][START_REF] Ning | Mitigating dos attacks against broadcast authentication in wireless sensor networks[END_REF], protecting IP networks [START_REF] Dean | Using client puzzles to protect TLS[END_REF][START_REF] Feng | The design and implementation of network puzzles[END_REF][START_REF] Fraser | Using client puzzles to mitigate distributed denial of service attacks in the tor anonymous routing environment[END_REF][START_REF] Juels | Client puzzles: A cryptographic countermeasure against connection depletion attacks[END_REF][START_REF] Wang | Defending against denial-of-service attacks with puzzle auctions[END_REF][START_REF] Wang | Mitigating bandwidth-exhaustion attacks using congestion puzzles[END_REF][START_REF] Wang | A multi-layer framework for puzzle-based denial-of-service defense[END_REF], protecting wireless networks [START_REF] Lei | Client puzzles based on quasi partial collisions against DoS attacks in UMTS[END_REF][START_REF] Martinovic | Wireless client puzzles in IEEE 802.11 networks: Security by wireless[END_REF], and preventing Sybil attacks [START_REF] Borisov | Computational puzzles as sybil defenses[END_REF]. With a computational puzzle scheme implemented, the server (playing the role of verifier) can mitigate an attack by asking every client (playing the role of prover) to solve a puzzle before allocating any required resource. The rationale is that, the number of "valid" requests from a malicious client will drop to some extent because the client has only limited resource to find puzzle solutions. In the literature, there has been some debate on whether computational puzzle schemes are really helpful to defect DoS attacks. Based on the collected results from ISPs in UK, Laurie and Clayton [START_REF] Laurie | Proof-of-Work Proves not to Work[END_REF] claimed that computational puzzle schemes are hardly effective in combating junk emails in practice. While, Liu and Camp [START_REF] Camp | Proof of Work can Work[END_REF] argued that computational puzzle schemes could be helpful if such schemes are used in combination with reputation systems.

Chen et al. [START_REF] Chen | Security notions and generic constructions for client puzzles[END_REF] proposed a security model for computational puzzle schemes tailored for defeating DoS attacks. In particular, their definition has focused on the unforgeability and difficulty of puzzles. However, their definition of puzzle difficulty implies neither the determinable difficulty property nor the parallel computation resistance property.

Definition of Computational Puzzle Schemes

In the design and analysis of computational puzzle schemes, we typically do not directly talk about how many CPU cycles or memory lookups are required for puzzle generation, puzzle solving, and puzzle solution verification. Instead, we often use the number of some generic operations (such as a hash or a multiplication in some group) as a puzzle complexity metric. For the notation purpose, we assume the generic operation to be denoted as Func.

A computational puzzle scheme consists of four (probabilistic) polynomialtime algorithms (Setup, PuzzleGen, PuzzleSol, PuzzleVer).

• Setup( , D): Run by the verifier, this algorithm takes a security parameter and an upper bound D on the puzzle hardness (see the definition below) as input, and outputs the public system parameter params and a private key mk. The public system parameter params should also include a specification about the metric function Func. This system parameter is implicitly part of the input to other algorithms, and we omit it in the descriptions for simplicity reasons. • PuzzleSol(puz): Run by a prover, this algorithm takes a puzzle puz as input and outputs a puzzle solution sol.

• PuzzleVer(mk, in f o, sol): Run by the verifier, this algorithm takes the private key mk, the puzzle information in f o, and the solution sol as input, and outputs 1 if sol is correct or 0 otherwise.

It is worth noting that a prover may not follow the PuzzleSol algorithm to find a solution, and in fact it can use any means to find the solution. But the properties determinable difficulty and parallel computation resistance, which will be defined in next subsection, imply that a prover need to perform d operations of Func in order to find a solution whether or not it follows the PuzzleSol algorithm. In addition, in some schemes, a long-term private key mk may be unnecessary or the additional information req may not be required as part of the input. In some application scenarios such as timed cryptography [START_REF] Boneh | Timed commitments[END_REF][START_REF] Garay | Timed release of standard digital signatures[END_REF][START_REF] Rivest | Time-lock puzzles and timed-release crypto[END_REF], a puzzle verification algorithm PuzzleVer is not explicitly required. Nonetheless, the above definition is meant to be general enough to cover most existing computational puzzle schemes.

Property Formulations for Computational Puzzles

Similar to other types of cryptographic schemes, soundness is a basic requirement for a computational puzzle scheme. In our case, soundness is defined to be that, given a puzzle with the hardness parameter d, the prover can find a correct solution by performing a total number of d operations of Func.

In the following, we formalize the properties determinable difficulty and parallel computation resistance. It is worth noting that, as in the case of most cryptographic formulations, we only consider a polynomial time adversary in our formulation and define that a property is achieved if an adversary has only a negligible advantage in the attack. The notion of negligibility is defined as follows.

(k) < 1 f (k) for all k ≥ N f .

The Determinable Difficulty Property

Informally, the determinable difficulty property implies that the verifier can precisely determine the required resource from the prover in solving a puzzle. Basically, this has two implications.

• One is that puzzles, generated by the verifier, are independent from each other in the sense that solving one puzzle does not help solve another puzzle. In other words, any puzzle generated by the verifier will require the prover to perform d Func operations.

• The other is that puzzles, generated by the verifier, are unpredictable in the sense that any puzzle generated by the verifier looks fresh so that the prover is unable to pre-compute the solution.

An attack succeeds if, after receiving a puzzle of hardness d, the adversary finds a solution by performing less than d Func operations. As a standard practice, attacks against the determinable difficulty property is simulated through the following three-phrase game between an adversary and a challenger, as depicted in Fig. 1. In the attack game, the adversary plays the role of a malicious prover and the challenger plays the role of the verifier.

In practice, the only information to a prover is the puzzles generated by the verifier and the public system parameter. It is straightforward to check that the adversary has all the same privileges as a malicious prover in practice. Formally, the determinable difficulty property is defined as follows.

Definition 2 Suppose that the adversary has performed d † Func operations in the Response phase of the attack game (defined in Fig. 1). A computational puzzle scheme achieves the determinable difficulty property if the probability d † < d * is negligible with respect to the security parameter . The probability is computed based on the (random) coin tosses of the challenger and the adversary.

Response phase:

The adversary can issue queries to the PuzzleGen oracle as in the challenge phase. At some point, the adversary terminates by outputting a correct solution sol * to the puzzle puz * .

Figure 1: The Attack Game

The Parallel Computation Resistance Property

As illustrated by the example schemes in Section 1, if a prover has access to a number of computers, it is able to solve a puzzle much faster than others by dividing the workload and letting the computers work in parallel. Informally, the parallel computation resistance property implies that a prover cannot accelerate the puzzle solving process by exploiting the parallelism. More concretely, the parallel computation resistance property requires that, after receiving a puzzle of hardness d from the verifier, a prover needs to sequentially perform d Func operations to find a solution. It means that the prover is unable to speed up the process by letting more than one computer work in parallel. This also implies that the best strategy for the prover is to use its fastest computer to solve the puzzle. We first introduce the notion of sequential computing time as follows.

x = {q 1 , q 2 , • • • , q x },
where, 1 ≤ x ≤ N and for all 1 ≤ i ≤ x -1, the q i+1 -th Func evaluation is started by a computer after the q i -th Func evaluation has ended (by any computer). The sequential computing time of the computation task is the maximum value of the sequence indexes, namely max x .

Formally, the parallel computation resistance property is defined as follows.

Definition 4 Suppose that the sequential computing time of the computation task of finding a correct solution sol * is d † in the Response phase of the attack game (defined in Fig. 1). A computational puzzle scheme achieves the parallel computation resistance property if the probability d † < d * is negligible with respect to the security parameter . The probability is computed based on the (random) coin tosses of the challenger and the adversary.

According to our definition, the parallel computation resistance property implies the determinable difficulty property. However, the reverse is clearly not true.

Parallel computation resistance is a very desirable property when we want to limit the disparity of computing powers of potential provers of a computational puzzle scheme. Moreover, in some cases, this property is crucial to mitigate attacks from malicious provers which control a cluster of Zombie computers. However, the determinable difficulty property is of independent interest in evaluating the security of computational puzzle schemes when the property parallel computation resistance is unnecessary. There are straightforward ways to convert a scheme, which achieves the parallel computation resistance property, into a scheme which achieves determinable difficulty and yet supports parallelism. For example, to generate a puzzle of hardness d, the verifier generates x independent sub-puzzles {puz 1 , • • • , puz x } such that each of the sub-puzzle has the hardness d x . By doing so, the total hardness is still d = d x • x, but a prover can solve the sub-puzzles in parallel.

Further Remarks

When designing computational puzzle schemes, a naturally desirable property is the computation disparity between the verifier and a prover: it should take very little resource for the verifier to generate a puzzle and verify a solution, while it should take the prover a certain amount of resources (determined by the verifier) to find a solution. Formally, this property can be evaluated by a ratio parameter v d , where v is the total number of Func operations required in generating a puzzle and verifying a solution, and d is the puzzle hardness.

With respect to puzzle hardness d, it is desirable that a computational puzzle scheme provides finer granularity support. Take the hash-based client puzzle scheme as an example 1 : a puzzle is in the form of (h, x 2 ), where h = H(x 1 ||x 2 ), H is a collision-resistant hash function, and x 1 has bit-length k; a solution is a k bit x 1 such that h = H(x 1 ||x 2 ). Clearly, the puzzle hardness can only be set exponentially with respect to k in this case. However, in the RSW computational puzzle scheme [START_REF] Rivest | Time-lock puzzles and timed-release crypto[END_REF] and the variant in Section 4, the puzzle hardness can be set linearly in the number of multiplications. Clearly, with the latter, the verifier can more flexibly set d according to the requirements of the underlying applications.

According to the definition of a computational puzzle scheme given in Section 2, the verifier may store some state information in f o for each puzzle. In this case, the scheme is said to be stateful. In practice, a stateful computational puzzle scheme may be considered to be inefficient because the verifier needs to store a piece of information for every unverified puzzle. Nevertheless, it is straightforward to turn a stateful scheme into a stateless one by sending both puz and in f o to the client, which should send sol and in f o back to the server for verification. In the case that in f o compromises the determinable difficulty property or the parallel computation resistance property, the verifier can protect in f o with a symmetric key encryption algorithm and a message authentication code algorithm. It is worth noting that a completely stateless computational puzzle scheme may be prone to replay attacks, in which an adversary replays a puzzle and its solution to the verifier. A more detailed discussion of this issue is beyond the scope of this paper, and should be addressed in the deployment of a computational puzzle scheme.

Analysis of a Variant of the RSW Scheme

In this section, we analyse a variant of the RSW computational puzzle scheme and show that it achieves the parallel computation resistance property. In the literature, Tritilanunt et al. [START_REF] Tritilanunt | Toward non-parallelizable client puzzles[END_REF] proposed a scheme towards achieving the parallel computation resistance property, but they did not give a rigorous security proof. It is unclear whether their scheme can be proven in our security model.

Description of the Scheme

The algorithms of the new scheme are defined as follows.

• Setup( , D): This algorithm selects two random large primes p, q and a cryptographic hash function H : {0, 1} * → Z * pq . The public parameter is ( , D, pq), and the master key is mk = (p, q). In addition, the metric function Func is multiplication in Z * pq . Given , the prime numbers p, q should be chosen in such a way that the required assumptions hold, as stated in next subsection.

• PuzzleGen(mk, d, req): This algorithm chooses r ∈ R Z pq and computes g = H(r||req), and outputs the puzzle puz = (g, d). The related puzzle information is in f o = (r, d, req).

• PuzzleSol(puz): This algorithm outputs sol = g 2 d mod pq.

• PuzzleVer(mk, in f o, sol): This algorithm returns 1 if sol ≡ (H(r||req)) 2 d mod φ(pq) (mod pq), and returns 0 otherwise.

In this scheme, g is computed as g = H(r||req), which can be regarded as a randomly chosen element of Z * pq if H is modeled as a random oracle (as shown in our security proof). By doing so, if needed, this element can be bound to situational information (such as the identity information of the prover) contained in req. Furthermore, the modification also facilitates the proof in our security model. Note that, in [START_REF] Rivest | Time-lock puzzles and timed-release crypto[END_REF], g is directly assumed to be a randomly chosen element of Z * pq . We further note that the hah function H can be replaced with any function that maps r||req to a random element in Z * pq .

The puzzle hardness d can be set linearly in the number of multiplications in the group Z * pq . It is clear that the workload of the verifier is constant regardless of the puzzle hardness. With respect to the verification complexity of the verifier, we omit the computation 2 d mod φ(pq) for two reasons. One is that it could be pre-computed and stored by the verifier. The other is that, in many cases, multiple puzzles might share the same hardness so that the computation only needs to be done once. As a consequence, it is straightforward to calculate that the average verification complexity for the verifier is approximately 3L 2 multiplications in Z * pq , where L is the bit-length of φ(pq). The computation disparity ratio parameter is 3L 2d .

Security Analysis

In [START_REF] Leander | On the equivalence of rsa and factoring regarding generic ring algorithms[END_REF], Leander and Rupp model Z pq as a generic ring structure in order to prove the equivalence of RSA and factoring. In this paper, we assume the same structure. For our purpose, we only describe a generic algorithm according to the multiplication operation in Z * pq in the following. It is worth noting that there are different ways of giving the definition (e.g. [START_REF] Nechaev | Complexity of a determinate algorithm for the discrete logarithm[END_REF][START_REF] Shoup | Lower bounds for discrete logarithms and related problems[END_REF]), and we follow that of Shoup [START_REF] Shoup | Lower bounds for discrete logarithms and related problems[END_REF].

Let g ∈ R Z * pq and σ be an encoding function of G = {g i | i ∈ N} on {0, 1} |pq| , where N is the set of integers and |pq| means the bit-length of pq. Suppose O is a multiplication oracle, which, for any r ∈ G, computes σ(r) as follows: if r has been calculated before, then the same value of σ(r) is returned; otherwise, it sets σ(r) to be a random value from {0, 1} |pq| \ S, where S is a set initialized to be {σ(g)}. A generic algorithm A is a probabilistic algorithm, which takes S and pq as input, and behaves as follows. At any time, A can send a query with the input (x, y, b), where x, y ∈ S and b ∈ {1, -1}, to O, and will receive σ(x • y b ) as the reply. After every query, the result is added to the set S.

Remark 1

The concept of generic algorithms, which traces back to Babai and Szemeredi [START_REF] Babai | On the complexity of matrix group problems i[END_REF], has been extensively used to model algebraic objects, such as groups and rings. It has made it possible to prove the security of many cryptographic protocols. Note that this methodology does possess an inherent limitation, namely a lot of attacks against cryptographic protocols have actually been subject to the successful exploit of the structures of the underlying algebraic objects. Nevertheless, a security proof can still provide a certain level of confidence in the security of the studied protocols.

To complete our proof, we further need the following computational assumption associated with the multiplication group Z * pq . This assumption is similar to the representation assumption in groups of prime order made by Brands [START_REF] Brands | Untraceable off-line cash in wallets with observers (extended abstract)[END_REF]. Definition 5 Let p, q be two prime numbers. The extended discrete logarithm assumption holds for Z * pq if the following event only occurs with a negligible probability with respect to a security parameter : Given (pq, g, g i (1 where g, g i (1 ≤ i ≤ V) are chosen uniformly at random from Z * pq and V is any polynomial in the security parameter, a polynomial-time adversary finds x 0,

≤ i ≤ V)),
x i (1 ≤ i ≤ V) ∈ N such that g x ≡ 1≤i≤V x i ∈N g x i i (mod pq).
In the following theorem, we prove that the variant scheme achieves the parallel computation resistance property (and certainly also the determinable difficulty property).

Theorem 1 If the adversary is modeled as a generic algorithm as defined above, then the variant scheme achieves the parallel computation resistance property based on the extended discrete logarithm assumption (given in Definition 5) in the random oracle model. Proof sketch. In order to prove the theorem, we need to show that the adversary's advantage in the attack game (defined in Fig. 1) is negligible according to Definition 4. Since the adversary is modeled as a generic algorithm, there are two types of oracle queries it may issue. One is the PuzzleGen query which can be issued to the challenger, and the other is multiplication oracle query which can be issued to the generic group oracle O.

We first consider a simple situation where the adversary does not issue any PuzzleGen query in the game. In this situation, the best strategy for the adversary to issue oracle queries to O in the Response phase is the following.

1. At the beginning, the adversary issues two queries with the inputs (σ(g * ), σ(g * ), 1) and (σ(g * ), σ(g * ), -1). Clearly, until it receives the replies, namely σ((g * ) 2 ) and σ(1), from the oracle O, it does not make sense for the adversary to send any other query.

2. After obtaining the replies, the set S becomes {σ(1), σ(g * ), σ((g * ) 2 )}.

Then the adversary issues queries with the following inputs (σ(g * ), σ((g * ) 2 ), 1), (σ((g * ) 2 ), σ((g * ) 2 ), 1), (σ(1), σ(g * ), -1), (σ(1), σ((g * ) 2 ), -1).

Clearly, until it receives the replies, namely σ((g * ) 3 ), σ((g * ) 4 ), σ((g * ) -1 ), and σ((g * ) -2 ) from the oracle O, it does not make sense for the adversary to send any other query. ), σ((g * ) 3 ), σ((g * ) 4 ), σ((g * ) -1 ), σ((g * ) -2 )}.

Then the adversary issues queries with (x, y, b), where x, y ∈ S, b ∈ {1, -1}, and σ(x • y b ) S. Clearly, until it receives the replies, it does not make sense for the adversary to send any other query.

4. The adversary continues the above process until it sends the response to the challenger.

Note that, we have enlarged the adversary's ability here, by allowing it to issue potentially an exponential number of queries to the oracle O. In fact, a polynomial-time adversary can only issue a polynomial number of queries. Nonetheless, if such an extensively empowered adversary cannot win the game, then a polynomial-time adversary can neither.

Suppose, at the end of the game, the adversary performs d * -1 steps as above (which also means the sequential computing time is d * -1), then the set S is a subset of

{σ((g * ) -2 d * -1 ), σ((g * ) -2 d * -1 +1 ), • • • , σ(1), • • • , σ((g * ) 2 d * -1 -1 ), σ((g * ) 2 d * -1 )}
In the game, the adversary could choose to submit a value σ((g * ) z ) from S or a value r from {0, 1} |pq| \ S as the response to the challenger.

• In the first case, if σ((g * ) z ) = σ((g * ) 2 d * ), then the adversary has actually found 2 d * -z such that (g * ) 2 d * -z ≡ 1 (mod pq). Based on the extended discrete logarithm assumption, the probability is negligible.

• In the second case, if (g * ) 2 d * has not been queried to the oracle, the probability r = σ((g * ) 2 d * ) is Q pq where Q is the total number of queries issued to the oracle O. Clearly, this probability is negligible. Otherwise, if (g * ) 2 d * has been queried to the oracle O, the probability r = σ((g * ) 2 d * ) is 0.

To summarize, in this simple situation, the adversary's advantage is negligible.

Next, we evaluate a more complex situation where the adversary is free to issue PuzzleGen oracle queries. This reflects the practical situation that a prover can ask the verifier for puzzles in a computational puzzle scheme. Our analysis will show that such a privilege does not help the adversary gain any benefit. Let σ(g 1 ), σ(g 2 ), • • • , σ(g V ) be the elements resulted from The set S 2 contains a polynomial number of encodings of the form σ( 1≤i≤V x i ∈N g x i i ), and the set S 3 contains a polynomial number of encodings of the form σ(A • B), where σ(A) ∈ S 1 and σ(B) ∈ S 2 .

In the game, the adversary could choose to submit a value σ((g * ) z ) from S 1 , S 2 , S 3 , or a value r from {0, 1} |pq| \ S as the response to the challenger.

• In the first case, from the analysis in the simple situation, the adversary's advantage is negligible.

• In the second and the third cases, if σ((g * ) z • 1≤i≤V x i ∈N g x i i ) = σ((g * ) 2 d * ), then the adversary has actually found 2 d * -z such that (g * ) 2 d * -z ≡ 1≤i≤V

x i ∈N g x i i (mod pq). Based on the extended discrete logarithm assumption, the probability is negligible.

• In the forth case, from the analysis in the simple situation, the adversary's advantage is negligible.

As a result, in this situation, the adversary's advantage is negligible. The theorem follows.

Conclusion

In this paper, we have revisited the concept of computational puzzle schemes, proposed a security model and presented formal definitions for two important properties, namely the determinable difficulty property and the parallel computation resistance property. We have proven that a variant of the RSW client puzzle scheme achieves both properties. To our knowledge, this is the first scheme which has been proven possessing the parallel computation resistance property. Compared with many other schemes such as the hash-based ones [START_REF] Back | Hashcash -amortizable publicly auditable cost functions[END_REF][START_REF] Juels | Client puzzles: A cryptographic countermeasure against connection depletion attacks[END_REF], the RSW scheme is computationally more expensive for the verifier, which needs to perform one exponentiation in verifying

•

  PuzzleGen(mk, d, req): Run by the verifier, this algorithm takes the private key mk, a hardness parameter d, and some additional information req as input, and outputs a puzzle puz and some relevant information in f o. The hardness parameter d is an integer which indicates the total number of Func operations required. The verifier sends puz to the prover, and keeps in f o for verifying the solution.

6 Page 6 Definition 1

 661 of 19 URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.The function P(k) : Z → R is said to be negligible with respect to k if, for every polynomial f (k), there exists an integer N f such that P

7

 7

1 . 2 .

 12 Setup phase: The challenger takes the security parameter and an upper bound D on the puzzle hardness as input, and runs the Setup algorithm to generate the master key mk and the public system parameter params. Challenge phase: The adversary issues queries to the PuzzleGen oracle with its chosen d and req, and obtains the reply puz = PuzzleGen(mk, d, req) from the challenger. At some point, the adversary sends d * and req * to the challenger, which sends puz * = PuzzleGen(mk, d * , req * ) back as a challenge.

11

 11

12 Page

 12

13

 13

3 .

 3 After obtaining the replies, the set S becomes {σ(1), σ(g * ), σ((g * ) 2

14

 14

S 1 =

 1 PuzzleGen oracle queries, where g 1 , g 2 , • • • , g V are random elements from Z * pq . Certainly, V is a polynomial in the security parameter. Suppose, at the end of the game, the adversary performs d * -1 steps in the Response phase, then the set S is a subset of S 1 S 2 S 3 . The set S 1 is the following. {σ((g * ) -2 d * -1 ), σ((g * ) -2 d * -1 +1 ), • • • , σ(1), • • • , σ((g * ) 2 d * -1 -1 ), σ((g * ) 2 d * -1 )}

  Let H be a one-way hash function. The verifier selects d random numbers r i (1 ≤ i ≤ d) as the puzzle, while the prover computes H(r i )(1 ≤ i ≤ d) as the puzzle solution. The verifier verifies a submitted puzzle solution h i (1 ≤ i ≤ d) through checking the equalities

	International Journal of Computer Mathematics
	verifies a submitted puzzle solution x through checking the equality H(x) = H(x ).
	r o F 2.
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

1 Page 1 of 19 URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.com

Definition 3

 3 Suppose that a prover has access to a number of computers which can work in parallel, and the Func function is implemented as a program in every computer. Let the each evaluation of the Func function be denoted as an event, marked by the starting time and ending time of the evaluation, and the set of events be denoted by g t (1 ≤ t ≤ N) in a computation task. Consider all sequences of the

	events of the following form
	Seq
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	8

  puzzle solution. Therefore, an interesting direction is to investigate how to improve the verifier's efficiency in pracice. Another interesting direction is to analyse the memory-bound computational puzzle schemes in our security model.

	F o r F o r F o r F o r
	P P P P
	e e e e
	e r e r e r e r
	R R R R
	e e e e
	v i e v i e v i e v i e
	w w w w
	O n l O n l O n l O n l
	y y y y
	15

a


Bartek Wierzba
email: abwierzba@agh.edu.pl

Marek Danielewski

Entropy Production during Inter Diffusion under Internal Stress

Keywords: diffusion, alloys, elasticity user supplied): entropy, entropy production, linear irreversible thermodynamics entropy, entropy production, linear irreversible thermodynamics, interdiffusion, mechano-chemistry, volume continuity equation PACS. 05.70.-a, 05.70.Ln

A theory of stress and elastic deformation during interdiffusion is presented. We show its consistency with the general Darken analysis and linear irreversible thermodynamics. Special consideration is given to the local entropy production. We derive the expression for the rate of entropy production during interdiffusion and present the practical computations of the internal energy and entropy densities and the entropy production. It is shown, that the entropy produced by the diffusion of mass is positive, 0 s A ≥ , and it does not depend on the frame of reference. The paper spans the gap between the Darken method, the linear irreversible thermodynamics and treatments by Larché, Cahn and Stephenson.

its equilibrium value in mechanical systems obeying Newton's equations [START_REF] Groot | Non-Equilibrium Thermodynamics[END_REF]. The entropy production fascinates many physics researches, moreover in the past ten years there has been renewed interest in thermodynamics of heat engines; many papers address issues of maximum power, maximum efficiency and minimum entropy production both from practical and theoretical point of view [2, 3,4,-5,6].

Providing a microscopic expression for the entropy production has been one of the grand aims of statistical mechanics, going back to the seminal work of Boltzmann. However, both the range of validity of the second law and of its proposed derivations have, from the very beginning, generated discussion and controversy [7].

The deformation during diffusion was investigated by Stephensohn [8]. He noted that when the components have different partial molar volumes or diffusion coefficients their simultaneous diffusion produce the imbalance in the volume transport which generate the elastic strain and internal stress. His theoretical framework shows the evolution of the stress fields during diffusion despite the implicit non-compressible conditions. Recently, the special consideration is given to the generation of internal stress and vacancy chemical potential [9]. The relaxation of such potential occurs via plastic deformation and vacancy formation and annihilation. Svoboda and Fischer distinguish between interstitial and substitutional components and vacancies, including the role of non-ideal vacancy source and sink activity. This concept considers the interaction of the diffusion process with the deformation process. The stress tensor and corresponding elastic strain energy represent the additional mechanical driving force for diffusion and generation/annihilation of vacancies.

The purpose of this paper is to develop the method that takes into account the both effects of stress and the component chemical potential terms in the diffusion potentials as well as relaxation of these terms, the resulting convective mass transport and the entropy production. Following

Stephenson we do not postulate the site conservation and consider multicomponent system. The method differs from the already published treatments [8], particularly in allowing the compressible elastic deformation and the effective use of the volume continuity equation [START_REF] Spalding | A Method for Computing Steady and Unsteady Flows Possessing Discontinuities of Density[END_REF][START_REF] Danielewski | [END_REF]. We show the general and exact mathematical description of the entropy production rate in the multicomponent, compressible solution which is valid for arbitrary initial conditions. The results are presented describing the entropy production in arbitrary binary diffusion couple. It is shown that the entropy production rate during diffusion is positive, does not depend on the frame of reference and depends on the initial conditions. Moreover, we show the deformation of the sample and the hydrostatic pressure field within the diffusion couple.

Model.

The Liouville transport theorem is a key mathematical tool used in this work to obtain the volume continuity equation and other conservation laws. It generalizes the Gauss-Ostrogradzki theorem and consequently can be applied for the vector fields (momentum, energy and entropy conservation in this work) [START_REF] Chorin | A Mathematical Introduction to Fluid Mechanics[END_REF].

We consider a solid solution at constant temperature, T const = , and in a closed system, 1 ,..., r m m const = , where m i represent a constant mass of the components.

Entropy -general description.

By the second law of thermodynamics the entropy, s, of the mixture contained by Ω ,

3 i i Ω = Ω ⊂

U

denote the volume occupied by the mixture at time t, is affected by the local entropy sources, s A , (due to the friction, diffusion, heat transport etc.) [START_REF] Groot | Non-Equilibrium Thermodynamics[END_REF][START_REF] Balluffi | Kinetics of Materials[END_REF]:

( ) We assume, that 3 R i Ω ⊆ represents any domain of the three-dimensional space and i υ is the velocity field defining the set i Ω . From the Liouville theorem the following identity holds ( )

d d d d div d d d i s s i i i i s s A J t t ρ ρ Ω Ω Ω Ω = Ω = - Ω ∑ ∫ ∫ ∫ (1)
d d div d . d i i i i i i i i i i i i i i s s s t t Ω Ω ∂   Ω = + Ω   ∂   ∑ ∑ ∑ ∫ ∫ ρ ρ ρ υ (2)
Combining Eqs. ( 1) and ( 2) and omitting the integrals the differential form of general equation for the entropy density follows:

( )

div div s s i i i i i i i s s J A t ρ ρ υ ∂ + + = ∂ ∑ ∑ (3)
2.2. Entropy when diffusion is not negligible. The thermodynamics postulates the following relation defining the entropy per mass unit, s [START_REF] Groot | Non-Equilibrium Thermodynamics[END_REF]:

( ) 1 2
, , , ,...,

m m m m r s s u N N N = Ω ,
where , ,

m m i u N Ω
denote the internal energy, the mixture specific volume per mass unit and the mass fraction (

m i i N ρ ρ =
), respectively. Multiplying the entropy by mixture density and splitting over the components results in partial Gibbs relation [START_REF] Danielewski | [END_REF]:

i i i i i i T s u ρ ρ ρ µ = - , (4)
where T is temperature and i µ denotes the diffusion potential, [START_REF] Groot | Non-Equilibrium Thermodynamics[END_REF][START_REF] Danielewski | [END_REF]. The differential of Eq. ( 4) is given by the relation:

ch m i i i p µ µ = + Ω
i i i i i i i i i i i i i i s u T Ts T s u t t t t t t t ρ ρ ρ µ ρ ρ ρ µ ρ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + + = + - - ∂ ∂ ∂ ∂ ∂ ∂ ∂ . (5)
Rearranging and combining with Gibbs equation, Eq. (4), we get:

i i i i i i i i s u T T s t t t t µ ρ ρ ρ ρ ∂ ∂ ∂ ∂ = - - ∂ ∂ ∂ ∂ , (6)
summing Eqs. ( 6) over all of the components and multiply by 1 ρ -, becomes:

m m m i i i i i i i i i s u T T N N N s t t t t µ ∂ ∂ ∂ ∂ = - - ∂ ∂ ∂ ∂ ∑ ∑ ∑ . ( 7
)
Using the Gibbs-Duhem relation, d d 0

m i i i s T N + =
∑ µ , we get:

m m i i i i i i s u T N N t t ∂ ∂ = ∂ ∂ ∑ ∑ (8)
which satisfies the second law of thermodynamics [START_REF] Groot | Non-Equilibrium Thermodynamics[END_REF].

F
Ω Ω Ω Ω = Ω = - Ω ∑ ∫ ∫ ∫ ρ ρ , (9)
where e J is the total energy flux through the mixture boundary. From the Liouville theorem the following identity follows:

( )

d d div d . d i i i i i i i i i i i i i i e e e t t Ω Ω ∂   Ω = + Ω   ∂   ∑ ∑ ∑ ∫ ∫ ρ ρ ρ υ (10)
Introducing Eq. ( 10) into ( 9) and omitting the integrals the differential form of energy conservation in the mixture follows:

( )

div div . e i i i i i i i e e J t ∂ + = - ∂ ∑ ∑ ρ ρ υ (11)
The component overall energy per mass unit, i e , includes the kinetic, internal and potential contributions:

2 1 2 = ext i i i i i i i e u V ρ ρ υ ρ ρ + + , (12)
where the internal energy, i u , follows from the Gibbs relation, (

where

m i i i i i i u N u u Ts µ = = + ∑ . ( 13
)
The total energy flux, e J , includes the flux due to the mechanical work performed on a mixture, p υ ⋅ I , and a heat flow, q J :

e q J p J υ = ⋅ + I , ( 14
)
where p is hydrostatic pressure, I represents the identity tensor and υ denotes the overall volume velocity,

m i i i i υ ρ υ = Ω ∑ .
Using Eqs. ( 11), ( 12) and ( 14) the balance equation for the overall energy becomes:

( ) ( ) ( ) 2 2 1 2 1 2 div div . I ext i i i i i ext q i i i i i i i i i i u V u V p J t ∂ + + + + + = - ⋅ + ∂ ∑ ∑ ρ υ ρ ρ ρ υ υ ρ υ ρ υ υ (15)
Assuming here the time independent external forcing ( 0) = ∑ , acting on the mass in Ω [1]:

ext V t ∂ ∂ = , the internal energy conservation law becomes: ( ) ( ) ( ) 2 2 1 1 2 2 div div div q i i i i i i i i i i i i i i u u p J t t ρ ρ υ ρ υ υ ρ υ υ ∂ ∂ + = - ⋅ + - - ∂ ∂ ∑ ∑ ∑ ∑ I (16)
d d d d d d d i i m i i i i i i i F t t ρυ ρ υ Ω Ω Ω Ω = Ω = Ω ∑ ∑ ∫ ∫ ∫ . ( 17
)
The rate of momentum of the mixture contained by Ω is affected by the force of elastic stress ( )

div m i i i p ρ Ω ∑ I and external force grad ext i i V ρ -∑ where ext V
is external potential (i.e.

gravitational field):

( )

div grad m ext i i i i i i i F F p V ρ ρ = = Ω - ∑ ∑ ∑ I , (18)
where

[ ) ( ) 3 3 3 : 0, , p L = ∞ × → I σ o is the elastic stress tensor and ( ) 3 3
,

L
is the space of linear mappings from 3 to 3 . In Eq. ( 18) we postulate that the ratio of the partial elastic stress tensor of the i-th component ( ) i σ to the overall elastic stress tensor ( ) σ equals the ratio of volume occupied by i-th component in one mole of the mixture to the overall volume of the mixture:

1 m i i i ρ -= Ω σ σ I . Thus we have ( 1 m m i i i ρ ρ Ω = Ω = ∑ ): m i i i ρ = Ω σ σ (19)
and consequently:

m i i i i i ρ = = Ω ∑ ∑ σ σ σ .
The Liouville theorem when the momentum is considered implies the following identity:

( )

d d div d . d i i i i i i i i i i i t t Ω Ω ∂   Ω = + Ω   ∂   ∫ ∫ ρ υ ρ υ ρ υ υ (20)
Combining Eqs. ( 17), ( 18) and (20) and omitting the integrals the differential form of partial equation of motion follows:

( ) ( ) div div grad m ext i i i i i i i i p V t ρ υ ρ υ υ ρ ρ ∂ + = Ω - ∂ I (21)
By summing Eqs. ( 21) for all components we get overall equation of motion in all media in which diffusion is non-negligible:

( )

div div grad ext i i i i i i i p V t ρ υ ρ υ υ ρ ∂ + = - ∂ ∑ ∑ I . ( 22
)
2.5. Pressure definition. We take into account the case when overall dilatation, t e , is a result of stress, e σ , and of the interdiffusion 1 (Darken velocity): The overall strain tensor, t ε , and the overall displacement vector u related by:

( )

1 grad grad 2 t T ε = + u u . ( 24
)
From Eq. ( 24) it follows that the overall dilatation is defined by:

div t e = u . ( 25
)
The elastic dilatation, e σ , equals [START_REF] Landau | Fluid Mechanics[END_REF]:

( )

3 1 2v e p E σ - = - . ( 26
)
The Darken strain tensor, D ε , and the Darken displacement vector u D are related by: ( ) ( )

1 grad grad 2 T D D D ε = + u u . ( 27
)
From Eq. ( 27) the Darken dilatation is defined by: div

D D e = u . ( 28
)
From Eq. ( 23) it follows:

t D e e e t t t σ ∂ ∂ ∂ = + ∂ ∂ ∂ . ( 29
)
Combining Eqs. (25), ( 26), ( 28) and (29):

( )

3 1 2 div div D v p E t υ υ - ∂ = - ∂ , ( 30
)
where

t υ = ∂ ∂ u and D D t υ = ∂ ∂ u .

Volume Continuity Equation.

The volume is not conserved, however the partial molar volumes obey the Euler relation:

1 m m m i i m i N ρ Ω = Ω = Ω ∑ . ( 31
)
The mixture volume, is the property transported by the velocity field of every mixture component contained by Ω :

d d d d d d d d d i m m i i i i t t t ρ ρ Ω Ω Ω Ω Ω = Ω Ω = Ω ∑ ∫ ∫ ∫ (32)
During an arbitrary transport process when volume is affected by the distribution of every mixture component and the stress field (the temperature is not considered here), from the Liouville theorem it follows:

( ) ( )

d d div d d m m m t t Ω Ω   ∂ Ω Ω Ω = + Ω Ω   ∂   ∫ ∫ ρ ρ ρ υ . ( 33
d div d . d m m t t Ω Ω   ∂ Ω = + Ω Ω   ∂   ∫ ρ ρ υ (34)
We consider here a multicomponent solution in a closed system, 1 ,..., r m m const = , where partial molar volumes and elastic properties do not depend on composition, ( )

1 ,..., i r f N N Ω ≠
. Thus, the volume can be only affected by the external stress:

d d d t x σ υ ∂ Ω Ω = ∫ . ( 35
)
Thus, introducing Eq. ( 35) into (34) the differential form of VCE follow:

( )

div div m m t σ ρ ρ υ υ ∂ Ω + Ω = ∂ . ( 36
)
Introducing Euler relation, Eq. ( 31) we can calculate the Darken velocity, D υ , of the mixture:

m D m d i i i i σ σ ρ υ υ υ ρ υ υ Ω = + + Ω = ∑ , ( 37
)
where d i υ is diffusion velocity defined by proper constitutive equation (section 2.8) and σ υ is a velocity generated by external elastic deformation.

Mass conservation law.

The mass of the mixture component, m i (t), contained in i Ω at the moment t equals:

( )

d i i i i m t ρ Ω = Ω ∫ (38)
The principle of conservation of mass states that the mass contained by i Ω is conserved2 . Thus,

d d 0 d i i i t ρ Ω Ω = ∫ . ( 39
)
By the Liouville theorem:

( )

d d div d d i i i i i i i i t t Ω Ω ∂   Ω = + Ω   ∂   ∫ ∫ ρ ρ ρ υ . ( 40
)
Combining Eqs. (39) and (40) and omitting the integrals results in local mass conservation laws for all components:

( ) ( )

div 0 i i i t ρ ρ υ ∂ + = ∂ . ( 41
grad grad d ch m i i i i i i B B p υ µ µ = - = - + Ω , ( 42
)
where i µ is the diffusion potential and B i denotes the mobility of the i-th component.

We consider mixture were the drift velocity is a result of the external elastic deformation ( ) σ υ and diffusion (the Darken drift, D υ ). Thus, the component velocities are defined by: where ( , ) ( , ) ( , )

d drift drift D i i t x t x t x σ υ υ υ υ υ υ = + = + . ( 43
)
2.9. Entropy production and entropy flux. Entropy production and entropy flux can be calculated by comparing Eqs. ( 3) and ( 8). Introducing the mass conservation law, Eq. ( 41) into overall entropy formula, Eq. ( 3), we get:

grad div s s i i i i i i i s s J A t ρ ρ υ ∂ + + = ∂ ∑ ∑ . ( 44
)
On the other side, upon introducing the internal energy conservation law, Eq. ( 16) and equation of motion, Eq. ( 21) into Eq. ( 8) we get:

( ) ( ) 1 1 1 grad div div q m i i i i i i i i i i i i T s u p J p t T T ρ ρ υ υ υ ρ ∂ = - - ⋅ + - Ω ∂ ∑ ∑ ∑ ∑ I I . ( 45
)
Rearranging the final entropy balance equation can be written in terms of fluxes and forces as:

div 1 1 . q i i i i i i i m q d d i i i i i i i i i i i i s J p s t T p J u T T T T   ∂ + + ∇ + =   ∂     Ω = ∇ + ∇ - ∇ + ∇     ∑ ∑ ∑ ∑ ∑ υ ρ ρ ρ µ υ ρ ρ υ υ υ (46)
Comparing Eqs. ( 44) and ( 46) the local entropy production rate and the entropy flux equal:

1 1 m s q d d i i i i i i i i i i i i p A J u T T T T ρ µ υ ρ ρ υ υ   Ω = ∇ + ∇ - ∇ + ∇     ∑ ∑ ∑ , ( 47
) q s J p J T υ + = . ( 48
)
Above relations allow to formulate the initial boundary-value problem of diffusion.

Formulation of the free-boundary value problem of interdiffusion under the stress field.

Data: ( )

1) ( ) 0, , 1
0, t ∈ ∞ duration of the process 6) ( ) 0, T ∈ ∞ the constant process temperature 7) ( ) 0, E ∈ ∞ the Young modulus 8) ( ) 0, v ∈ ∞ the Poisson number
The unknowns:

1) { } [ ) * * 0 : 0, i t t t ≤ ≤ Ω → ∞ × U ρ
( ) div , i d drift i i i t ∂ = - ∂ + ρ ρ υ ρ υ the mass conservation law for the i-th component. 2) d div d div d d x x t σ υ υ Ω Ω Ω = = ∫ ∫
, The volume continuity equation (VCE)

3) ( )

3 1 2 div div D v p E t υ υ - ∂ = - ∂ evolution of pressure in the mixture 4) grad : grad m i i i i i i i i u u p t ρ ρ ρ υ υ ∂ + = -Ω ∂ I the internal energy conservation of the i-th component,
where the overall internal energy of the mixture is Ts u µ =partial Gibbs-Duhem relation. The overall entropy of the mixture is calculated from:

i i i u u ρ ρ = ∑ .
i i i s s ρ ρ = ∑ . 6) grad grad m s d i i i i i i i i p A T T ρ µ υ ρ υ   Ω = -     ∑ ∑
, the entropy production rate in case of isothermal process (i.e. the temperature gradient is negligible).

Initial conditions:

1)

[ ) 0 : 0, i ρ Ω → ∞ the initial density distribution of the i-th component 2) 0 0 i i i c M ρ = ∑
the initial overall molar concentration of the mixture where n represents the unit vector normal to the boundary i.e., in the closed system analyzed here, the mass flow through the mixture boundaries does not occur.

2) grad 0 on

p p n n ∂ = = ∂ Ω ∂
the pressure derivative at the boundaries equals zero, i.e., there is no external forcing.

Results.

In this section we present the internal energy evolution and the entropy generated when diffusion is not negligible. The data used to simulate the mechano-chemical diffusion process in the arbitrary binary alloy are shown in Table 1. We consider an ideal solid solution ( ) Stress effect on interdiffusion. In Fig. 1 the impact of the pressure field on diffusion is presented.

i i a c = in
The solid line shows the concentration profile when the stress (pressure) field is considered The method allows calculating the entropy production during diffusion process, Fig. 3.

The entropy produced by the mass diffusion is positive, 0 s A ≥ , and it does not depend on the frame of reference, Eq. ( 47).

The initial conditions markedly affect the entropy produced during diffusion process. Figure 4 presents the evolutions of concentration for four different diffusion couples shown in Tables 1 and2.

Figure 5 shows the comparison of the drift velocity, pressure, internal energy and entropy in case of different initial composition. The different behavior of the evolution e.g. pressure profile can be observed.

Finally in Fig. 6 we show that the entropy produced by the mass diffusion in the case of different initial conditions is positive, 0 s A ≥ .

Summary.

We have shown the consistency of the generalized Darken method (bi-velocity method)

with the linear irreversible thermodynamics. We defined the material fixed frame of reference and effectively used the volume continuity equation. The resulting system of physical laws, initial and boundary conditions (i.e., the initial boundary value problem) allows using the bi-velocity method to model the complex mechano-chemical processes. We derived the expressions for the rate of entropy production during interdiffusion and presented the evolution of the internal energy, entropy and entropy production. It is shown, that the entropy produced by the diffusion of mass is positive, 0 s A ≥ , and it does not depend on the frame of reference.

The presented method spans the gap between the Darken method, linear irreversible thermodynamics and treatments by Larché, Cahn and Stephenson. We have shown its effective use to predict 1) the time dependence of the mixture concentration, 2) the internal pressure field generated during diffusion, 3) the energy and entropy of the mixture and 4) the entropy production term.

In this work we considered the binary solid solution. It ternary and higher systems the multiple diffusion paths are possible and often observed. The minimization of the entropy production can allow to solve this long lasting problem, i.e., can in a unique way define the evolution of the diffusion paths.

  , where s J is entropy flux, and i ρ ρ denote the overall and components mass densities in the mixture and the entropy density is a sum of all components

2 . 4 .

 24 Equation of motion. By Newton's law the rate of momentum change equals the overall force,

1

  It can be visualized as a spatially non-uniform Darken drift velocity, e.g., in a case of the stress free interdiffusion when

  ,..., the molar volume of the i-th component, i=1,2,…,r

  density of the i-th component, where i = 1,…,r.

  of the i-th component at the mixture boundary,

.

  It can be seen that stress field slows down the diffusion process. Spatial distribution of the: a) drift velocity, b) pressure field, c) internal energy and d) entropy in the binary alloy is shown in Fig. 2.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .

 123456 Figure 1. The influence of the internal stress field on diffusion in the binary alloy. 289x202mm (150 x 150 DPI)

2.3. Energy conservation law.

  According to the principle of the energy conservation, the overall energy of the mixture contained by Ω is affected by the energy flow into the volume Ω through its boundary ∂ Ω[START_REF] Groot | Non-Equilibrium Thermodynamics[END_REF]:

	d d t	e	d	d d t	i	i	i i e	d	i	div d e J
				o r
					P
						e
							e r
									R
									e
										v i e
										w
										O n l
										y

  Combining Eqs. (32) and (33) the Volume Continuity Equation (VCE) follows:

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	)

Table 1 .

 1 The data used to simulate the diffusion under the stress field in arbitrary binary alloy.

	Component Component diffusivity,	Molar mass,	Molar volume at p = 10 5 Pa
		D i [cm 2 s -1 ]	M i [g mol -1 ]	Ω i [cm 3 mol -1 ]
	1	1.0 •10 -9	100	10
	2	1.0 •10 -10	200	20

Table 2 .

 2 Initial concentrations in the diffusion couples shown in Figs.4-6.

	Philosophical Magazine & Philosophical Magazine Letters	Page 15 of 21
		Component concentration, at. %
	Diffusion couple	[ x ∈ -0.5, 0	)	[cm]	x ∈	[	] 0, 0.5	[cm]
	1	0					100
	2	10					100
	3	50					100
	4	75					100

http://mc.manuscriptcentral.com/pm-pml

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

To simplify the relations, in this work we do not consider the chemical and/or nuclear reactions in the mixture, i.e., we do neglect the local sources and sinks of mass.

Acknowledgments. This work has been supported by the Ministry of Higher Education and Science in Poland, project NR15-0121-10. The software CADiff2D available from authors.


Gboyega A Adebayo

Bede C Anusionwu

Optical Properties of -Se Revealed by ab initio Calculations

Keywords: Density Functional Theory, Electronic Structure, Infrared Spectra, Amorphous Semiconductors, First Principles Theory, ab initio simulation, Elemental Semiconductors

Detailed optical functions and optical properties of a α-monoclinic phase of Selenium (α-Se) are presented over a wide range of angular frequencies ω by first-principles calculations. An optical band gap of 2.05eV and a mobility gap of 2.17eV ± 0.05 are observed in the present calculations. The imaginary part of the dielectric function has been obtained from the frequency dependent dielectric matrix and the use of the Kramers-Kronig (KK) relation allows us to determine the real part of the dielectric function. Optical functions and optical properties revealed the occurrence of a long range order in Se. For the first time, infrared and absorption spectra of monoclinic-Se were studied using the Kramers-Kronig (KK) relation by ab initio method. And lastly, detailed studies of other dielectric related optical constants such as reflectivity, transmission, reflection and absorption coefficients were also presented.

Introduction

Selenium is a semiconductor possessing unusual properties which are not found in other single elemental semiconductors; aside from this fact, Se is also used in photovoltaic devices such as photoreceptors and photoconductors where the unusual behavior is an advantage in reducing the band gap of Se. Up till now, information on the optical properties of different phases of this element is still very scanty. Several studies on Se concern only the structures or quench behavior of this element [1] - [5], so detailed optical functions calculations are very useful ways of interpreting the anomaly in Se. Monoclinic Se and other phases of Selenium are important because of their vast applications in electronic and optoelectronic devices. Especially, amorphous Se is an excellent photoconductor used in active matrix flat-panel imagers, digital x-ray imaging and many more devices [6] - [13]. A proper understanding of the optical constants of Se is therefore of great importance in the design of modern-day optoelectronics based technologies.

That Se exhibits many photo-induced phenomena is well known [14] - [START_REF] Tauc | Amorphous and Liquid Semiconductor[END_REF] and although, there seems to be some attention focused on the experimental and theoretical studies of the optical constants of Se [START_REF] Innami | [END_REF], [19] there seems to be a need for more insight based on quantum mechanical calculations. This is because many of the observed effects in Se have been, practically, studied either over a small range of the optical spectrum or the majority of available calculations are based on vibrational density of state (DOS). This calls for a more accurate method to determine the optical properties of Se. First principle or ab initio calculations using the Density Functional Theory (DFT) approach coupled with the use of the pseudopotential provide the exact approach for obtaining accurate information on electronic band structures and optical properties over the entire electromagnetic spectrum.

In this article, we present the optical properties of monoclinic α-Se from the frequency-dependent dielectric function over a wide photon-energy range. The present calculations explores the frequency dependent dielectric matrix to determine the imaginary part of the dielectric function from which the Kramers-Kronig relation has been applied to obtain the real part. This allows us to obtain better optical properties of Se over the fundamental lattice vibration DOS approach. To our knowledge, no frequency-dependent dielectric matrix based ab initio calculation exists on monoclinic α-Se.

The present work is organized in the following sections: Section 2 gives the important theoretical account of the computational procedure used in our work, while Section 3 states the main results and discussion. Finally, in Section 4, concluding remarks are provided.

Computational Technique

The absorption coefficient α and the complex refractive index ñ are among the most important optical constants of a material. The extinction coefficient, k(ω) which is the imaginary part of ñ and α are related through the equation

α(ω) = ω (2) (ω) ck(ω) (1)
where "c" is the speed of light. We obtained the reflectivity from the equation

R(ω) = (n(ω) -1) 2 + k(ω) 2 n(ω) + 1) 2 + k(ω) 2 (2)
while the optical conductivity is given as

σ(ω) = ω (2) (ω) 4π (3)
From a knowledge of the refractive index, it is easy to determine the transmission and the reflection coefficients by employing

T coef = 4n 1 (ω)n 2 (ω) (n 1 (ω) + n 2 (ω)) 2 (4)
and

R coef = n 1 (ω) -n 2 (ω) n 1 (ω) + n 2 (ω) (5)
Here, the subscripts are the indices of refractive in air and Se.

Using plane wave basis sets, electron-ion interactions are described by the PAW pseudopotentials in the frame work of the Vienna Ab initio Simulation Package [20][21].

After optimizing the system, a 16 × 16 × 16 mesh with a Monkhorst Pack was used to implement a DFT method. In the determination of the partial occupancies of the wave functions, we employed the tetragonal method by setting the width of the smearing to 0.1eV. A quasi-Newton algorithm [21] was employed to achieve instantaneous ground state ionic relaxation. In the present calculations, an appreciable number of empty states in the conduction band are required to estimate the dielectric function and optical properties, for this reason, a two-step DFT calculation is implemented, we have used 50 bands from which convergence with size of the system used has been observed.

Band transitions calculations were performed by determining the joint density of state (JDOS) based on the calculations of DOS functions

DOS = (N (E i ) -N (E i-1 ))/ΔE ( 6
)
Here ΔE is the energy difference between two grid points, where

N (E i ) = Ei -∞ n(E)dE (7)
In terms of photon energy, JDOS is then written as

JDOS(hω) = 4 V 2 v,c δ(E c -Ev -hω) ( 8 )
V is the volume of a sample, E v and E c are energies of the valence and conduction bands respectively. Therefore,

JDOS(hω) = ∞ -∞ N v (E)N c (E + hω)dE (9)
The spectra-dependence of the imaginary part of the dielectric function was obtained through a summation over empty states

(2)

αβ (ω) = 4π 2 e 2 Ω lim q→0 1 q 2 c,v,k 2w k δ(E ck -E vk -ω) × < u ck+ αq |u vk >< u ck+ β q |u vk > * ( 10
)
The function u ck refers to the periodic part of the wave functions at the k-point k.

Using the KK relation, we determined the real part (1) (ω) of (ω) thus [20] (1)

αω (ω) = 1 + 2 π P ∞ 0 (2) αβ (ω )ω ω 2 -ω 2 + iη dω (11)
Here, ω is an integration variable, while P is the principal value of integration which is needed to avoid singularity at ω = ω and η is a complex shift.

Results and Discussion

The absorption coefficient is an important aspect in the technology of optical waveguide engineering as it gives an insight into absorbed photons by a material. In Fig. 1 we show our calculated absorption coefficient α(ω) for α-Se. In the infrared region, there is an occurence of a weak "symmetric peak" centred at around 0.25eV and another peak situated at 16eV with minima depressions on both sides; furthermore, a split third peak is noticed at around 21.25eV with a single minimum depression to the right. Lastly, there are several oscillational small peaks between 26eV and 50eV.

A major deduction from Fig. 1 is that α-Se has an optical band gap of 2.05eV, and a mobility gap of 2.17eV ± 0.05 which corresponds to the observed value in amorphous Selenium [22], [23] at ambient condition. The small peaks region could be a region of multiple reflection where the absorption coefficient is very small.

Generally, optical properties of semiconductors are exhibited at a photon energy below 4eV and from band theory the density of states (DOS) diagram is a key to understanding the optical properties of an amorphous semiconductor like Selenium.

Consequently, density of states and absorption coefficient curves can be used to explain band and mobility gaps in Selenium. A quick glance at Fig. 1 and Fig. 4 show clearly that the only observed peak satisfying the condition mentioned above has an energy gap that approximately equals 2.05eV. This value is relatively close to measured values of 2.53eV and 2.3eV by Knights and Davis [24], and Caywood

and Taynai [25] respectively. Also, a perusal of these figures shows localized states in the energy band. A major difference between crystalline and amorphous states is the existence of localized states called mobility gap in the band gaps. Mobility gap in semiconductors can be attributed to the existence of a high degree of short range spatial order in atomic structures. In the present calculations, the mobility gap of (2.17 In Fig. 3, we present the transmission and reflection coefficients (F ig.3a) alongside the reflectivity (F ig.3b) in α-Se.

In the low frequencies region, there are few noticeable peaks (F ig.3b) around 16eV, 21eV and 36eV. The region of high reflectivity corresponds to points in the lowest transmission and highest reflection coefficients. These graphs depict a long range order in Selenium at high frequencies as can also be seen in Fig. 4 where a plot of JDOS is presented. The unusual behaviour of Selenium as photoreceptor and photoconductor which tends to reduce the band gap in Selenium can be explained in terms of decreasing optical band gap energy which may be induced as a result of multiple optical reflections like the one observed in the tail of Fig. 2 A mobility gap of 2.17ev ± 0.05 is observed between -2eV and 2.5eV. There are several other mobility gaps, some of which have overlapping electron-hole state interactions. One can infer from these overlaps therefore, superpositions of molecular orbital states and bearing in mind that linear combinations of these orbitals will result in both anitbonding and bonding states.

Although, the DFT method employed in the present calculations is highly accurate since it relies on quantum ab initio calculations, however, only calculated energies up to 15eV of the electromagnetic spectrum can be considered accurate. The above statement is borne out by the fact that Selenium is a semiconductor and at high energies, due to artifacts in calculations, bands may not be calculated correctly while at energies exceeding 15eV a calculated refractive index may become smaller than unity.

Conclusion

We

Alpha-Se

Fig. 2

 2 Fig.2 shows the imaginary part of the dielectric function due to α-Se. The curve shows typical features of amorphous materials. We see clearly the four-region divisions which we believe correspond to transmittance values according to which α-Se could be transparent, weak absorptive, strong absorptive and reflective. The oscillatory tail could be attributed to region of multiple reflections. The Kramers-Kronig relation of eqn.11 was derived under general conditions and is therefore valid for semiconductors as well as for normal crystals, while refractive indices are derived from dielectric functions, any large change in the imaginary part of the dielectric function and refractive index are very often indications of some differences in the bonding and thus the structure of a material. Therefore, in contrast to crystalline structures, the abrupt change in Fig.2 can be attributed to the amorphous nature of Selenium.

  while formation of stable glasses of Se in chalcogenide are made possible by excellent infrared transparency, and a high refractive index among others which are inferred from Figures1 and 3a. Fig.4portrays the existence of several localized states in Selenium, the electron-hole interaction effect in the energy distribution of those states can not be underestimated at a large hω. This is most noticeable between hω ≈ 33eV and 42.5eV where there are overlaps of electron traps and hole traps. There are deep shallow traps both in the valence and conduction bands of Fig.4, worthy of mention are the electron traps at ≈ 23eV which extend into a well ≈ 0.4 in arbitrary units.

  have obtained the optical functions and optical properties of Se by applying quantum mechanical-based ab initio calculations which give essential features of the calculated quantities over an extended photon-energy range. An optical band gap of 2.05eV has been obtained in the present work, while the mobility gap from the density of state calculation shows a value of 2.17eV ± 0.05. Four regions corresponding to transmittance values in the extinction coefficient of α-Se, a transparent region, weak and strong absorption regions and a reflective region. Other optical constants which are related to the dielectric functions such as absorption coefficient, transmission and reflection coefficients as well as the refractive index of α-Se were also presented.

Figure 1 :

 1 Figure 1: Absorption coefficient in Alpha Selenium at high pressure

Figure 2 :

 2 Figure 2: Imaginary part of the dielectric function due to α-Se

Figure 3 :

 3 Figure 3: Optical coefficients and Reflectivity of Selenium at high pressure

Figure 4 :

 4 Figure 4: The joint density of state due to α-Se

Parmanent address: Department of Physics, UNAAB. Email: gadebayo@ictp.it

Acknowledgments

One of us, GAA acknowledges financial support from the Royal Society of London under the International Incoming Visitor Grant.


Volkmar Leute
email: leutev@uni-muenster.de

The influence of the ideal mixing entropy on concentration profiles and diffusion paths in ternary systems

Keywords: Interdiffusion, Ternary systems, Thermodynamics, Computer simulation

The diffusion profiles and the reaction paths in ternary solid solutions are determined by both thermodynamics and kinetics. The matrix of the diffusion coefficient can be described as the product of the Hessian matrix for the thermodynamic influences and the Onsager matrix for kinetic influences.

In this paper the interest is focused on the influence of the ideal part of the Hessian matrix, i.e. the ideal mixing entropy on interdiffusion. The ideal diffusion profiles are calculated by a computer simulation and they are compared with experimental results from the literature. These comparisons reveal that in most cases the qualitative shape of the diffusion profiles and of the reaction paths can be considered as caused by the ideal mixing entropy. Surprisingly the shape of the diffusion profiles turns out to depend on the component that was chosen as the so-called solvent of the ternary mixture. This means that the ideal reaction paths do not show the triangular symmetry expected for an ideal ternary system. Especially reaction paths between starting positions showing the same concentration of one of the three components do not run along straight lines.

Introduction

By definition ternary systems are systems in which the composition can be described by two independent variables and the diffusion by a 2x2 diffusion matrix. Diffusion paths k(l) in ternary systems of type A k B l C m with k + l + m = 1 are calculated from two mole fraction profiles e.g. k(z) and l(z) by eliminating the space coordinate z for diffusion in one dimension. These mole fraction profiles are usually measured by an electron microprobe or they can be determined from computer simulations realized by solving the diffusion equations with a finite difference algorithm. [1,[START_REF] Kirkaldy | Diffusion in the Condensed State[END_REF].

The matrix D of the diffusion coefficients can be written as the product of the Onsager matrix L and of the Hessian matrix M of the Gibbs free energy [3,4]. In matrix notation one obtains:

D = 1 C (L • M) ( 1 )
C is here a scalar, namely the mean molar concentration of the components. The Onsager matrix considers the mobility part of the diffusion matrix, whereas the Hessian matrix describes the thermodynamic contribution to the diffusion process. Generally it is payed much too little attention to the contribution of the Hessian matrix and in the literature in some cases the Hessian matrix is even not correctly formulated. Therefore we will discuss in this paper the influence of the Hessian matrix in a thermodynamically ideal system on the diffusion profiles of ternary systems. For this purpose the Onsager matrix will be taken as a scalar and all thermodynamic interaction parameters as zero.

Theoretical considerations 2.1 The Onsager matrix L

In a ternary system A k B l C m with two independent mole fractions the flux of each component i ∈ {A, B, C} is principally a linear function of the driving forces X i [START_REF] De Groot | Grundlagen der Thermodynamik Irreversibler Prozesse[END_REF]:

J A = l AA X A + l AB X B + l AC X C J B = l BA X A + l BB X B + l BC X C (2) J C = l CA X A + l CB X B + l CC X C ,
but, because of i J i = 0, these fluxes are not independent. To get a description with independent fluxes and forces we have to reduce this equation system by taking one component -being called the solvent [6] -as depending on the others. Choosing C as the solvent we obtain [6,[START_REF] De Groot | Grundlagen der Thermodynamik Irreversibler Prozesse[END_REF]:

J A = L C AA (X A -X C ) + L C AB (X B -X C ) ( 3 ) J B = L C BA (X A -X C ) + L C BB (X B -X C )
Depending on the component that was chosen as solvent the matrices describing these equation systems are:

L C AA L C AB L C BA L C BB L A BB L A BC L A CB L A CC L B CC L B CA L B AC L B AA ( 4
)
The upper indices indicate, which component was chosen as the solvent. The term solvent is only a conventional designation for that component whose composition was chosen as depending on those of the other two components. It is used for the whole system covering all compositions including even such compositions, where this 'solvent' is the minor component.

The Onsager coefficients l ij of the equation system (2) can be described as functions of the new Onsager coefficients L ij in the matrices (4) [START_REF] De Groot | Grundlagen der Thermodynamik Irreversibler Prozesse[END_REF]. As the reduced 2x2 matrices belong to equation systems with independent fluxes and independent forces, the coefficients L ij can take any values, provided the restrictions L ii ≥ 0, L jj ≥ 0, L ii L jj ≥ 1/4(L ij + L ji ) [START_REF] De Groot | Grundlagen der Thermodynamik Irreversibler Prozesse[END_REF] and the Onsager relations L ij = L ji [8] are considered.

As our principal object in this paper will be to show the influence of thermodynamics on diffusion processes, we will simplify the Onsager matrices (4) as scalar, i.e. we will take the off diagonal elements as zero and all diagonal elements as constant and equal. Thus, such an Onsager matrix L will only affect the length of the diffusion profiles, but not the shape of the diffusion path.

The ideal mixing term of the Hessian matrix M

To calculate the elements of the Hessian matrix for a ternary system of type A k B l C m we need the mean molar Gibbs energy g as a function of two independent mole fractions, e.g. k and l. The Gibbs energy

g = kμ A + lμ B + mμ C (5)
can be split in the standard term g 0 , the ideal mixing term g id and the excess term g E :

g = g 0 + g id + g E → M = M 0 + M id + M E . ( 6
)
With component C as solvent we have to take the mole fractions k and l of the components A and B as the independent variables. Considering that k + l + m = 1 it follows from eq.( 5) that

∂g ∂k = μ A -μ C = μ A-C ∂g ∂l = μ B -μ C = μ B-C (7)
These are the chemical potentials of the difference building units (A-C) and (B-C) describing the exchange of a solvent atom C for one of the independent components A or B. The negative gradients of the chemical potentials of these building units are the independent forces according to eq.( 3). The partial differentials of these chemical potentials with respect to the independent variables k or l are the second derivatives of the mean molar Gibbs energy, i.e. the elements of the Hessian matrix with C as solvent:

∂μ A-C ∂k = ∂ 2 g ∂k 2 = M C kk , ∂μ A-C ∂l = ∂ 2 g ∂k∂l = M C kl , ( 8
)
∂μ B-C ∂l = ∂ 2 g ∂l 2 = M C ll , ∂μ B-C ∂k = ∂ 2 g ∂l∂k = M C lk . ( 9
)
As the chemical standard potentials of the pure components are independent of composition, the standard term M 0 vanishes, but nevertheless also in thermodynamically ideal systems the three components are still energetically different. The ideal mixing term g id is restricted to the ideal mixing entropy, because the ideal mixing enthalpy is zero per definition.

g id = -T s id = RT • ln k k • l l • m m with m = 1 -k -l. ( 10
)
Thus, for k and l as the independent variables, the elements of the ideal mixing term M id of the Hessian matrix become

M C,id kk = RT 1 k + 1 m , M C,id ll = RT 1 l + 1 m , M C,id kl = M C,id lk = RT • 1 m .
(11) It is astonishing that in the literature relevant to this subject [3] the elements of the ideal mixing term of the Hessian matrix are given as

M C,id kk = RT 1 k , M C,id ll = RT 1 l , M C,id kl = M C,id lk = 0. ( 12
)
The reason for this error is probably that, in differentiating the chemical potentials μij = ∂(μ iμ C )/∂x j with i, j ∈ {A, B}, it was not considered that the chemical potential μ C = μ 0 C + RT ln m of the dependent component C depends also on the mole fraction k and l via m = 1kl.

The appearing of the term 1/m in eq.( 11) has severe consequences. This term effects that even in ideal mixtures a thermodynamic coupling of fluxes occurs because of M C,id kl = M C,id lk = 0. Moreover, also the diagonal elements of the Hessian matrix depend on 1/m.

If -as in developing eq.( 11) -the mole fractions k and l of the components A and B were chosen as the independent variables, the component C belonging to the mole fraction m = 1kl is called -as already mentioned above -the solvent of this ternary system. The superscript C in eq.( 11) means that for this case C was chosen as the solvent. The corresponding cases for the solvents A and B read as follows:

M A,id ll = RT 1 l + 1 k , M A,id mm = RT 1 m + 1 k , M A,id lm = M X,id ml = RT • 1 k . ( 13
)
M B,id mm = RT 1 m + 1 l , M B,id kk = RT 1 k + 1 l , M B,id mk = M Y,id km = RT • 1 l .
(14) In order that calculated mole fraction profiles or diffusion paths for ideal systems were independent of the choice of the solvent, both the determinant and the trace of the ideal Hessian matrix would have to be invariant with respect to the choice of the solvent S [6]. The determinants d S of the matrices with the elements of eqs. [START_REF] Ho | Interdiffusion fluxes and transport coefficients under multiple gradients in selected ternary systems[END_REF], (13) or [START_REF] Dayananda | Diffusion in solids: recent developments[END_REF] are indeed invariant:

d A = d B = d C = 1 k • l • m , ( 15
)
but their traces are not invariant. Thus the simulation of mole fraction profiles and diffusion paths in ideal ternary systems must yield different results for different choices of the solvent. By the way, for the wrong matrices corresponding to eq.( 12) neither the determinants nor the traces are invariant with respect to the solvent choice, and therefore mole fraction profiles calculated on this basis would also have to depend on the choice of the solvent.

The usual expectation that diffusion paths in ternary ideal systems would run along straight lines is generally not true and the statement that in an ideal solution all elements of the diffusion matrix reduce to the purely kinetic form D ij = RT Lij/x j [9] does not hold either. Instead we have to consider that even in ideal systems the Hessian matrix is not diagonal and that their off-diagonal elements, M S,id ij with (i = j), must always depend on composition (eq.11,13,14).

Whereas the function g id (eq.10) shows the threefold symmetry of the Gibbs triangle, the matrix of its second derivatives (eq.11) does not so. Thus, the whole field of possible reaction paths in the thermodynamically ideal system does not yield the expected triangular symmetry, but only a mirror symmetry related to the line x = y for the case that x and y are chosen as the independent mole fractions. Moreover, the simulation of reaction paths for thermodynamically ideal ternary triangular systems shows that no reaction path yields a straight line except the path for the reaction between the component chosen as solvent and the equimolar mixture of the other two components.

Comparison of calculated diffusion profiles

with experiments from the literature

Cu-Ag-Au

Fig. 1a shows the phase diagram with the miscibility gap at 1000 K for the system Cu k Ag l Au m according to Ziebold and Ogilvie [10]. This figure contains a series of diffusion paths. For comparison we have calculated these diffusion paths for the ideal system X k Y l Z m as shown in Fig. 1c,d. Fig. 1c belongs to Z as solvent and Fig. 1d to X as solvent. One can clearly see that the shapes of these calculated diffusion paths depend on the solvent choice.

Comparing the calculated diffusion paths starting at pure X or Y with the experimental ones starting at Cu or Ag, one can see that the diffusion paths of Fig. 1c correspond qualitatively much better with the experimental paths in Fig. 1a than those of Fig. 1d. For the two diffusion paths with the starting points (l = 0, m = 0.55) and (k = 0, m = 0.45) or (l = 0, m = 0.55) and (k = 0, m = 0.7) one cannot decide so easily, because these calculated diffusion paths for the two solvent cases differ only quantitatively but not qualitatively. Nevertheless, it looks rather clear that in the system Cu-Ag-Au the component Au acts as solvent. This result is consistent with the experience that in ternary systems with a spinodal miscibility gap near the edge XY [4] the component on the corner Z opposite to the miscibility gap acts most probably as solvent. Though the given system deviates from the thermodynamically ideal behavior so much that it forms even at 1000 K a spinodal miscibility gap, one can understand the qualitative shape of the diffusion paths just from the ideal mixing entropy without considering an influence of the Onsager matrix or any interaction parameters for real systems.

The diffusion paths starting with one of the pure components show the well known S-shaped course, i.e. they cross the straight line joining the two starting points of the path only one time and their Au profiles (cf.Fig. 1b) Figure 1: (a) Experimental diffusion paths and the miscibility gap at 1000 K in the system Cu-Ag-Au (After [10]). (b) The calculated mole fraction profiles of Au with the same starting compositions as the diffusion paths in (a), but for a thermodynamically ideally behaving system Cu k Ag l Au m with Au as solvent. (c) Calculated diffusion paths for a thermodynamically ideally behaving system X-Y-Z for Z and (d) for X as solvent.

don't show any region of uphill diffusion.

The two paths in Fig. 1c between (l = 0, m = 0.55) and (k = 0, m = 0.45) or between (l = 0, m = 0.55) and (k = 0, m = 0.7), however, cross the straight line between the two starting points two times. Such a behaviour was until now seldom experimentally observed [START_REF] Ho | Interdiffusion fluxes and transport coefficients under multiple gradients in selected ternary systems[END_REF], but the computer simulations of diffusion paths for an ideal system show that this is a usual behavior for diffusion couples with about equal starting content of the solvent component in both partners of the couple. Fig. 1b shows the calculated m-profiles for the diffusion paths of Fig. 1c. In the profiles for the diffusion path between (l = 0, m = 0.55) and (k = 0, m = 0.7) (Fig. 1b) the Au flux runs from right to left including two regions of uphill diffusion (the path regions with negative slope). In experiments the tiny minima at the ends of the diffusion path can easily be overlooked, because they correspond in the Au profile to the two extended regions with very slowly changing Au concentrations outside the two shallow minima of the profile. This could be the reason for the missing minimum on the left side of the corresponding experimental reaction path in Fig. 1a.

In the diffusion path between (l = 0, m = 0.55) and (k = 0.55, m = 0.45) and in the m profile we have an equivalent situation, but now the Au flux runs from left to right. Interdiffusion processes in the system Ni-Co-Fe were experimentally intensively investigated by Vignes and Sabatier [START_REF] Vignes | [END_REF]. Fig. 2a shows the Ni and Co profiles for one of the many diffusion paths they had measured. By Kirkaldy and Young [13] this system was taken as a thermodynamically ideal solution with Fe as solvent and it was 'emphasized that the diffusion interactions are of entirely kinetic origin'. As we have already shown in [4] this is wrong, because already the influence of the ideal mixing entropy suffices to explain the qualitative behaviour of the experimentally determined diffusion paths. Thus, one cannot say that the redistribution of Ni is 'due primarily to the diffusional cross-effect represented by D 12 ' [13]. The redistribution is rather a thermodynamic effect, caused by the ideal mixing entropy as shown in Fig. 2b. Moreover one must not set D 21 = 0, because the off-diagonal elements of the Hessian matrix are equal and do not vanish even in thermodynamically ideal systems (M id 12 = M id 21 = 0, eq. [START_REF] Ho | Interdiffusion fluxes and transport coefficients under multiple gradients in selected ternary systems[END_REF]13,[START_REF] Dayananda | Diffusion in solids: recent developments[END_REF]).

Ni-Co-Fe

Thus, the diffusion behaviour in the system Ni-Co-Fe is predominantly determined by thermodynamics, namely by the ideal part of the Hessian matrix (Eq.11). For example, for the diffusion couple shown in Fig. 2, the computer simulations of the mole fraction profiles for a thermodynamically ideal system -even with the scalar Onsager matrix -yield the experimentally observed redistribution of Ni. The assumption in the literature that in this system Fe acts as solvent was confirmed by our computer simulations. Fig. 3a shows the mole fraction profiles measured by Dayananda in the system Cu-Ni-Zn at about 1000 K. The Ni profile is characterized by a light uphill diffusion. If the profiles of the three components are calculated by computer simulation for a thermodynamically ideal system, the region of clear uphill diffusion for the Ni profile can only be reproduced if the component Zn is chosen as solvent. For Ni as solvent one gets only a very tiny maximum in the Ni profile, and for Cu as solvent there is not at all any uphill diffusion for Ni.

Cu-Ni-Zn

These results reveal that the shape of the calculated reaction paths depends on which component was chosen as solvent and that by comparison of the calculated and experimentally determined profiles one can decide which component acts really as solvent.

Fe-Al-Ni

Some diffusion paths within the ordered β 2 -phase of the system Fe-Al-Ni were investigated by Moyer and Dayananda [START_REF] Moyer | [END_REF]. The left plots in Fig. 4 show the experimentally determined mole fraction profiles of these diffusion paths. The right plots show the mole fraction profiles of all three components Fe, Al and Ni calculated using a scalar Onsager matrix and the assumption of a thermodynamically ideal solution with Al chosen as solvent. Figs. 4a,b belong to a diffusion path with the same Al content and Fig. 4c,d to one with the same Fe content in both starting compositions. For such situations in ideally behaving systems it is usually expected that the diffusion paths yield straight lines between the starting points, i.e. the Al profile in Fig. 4b and the Fe profile in 4d would then have to be straight lines. But both, the experimentally determined and the calculated profiles show distinctly a redistribution of Al or Fe by the diffusion process. Al in Fig. 4a,b diffuses from the right to the left side, Fe in Fig. 4c,d, however from the left to the right side. In Fig. 4d the calculated redistribution is much weaker than the experimentally measured one. Nevertheless, the computer simulations of the mole fraction profiles yield -in spite of the restrictive assumptionsqualitatively the same results as the experiments of Moyer and Dayananda [START_REF] Moyer | [END_REF].

Discussion

The calculations have shown that, if the full ideal part of the Hessian matrix (eq.11) is considered, the shape of a diffusion path depends on the choice of the solvent. Until now, we do not know how one could predict from physical arguments which of the three components acts really as solvent. However, by comparing the calculated mole fraction profiles or diffusion paths with the experimental data we can determine, which component has to be chosen as solvent. Kirkaldy has already written -with regard to the discussion of diffusion coefficients in ternary systems -that 'the choice of solvent is by no means a trivial matter ' [6]. If in an ideal ternary system all three components could be taken as identical, the thermodynamic properties must not depend on the arrangement of the three components, i.e. the collectivity of the diffusion paths should obey the triangular symmetry of the phase triangle independent of the solvent choice. But, generally the three components are not identical. They still differ energetically by their chemical standard potentials μ 0 i , and therefore it is possible that an exchange of two pure components -or equivalently a change of the independent variables -will influence the shape of the reaction paths. Such an influence of the chemical nature of the components on the behaviour of ideal mixtures can also be detected for quasiternary systems of type

(M k N 1-k )(X l Y 1-l ), in which all four quasibinary subsystems (M k N 1-k )X, (M k N 1-k )Y, M(X l Y 1-l )
and N(X l Y 1-l ) behave ideally. In these cases the essential diffusion paths also don't run along straight lines (cf. Fig. 5b). Really, one could assume that in this ideal square system the shape of the diffusion paths would have to obey a mirror symmetry, i.e. an exchange of M for N or of X for Y should have no influence on the shape of the reaction paths. But Fig. 5b reveals that reflecting the diffusion path at l = 0.5 or at k = 0.5 yields new diffusion paths.

The Hessian matrix of such an ideal quasiternary system reads :

M id 11 = RT 1 k(1 -k) , M id 12 = M id 21 = -Δ R G 0 , M id 22 = RT

Summary

The comparisons of the experimental mole fraction profiles for extended solid solutions in ternary systems with the profiles determined by computer simulations for the same starting compositions show that most probably the qualitative behaviour of many experimental profiles and diffusion paths is essentially determined by the ideal mixing entropy of the solid solution.

Especially we could show by the computer simulations that the usually expected quasibinary sections, i.e. straight lines for the reaction paths, in thermodynamically nearly ideal mixtures do generally not occur. The single exception is the diffusion path between the component being chosen as solvent and the equimolar mixture of the other two components of the ternary system [4].

Thus, I think that, for the interpretation of diffusion profiles and reaction paths in ternary systems, it is essential to consider the influence of the ideal mixing entropy. The simulation of diffusion profiles and reaction paths for thermodynamically ideal conditions can yield valuable information on the qualitative behaviour of the diffusion process. In cases, where it is known which component acts as solvent -as for example in systems with a spinodal miscibility gap detected at lower temperatures -the reaction paths for nearly ideal systems can in many cases qualitatively be predicted by the computer simulation, even if the diffusion coefficients are not known.

Figure 2 :

 2 Figure 2: (a) Experimentally determined mole fraction profiles for Co and Ni in the system Ni-Co-Fe (after [12]); (b) Calculated profiles with the same starting compositions as in (a), but for the thermodynamically ideally behaving system with Fe as solvent.

Figure 3 :

 3 Figure 3: (a) Mole fraction profiles for Cu, Ni and Zn in the system Cu-Ni-Zn (after Dayananda [14]); (b) The calculated profiles for the thermodynamically ideally behaving system with the same starting positions as in (a) and with Zn as solvent.

Figure 4 :

 4 Figure 4: (a) and (c): Mole fraction profiles for Fe, Al and Ni in the system Fe-Al-Ni (after [15]); (b) and (d): The calculated profiles with the same starting positions as in (a) and (c) for the thermodynamically ideally behaving system with Al as solvent.

Figure 5 :

 5 Figure 5: (a) Calculated k and l profiles for the reaction between M(X 0.5 Y 0.5 ) and N(X 0.5 Y 0.5 ) at 1000 K in a thermodynamically ideal system (M k N 1-k )(X l Y 1-l ) with a standard Gibbs energy of the reaction MX + NY → NX + MY of Δ r G 0 = +20 kJ/mol. (b) The diffusion path calculated from the profiles in (a).

l(1l)(16) In contrast to this formulation the off-diagonal elements of the ideal matrix are generally taken as zero[1], and Δ R G 0 is then treated as part of the offdiagonal elements of the excess term. But, as Δ R G 0 = ν i μ 0 i , this term does not vanish in ideal systems, and as it only depends on the standard chemical potentials of the components of the reaction (R), MX + NY → NX + MY[16,17], it can alternatively be identified with the off-diagonal elements of the ideal part of the Hessian matrix. That means that, even in thermodynamically ideal systems, the shape of the diffusion path must change, if the chemical nature of the components changes. The uphill diffusion, shown by the l-profile in Fig.5a, is caused by the thermodynamic coupling because of the off-diagonal elements of the Hessian matrix.Similarly the role of the solvent in ideal triangular systems could depend on the values of the standard chemical potentials of the components, and if this is the case, the choice of the independent composition variables or of the component acting as solvent is not at our free disposal. Moreover, if one describes a reaction path l(k) of a ternary triangular system in a rectangular Cartesian coordinate system with the composition of the solvent at its origin, the mirror symmetry related to the diagonal k = l for an ideal system is reasonable, but one is not automatically obliged to look also for a triangular symmetry.