Björn A Persson
email: bjorn.persson@teokem.lu.se

Mikael Lund

Jan Forsman

Dereck Chatterton

Torbjörn Åkesson

Molecular Evidence of Stereo-Specific Lactoferrin Dimers in Solution $

Keywords:

Gathering experimental evidence suggests that bovine as well as human lactoferrin self-associate in aqueous solution. Still, a molecular level explanation is unavailable. Using force field based molecular modeling of the proteinprotein interaction free energy we demonstrate (1) that lactoferrin forms highly stereo-specific dimers at neutral pH and (2) that the self-association is driven by a high charge complementarity across the contact surface of the proteins. Our theoretical predictions of dimer formation are verified by electrophoretic mobility and N-terminal sequence analysis on bovine lactoferrin.

Main

This combined theoretical and experimental study of the pair interaction of bovine lactoferrin (LF) shows that the protein self-assembles into dimers. Free energy calculations reveal a well defined, stereo-specific complex that is stabilized by a region of complementary amino acid residues. Historically, LF was first isolated from bovine milk but the glycoprotein is found in all mammals. It is secreted in various exocrine fluids such as bile, saliva, pancreatic fluid and tears [START_REF] Adlerova | Lactoferrin: a review[END_REF] . Inflammatory reactions and viral infections are often followed by increased plasma concentrations of LF and the protein is considered part of the acute-phase immune response. More specifically, LF has antibacterial, antiviral and anti-inflammatory activity, underlining that the protein is an important part of the immune system [START_REF] Brock | The physiology of lactoferrin[END_REF]. LF is also interesting from a technical perspective, for instance, cationic protein-stabilized emulsions have been prepared of the protein and LF is known to be a fouling component of contact lenses [START_REF] Sarkar | Colloidal stability and interactions of milk-protein-stabilized emulsions in an artificial saliva[END_REF][START_REF] Meagher | Interaction between adsorbed lactoferrin layers measured directly with the atomic force microscope[END_REF]. LF has also received attention as a nutritional additive in milk formulas [START_REF] Tomita | Twenty-five years of research on bovine lactoferrin applications[END_REF]. The protein consists of two major lobes and approximately 700 amino acids. The macromolecule has the shape of a dumbbell and is well described by a bi-axial ellipsoid with half-axis of 47 Å and 26 Å. Several experimental studies imply that the protein self-assembles: Dimers of bovine LF have been confused with IgG2, that has a molecular weight twice of the LF monomer and sedimentation equilibrium has revealed a concentrationdependent average molecular weight [START_REF] Schanbacher | The similarity of bovine lactoferrin dimer to igg2[END_REF][START_REF] Castellino | Structural studies on bovine lactoferrin[END_REF]. Further, LF is found to elute with molecular weights corresponding to monomers, dimers and trimers, according to chromatographic analysis of bovine milk [START_REF] Harmon | Changes in lactoferrin, immunoglobulin g, bovine serum albumin and alphalactalbumin during acute experimental and natural coliform mastitis in cows[END_REF][START_REF] Wang | Identification of lactoferrin complexes in bovine mammary secretions during mammary gland involution[END_REF]. Scattering studies have also revealed large aggregates of human LF and that these dissolve at high ionic strength [START_REF] Babina | Effect of nucleotides on the oligomeric state of human lactoferrin[END_REF]. These observations all suggest that LF is able to self-assemble into larger structures, although no molecular level explanation has been presented. Using Metropolis Monte Carlo simulations we investigate the angularly averaged free energy of interaction between two LF molecules in an aqueous salt solution using the Faunus open source software package [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF][START_REF] Lund | Faunus: An object oriented framework for molecular simulation[END_REF]. This software package has been applied with success for computation of the virial coefficient of lysozyme and binding constants of lysozyme and apo-α-lactalbumin. [START_REF] Lund | A mesoscopic model for protein-protein interactions in solution[END_REF][START_REF] Persson | Association and electrostatic steering of alphalactalbumin-lysozyme heterodimers[END_REF] The simulation procedure is described in detail in the Appendix as well as elsewhere and we shall give here only a brief overview [START_REF] Lund | A mesoscopic model for protein-protein interactions in solution[END_REF][START_REF] Lund | On the charge regulation of proteins[END_REF]: Each protein is treated as a rigid body consisting of spherical amino acid beads that may be either charged or neutral, depending on the type and protonation state; PDB entry 1BLF was used to generate the coarse holo-Lactoferrin structure. A rigid model is reasonable since the structure of the holo-form is considered stiff and is resistant to proteolytic degradation. [START_REF] Baker | A structural framework for understanding the multifunctional character of lactoferrin[END_REF][START_REF] Adlerova | Lactoferrin: a review[END_REF] The solvent is treated as a dielectric continuum and particles interact via a combined Lennard-Jones (LJ) and salt screened Debye-Hckel potential, where the total energy for each configuration, U , is given by the sum of the amino acid pair interactions.

U = N i<j q i q j 4π 0 r r ij e -κr ij + 4 σ ij r ij 12 - σ ij r ij 6 (1)
r ij is the separation between amino acids, q is the charge of an amino acid and κ is the inverse Debye length, related to the ionic strength. 0 is the permittivity of vaccum and r the dielectric constant of water, taken as 80 in this work. σ ij is the mean diameter of the amino acids and is the LJ interaction parameter (0.05 kT ). k is Boltzmanns constant and T is the temperature. During the MC simulation, proteins are translated as well as rotated to explore configurations in the canonical ensemble at 298 K. The free energy of interaction (or potential of mean force, PMF) is obtained by sampling the protein-protein radial distribution function as a function of the mass center separation, R. In contrast to typical docking methods, the present approach correctly takes into account thermal motion and hence the entropy of the protein-protein association process. Figure 1A shows the calculated lactoferrin-lactoferrin free energy of interaction at pH 5, 7 and 9 and at an ionic strength of 5 mM (κ = 0.023 Å-1 ). At these pH values the LF net charge is +23e, +14e and +4e. Titration simulations estimate pI to 9.4, slightly above the experimental range, 8-9 [START_REF] Belegrinou | ph-dependent immobilization of proteins on surfaces functionalized by plasma-enhanced chemical vapor depostion of poly(acrylic acid)-and poly(ethylene oxide)-like films[END_REF][START_REF] Wahlgren | The adsorption from solutions of beta-lactoglobulin mixed with lactoferrin or lysozyme onot silica and methylated silica surfaces[END_REF][START_REF] Superti | Involvement of bovine lactoferrin metal saturation, sialic acid and protein fragments in the inhibition of rotavirus infection[END_REF]. At pH 5 the interaction is purely repulsive due to a large net charge of the proteins. How- ever, as pH is increased, a deep narrow minimum appears at approximately 47 Å. The attraction broadens at pH 9, close to isoelectric conditions, but the minimum is, interestingly, very distinct. We now demonstrate that the narrow free energy minimum is due to a few highly specific residue-residue interactions that lock the entire complex into a unique constellation. Firstly, by inspecting the MC generated configurations (see animation in Appendix) it is clear that when the proteins are in contact, only few orientations are observed. Secondly, we identify a close to orthogonal set of vectors, spanned from the proteins mass center to the mass centers of E385 and E66 and correlate the E385 angles with the E66 dihedrals (inset of Figure 2). This allows us to estimate the orientational free energy and, as evident from Figure 2, the two proteins are aligned in a strongly stereo-specific manner, dominated by a narrow set of protein-protein orientations. The choice of E385 was motivated by the fact that it is parallel to the line connecting the mass centers at the free energy minimum. The total interaction between the two proteins

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

is influenced by both electrostatic and van der Waals interactions. Due to the non-zero molecular net charge, any attractive electrostatic contribution must be due to charge patches on the protein surface. To qualitatively investigate the effect of this patchiness, we identify a few amino acids with high inter-protein charge-complementarity. These are selected by inspecting the protein-protein contact surface of different configurations at the free energy minimum. Two such candidates are glutamate E535 and lysine K498. In Figure 1B we present the PMF between mutated forms of LF as well as for the wildtype at pH 8, where the charge of the wildtype is +11e. Reverting the negative charge of E535 has dramatic effects on the potential of mean force: this single point mutation leads to a complete loss of the attractive free energy minimum. The mutation increases the net charge of the protein by two units and an increased repulsion can be expected. Conversely, the protein net charge is maintained by swapping charges of K498 with E535both of which are in the contact region. This also significantly diminishes the free energy minimum, underlining that the contact surface is highly tailored to facilitate a dimeric complex, assembled via complementary electrostatic interactions. Coordinates of a typical dimer configuration is supplied in the online material. Lastly, we provide experimental data to back our theoretical predictions. LF was separated using 12% Tris-Glycine SDS-PAGE under non-reducing conditions followed by staining with Coomassie blue R-250. Results indicate a primary band at 80 kDa, consistent with the molecular weight of the protein in its monomeric form. Two distinct bands, consistent with the molecular weight of dimeric and trimeric LF are also clearly visible at 160 kDa and 240 kDa (Figure 3). To confirm identity of these bands, electroblotting was performed using PVDF membranes (Invitrogen). The 80 kDa and 160 kDa electroblotted bands were excised and amino acid sequence analyses were performed on an Applied Biosystems PROCISE HT protein sequencer with on-line identification of phenylthiohydantoin-derivatives using standard programmes. Lactoferrin was identified in both bands.

While experimental data and the present simulations suggest the formation of LF dimer, it is desirable to further investigate special orientation of the aggregates by experimental scattering techniques. It would also be of importance to evaluate if the model is able to predict trimers. We intend to continue with such studies and to further scrutinize the driving force in the simulated system. Initial examination of more detailed models shows consistent results with this work.

A C C E P T E D M A N U S C R I P T

Acknowledgment

We would like to express thanks to the Royal Physiographic Society in Lund, Organizing Molecular Matter (Lineaus Center of Excellence), and Esben S. Sørensen, University of Århus for the sequence analysis.

Appendix A.

To coarse grain the protein we assume that each amino acid has a density of 0.9 g/ml and use the molecular weight to compute the radius. The sphere is centered at the mass center in the atomistic amino acid structure. It has been shown that this provides a good topological description of proteins and reproduces atomistic results [START_REF] Lund | A mesoscopic model for protein-protein interactions in solution[END_REF]. To evaluate the average charge of each amino acid we perform Monte Carlo simulations with single proteins. In such calculations, each amino acid is allowed to titrate in a Markov chain, according to a Boltzmann weight that is modified by the pKa-value and pH of the solution. This provides us with an estimate of the charge distribution of the proteins that are being used as input when we evaluate the pairinteraction. A complete description of simulation details and parameters can be found in the reference [START_REF] Lund | A mesoscopic model for protein-protein interactions in solution[END_REF] and [START_REF] Lund | On the charge regulation of proteins[END_REF]. The radii of some of the different amino acids, along with pKa-values used to generate average charges, are provided in Table A.1. Lactoferrin is here modeled in its holo-form. This is achived by placing a unit charge at the site of each ferric ion in the original coordinates, approximating a ferric ion synergistically bound with a carbonate. The magnitude of the Lennard-Jones interaction parameter is obtained

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

from a Hamaker argument. Since the densities of proteins are fairly constant one may estimate the contribution of the interactions across a water medium on an amino acid level. The magnitude in this work is chosen so it reproduces that of a previously published model [START_REF] Lund | A mesoscopic model for protein-protein interactions in solution[END_REF]. Since the range of interaction depends on the size of the amino acid, large residues will effectively attract each other more than smaller which is reasonable. This does not capture any sort of specific hydrophobic interactions. What is being modeled is rather an averaged dispersion-hydrophobic interaction and the geometrical influences of the experimental structure. In the simulations we restrict the sampling to distances shorter than 200 Å and the radial distribution function is sampled with a bin width of 1 Å. Each PMF is generated with approximately 10 7 configurations, which was found sufficient to produce static results with low noise. The animation provided displays a representative trajectory of the Monte Carlo simulation at pH 8 and clearly shows how the two proteins binds in a specific orientation. Red colored particles are negatively charged amino acids and blue colored are positively charged.

Figure 1 :

 1 Figure 1: Free energy of interaction (PMF) between lactoferrin molecules as a function of the mass center separation, R, at (A) different pH values for the wildtype, (B) different mutations and wildtype at pH 8. The ionic strength is 5 mM in all cases. Note, one kT equals 2.48 kJ/mol at 298 K. The maximum error is about 0.2 kT.

Figure 2 :

 2 Figure 2: Angular free energy (kT) of the connection angle (red vectors) and dihedral (blue vectors) at pH 8 and at 5 mM ionic strength. Sampling is performed within a 60 mass center separation. The inset shows the two proteins and the reference frames. The scale is relative to a free orientation.

Figure 3 :

 3 Figure 3: SDS PAGE analysis of LF. Proteins were separated using 12% Tris-Glycine SDS PAGE under non-reducing conditions and pH 8.3 following by staining with Coomassie blue R-250. Lane 1: Molecular weight marker; Lane 2: LF. N-terminal sequence analysis of the electroblotted bands is indicated.

Figure A. 4 :

 4 Figure A.4: Lactoferrin animation, click to play.

Table A . 1 :

 A1 Particle properties.

	Particle Radius ( Å) pKa
	Ala	3.1	-
	Arg	4.0	12.0
	Asp	3.6	4.0
	Gly	2.9	-
	Glu	3.8	4.4
	His	3.9	6.3
	Lys	3.7	10.4

A C C E P T E D

M A N U S C R I P T


Andreas Löpker
email: lopker@eurandom.tue.nl

David Perry
email: dperry@stat.haifa.ac.il

The idle period of the finite G/M/1 queue with an interpretation in risk theory

Keywords: Finite G/M/1, Finite M/G/1, Workload, Idle period, Sample path analysis, Level crossing, Risk process, Deficit at ruin Mathematics Subject Classification (2000) Primary 60K25, 91B30, Secondary 90B22, 68M20

HAL is

Fig. 1 The workload process V t and the associated risk reserve process R t . The idle period of the finite G/M/1 queue coincides with the deficit at ruin of the risk process with restricted accessibility," "finite-buffer queues," "uniformly bounded virtual waiting time," "limited queueing waiting time," "finite dam," etc. [7, 13, 15, 16, 19, 25-29, 31, 36]. The current paper can be seen as a continuation of [START_REF] Adan | The G/M/1 queue revisited[END_REF], with the following extra feature. We investigate a G/M/1 queue with restricted accessibility in the sense that the workload (virtual waiting time) V t is bounded by 1. If V t plus the service time of an arriving customer exceeds 1, only 1 -V t of the service requirement is accepted.

This paper focusses on the study of I , the idle period of the finite queue, i.e. the duration of the period between the time when the queue becomes empty and the next customer arrives. The random variable I can also be interpreted as the deficit at ruin of a modified risk reserve process R t in the compound Poisson case with a constant barrier strategy (see Fig. 1). When the risk reserve process reaches level 1, dividends are paid out with constant rate equal to 1, so that R t stays constant until the next claim occurs. Risk models with similar dividend strategies have been intensively studied in the insurance risk literature-often with a focus on optimality of the chosen barrier. We refer to [START_REF] Albrecher | Risk theory with a nonlinear dividend barrier[END_REF][START_REF] Albrecher | On the distribution of dividend payments and the discounted penalty function in a risk model with linear dividend barrier[END_REF][START_REF] Avram | Ruin probabilities and deficit for the renewal risk model with phase-type interarrival times[END_REF][START_REF] Gerber | On the probability of ruin in the presence of a linear dividend barrier[END_REF][START_REF] Lin | The compound Poisson risk model with a threshold dividend strategy[END_REF], where risk models with different dividend strategies are treated. Models with constant dividend barrier have been investigated in [START_REF] Li | On a class of renewal risk models with a constant dividend barrier[END_REF] and [START_REF] Lin | The classical risk model with a constant dividend barrier: Analysis of the Gerber-Shiu discounted penalty function[END_REF]. In [START_REF] Lin | The classical risk model with a constant dividend barrier: Analysis of the Gerber-Shiu discounted penalty function[END_REF] an integro-differential equation for the Gerber-Shiu discounted penalty function is found. Using this result the moments of the deficit at ruin, which is identical to our idle period I , are determined. In our study we focus on the law of I and find expressions for the Laplace-Stieltjes transform (LST) and distribution, the latter in the case with set-up time. We introduce two different methods to find these quantities. In [START_REF] Li | On a class of renewal risk models with a constant dividend barrier[END_REF] the model [START_REF] Lin | The classical risk model with a constant dividend barrier: Analysis of the Gerber-Shiu discounted penalty function[END_REF] was extended to hypoexponential claim inter-arrival times. For a general introduction to risk models, see [START_REF] Asmussen | Ruin Probabilities[END_REF] or [START_REF] Rolski | Stochastic Processes for Insurance and Finance[END_REF].

After introducing notation in Sect. 2, we give a short overview of the G/M/1 queueing model and related useful results in Sect. 3. An auxiliary result for the conditional distribution of idle period, given the workload process reaches level 1, is derived in Sect. 4. In the following sections we present two new methods to derive relations for the LST and probability distribution function of I ; both methods are based on a sample path analysis. The first method, demonstrated in Sect. 5, uses the idea of collecting subsequent overshoots over level 1 of the original G/M/1 workload process to form the service times of a new queue that can again be identified as a G/M/1 queue. Using this construction, the LST of the idle period of the finite queue without set-up time is derived. For the special case of an M/M/1 queue, the result is then checked in Sect. 6. The second method is presented in Sect. 7. It is based on the observation that the idle period of the original queue can be seen as the overshoot of the workload process in a dual M/G/1 queue. We use this observation and regeneration theory to construct a modified M/G/1-type process and derive a formula for the distribution of the idle period in the case with set-up time.

Level-crossing and sample-path based arguments similar to the ones used in the present paper can be found also in [START_REF] Brill | Single-server queues with delay-dependent arrival streams[END_REF][START_REF] Kaspi | Dam processes with state dependent batch sizes and intermittent production processes with state dependent rates[END_REF].

Preliminaries and notation

We have already mentioned the workload process V t , the set-up time a and the idle period I of the finite G/M/1 queue. Whenever it is important to distinguish the successive idle periods, we write I 1 , I 2 , . . . for the first, second, and the subsequent idle periods; otherwise we use a generic random variable I . Note that I 1 , I 2 , . . . form an i.i.d. sequence of random variables.

Let S 1 , S 2 , . . . denote the inter-arrival times of the customers and let F S be the distribution function of S 1 , with 1/μ = E(S 1 ). Let Z 1 , Z 2 , . . . denote the exponential service times, having mean E(Z 1 ) = 1/λ.

We let ρ = λ/μ, so that 1/ρ is the traffic intensity of the G/M/1 queue. Since we are concerned with the finite queue, we can ignore all stability issues and investigate both ρ > 1, when the standard G/M/1 queue is stable, and ρ ≤ 1, when it is unstable.

If not otherwise stated, we denote the LST E(e -sX ) of a random variable X by φ X (s). The convolution of two probability distribution functions F and G or two functions f and g is denoted by F * G and f * g. F j * and f j * denote the j -fold convolution of F and f respectively.

The standard G/M/1 queue

It is instructive to first review the standard G/M/1 queue and introduce some known relevant results for the standard case with no set-up time, i.e. a = 0. Let V t denote the workload process of this queue and let I denote the idle period. The LST of I is then given by

φ I (s) = λ • z -φ S (s) s -λ(1 -z) , ( 1
)
where z is the smallest positive root of z = φ S (λ(1z)) (see [START_REF] Prabhu | Stochastic Storage Processes[END_REF], p. 35, and [START_REF] Adan | The G/M/1 queue revisited[END_REF]) and φ S is the LST of S. An inversion is possible for ρ ≥ 1, when Lagrange's theorem yields (see [START_REF] Takács | Introduction to the Theory of Queues[END_REF])

z = ∞ j =1 λ j -1 j ! ∞ 0 x j -1 e -λx dF j * S (x).
In the case that ρ < 1 we have z = 1, so that (1) reduces to

φ I (s) = λ • 1 -φ S (s) s . ( 2
)
The distribution function of I is then given by 4) in [START_REF] Gerber | On the probability and severity of ruin[END_REF] for the risk process context). Note that (1φ S (s))/μs is the LST of the limit of the residual lifetime in a renewal process with epochs having distribution F S and that φ I (s) is a transform of a defective probability distribution function; in particular, P(

F I (x) = λ x 0 1 -F S (u) du (cf. (
I < ∞) = φ I (0+) = ρ.
In the sequel the so-called dual M/G/1 queue is obtained by interchanging the inter-arrival and service times, so that Z 1 , Z 2 , . . . denote the inter-arrival times and the variables S 1 , S 2 , . . . become the successive service times of the dual M/G/1 queue. Note that when ρ > 1, the G/M/1 queue is stable, while the dual M/G/1 queue is not. On the other hand, if ρ < 1 then the M/G/1 queue is stable, but not the G/M/1 queue. In the latter case, the workload process W t of the dual M/G/1 queue has a stationary distribution F with LST

φ F (s) = 1 -ρ 1 -φ I (s) , (3)
which is the transform version of the Pollaczek-Khintchine formula.

Let π be the probability that V t up-crosses level 1 during a busy period and let η be the probability that starting in 1, the process hits 0 before it returns to 1. It has been shown that if ρ < 1 then

π = 1 -ρ F (1) (4)
and

η = 1 - F * F S (1) F (1) = f (1) λ F (1) , ( 5
)
where f (x) is the density of F (x) for x > 0 and the second step of (5) follows from the Pollaczek-Khintchine formula. These formulas can be found in [START_REF] Adan | The G/M/1 queue revisited[END_REF][START_REF] Cohen | Extreme value distribution for the M/G/1 and the G/M/1 queueing systems[END_REF][START_REF] Cohen | The Single Server Queue[END_REF][START_REF] Perry | A controlled M/G/1 workload process with an application to perishable inventory systems[END_REF][START_REF] Ross | Hitting time in an M/G/1 queue[END_REF][START_REF] Takács | Application of Ballot theorems in the theory of queues[END_REF]. Note that in the notation of [START_REF] Ross | Hitting time in an M/G/1 queue[END_REF] the r.h.s. of ( 4) is given by 1q while the r.h.s. of ( 5) is denoted by 1p. In the notation of [START_REF] Adan | The G/M/1 queue revisited[END_REF], θ 1 (0, 1) represents the π in (4) and P(T - 1 < T + x ) is η in (5) (also note that in [START_REF] Adan | The G/M/1 queue revisited[END_REF] the symbols μ and λ have a reversed meaning). For ρ > 1 the probability π , the probability that the cycle maximum is larger than 1, is given in Theorem 4 of [START_REF] Adan | The G/M/1 queue revisited[END_REF]. Further references and results about the standard G/M/1 queue can be found in [START_REF] Adan | The G/M/1 queue revisited[END_REF][START_REF] Asmussen | Applied Probability and Queues[END_REF][START_REF] Brill | Level Crossing Methods in Stochastic Models[END_REF][START_REF] Cohen | The Single Server Queue[END_REF][START_REF] Kleinrock | Queueing Systems[END_REF].

Conditional idle period

We derive a formula for the conditional distribution of the idle period of the standard G/M/1 queue, given the event that the workload process exceeds level 1 during a busy period; or equivalently: we present an expression for the distribution of the deficit at ruin of the risk reserve process R t , given that some dividends were paid out.

Let V max be the maximum of V t during the first busy cycle and let H be the conditional distribution function of the idle period I , given the event { V max ≥ 1}.

Theorem 1

The conditional distribution H is equal to the residual lifetime distribution at time t = 1 of a renewal process with renewal times having the same distribution as the idle period I .

Proof By tracing Fig. 2 for a typical sample path of V t , let T * and T 1 denote the last up-and down-crossing times of level 1 before the idle period starts. T 1 is the endpoint of an excess period of V t over level 1. The time X 1 from T 1 to the next arrival has the same distribution as that of the G/M/1 idle period, since the interval T 1 -T * can be seen as the busy period of a G/M/1 queue (indicated by a gray area). At time T 1 + X 1 another busy period of a G/M/1 queue starts; it ends at time T 2 . Again, the distribution of X 2 is the same as that of an idle period and we see that this property also holds for X 3 , X 4 , . . . . Thus, the sequence X 1 , X 2 , . . . is i.i.d. and forms a renewal process.

Let I = X 1 + X 2 + • • • + X κ -1, where κ = inf{k|X 1 + • • • + X k > 1};
I can be seen as its residual lifetime at time t = 1 of the renewal process.

Expressions for the distribution of the residual lifetime of a renewal process can be found in [START_REF] Belzunce | A note on stochastic comparisons of excess lifetimes of renewal processes[END_REF] and [START_REF] Cohen | The Single Server Queue[END_REF]. We note that H is a solution of the renewal equation ( [START_REF] Asmussen | Applied Probability and Queues[END_REF], p. 143)

H (x) = F I (1 + x) -F I (1) + F I * H (x). ( 6
)
Fig. 3 Construction used in the proof of Theorem 2

The idle period without set-up time

In this section we let the set-up time a = 0. We are interested in the LST of the idle period I of the finite G/M/1 queue.

Theorem 2

The LST of I is given by

φ I (s) = φ L (s) π + (1 -π)φ L (s) , ( 7
)
with

φ L (s) = z 0 -z 0 -1 + s λη φ D (s), ρ ≥ 1 1 -s λη φ D (s), ρ < 1 ( 8
)
where z 0 is the smallest positive root of z 0 = φ L (λη(1z 0 )) and

φ D (s) = πφ H (s) 1 -φ I (s) + πφ H (s) . ( 9
)
Here φ H (s) = ∞ 0 e -su dH (u), where H is the conditional distribution of the idle period I , given the event { V max ≥ 1}.

Before we prove Theorem 2, we note that the LST φ I is given in [START_REF] Adan | The G/M/1 queue revisited[END_REF] and that φ H can in principle be calculated from Theorem 1. For ρ ≥ 1, (8) is implicit in the sense that φ L (s) is given in terms of the root z 0 , which itself can be determined only when φ L is known.

Proof The dashed line in Fig. 3(a) shows the standard G/M/1 workload process and the solid line shows the finite G/M/1 queue. Let K 1 be the number of overshoots U 1,1 , U 1,2 , . . . , U 1,K 1 of level 1 during the first busy period R 1 of the finite queue. Note that π = P(K 1 > 0) and η = P(K 1 = 1|K 1 > 0). We construct a third queue as follows. We collect the overshoots U 1,1 , U 1,2 , . . . , U 1,K 1 during R 1 and let them be processed during the idle period I 1 by a new server. As soon as a new customer arrives, this new server stops working and collects the overshoots U 2,1 , U 2,2 , . . . , U 2,K 2 during R 2 of the finite queue. The server then continues working during the second idle period, and so on. Once all the overshoots are processed, the server becomes idle until a busy period with overshoots ends. Finally we remove the inactive periods R 1 , R 2 , . . . , so that we end up with a situation as shown in Fig. 3(b). Let L k , A k and D k denote the inter-arrival times, service times and the idle periods of this new queue. Note that all three variables can be represented as geometric sums. Indeed, we have

L 1 = M i=1 I i , ( 10
)
where M is a geometrically distributed random variable with support M ≥ 1 and P(M = 1) = π , representing the number of idle periods of the finite queue during a cycle of the standard queue. From [START_REF] Brill | Level Crossing Methods in Stochastic Models[END_REF] it follows that

φ L (s) = ∞ i=1 φ I (s) i π(1 -π) i-1 = πφ I (s) 1 -(1 -π)φ I (s) ,
and then relation [START_REF] Bekker | Finite-buffer queues with workload-dependent service and arrival rates[END_REF] follows immediately. The first service time of the new queue A 1 consists of the cumulated overshoots during the first busy period of the finite queue,

A 1 = K i=1 U 1,i , ( 11
)
where the variable K denotes the number of overshoots of level 1 during a cycle of the finite queue. K is also geometric with K ≥ 1 and P(K = 1) = η. Since the U i are exponential with LST φ U (s) = λ/(λ + s), we obtain

φ A (s) = λ λ + s η 1 -(1 -η) λ λ+s = λη s + λη ,
so that the A i are exponentially distributed with rate λη. It follows that the new queue is again a G/M/1 queue with service rate λη. We can apply formula (1) for the law of its idle period D 1 . Hence, replacing λ by λη, z by z 0 and φ S by φ L in (1), we obtain

φ D (s) = λη • z 0 -φ L (s) s -λη(1 -z 0 )
, and (8) follows. But D 1 can also be seen as the sum of the idle periods of the standard queue (see Fig. 3(b)),

D 1 = I 1 + M-1 i=1 I i+1 . ( 12
)
The LST of I 1 is φ H , the transform of the conditional distribution of the idle period of the standard G/M/1 queue, given level 1 is reached during the first busy period. The variables I 2 , I 3 , . . . are i.i.d. and have the same distribution as the idle period of the standard G/M/1 queue conditioned on the event that level 1 was not reached during the last busy period. Letting φ + denote the LST of I 2 , we thus have

φ D (s) = φ H (s) • π 1 -(1 -π)φ + (s)
.

From the law of total probability we obtain φ I (s) = (1π)φ + (s) + πφ H (s) and hence (8) follows.

The M/M/1 special case

In this section we check Theorem 2 for the M/M/1 case, where the inter-arrival distribution is exponential with mean 1/μ. We assume that ρ < 1. Then φ S (s) = μ/(μ + s) and the smallest positive root of z = φ S (λ(1z)) is simply z = 1. Consequently we have

φ I (s) = λ • 1 -φ S (s) s = ρ μ μ + s ,
for the idle period of the standard G/M/1 queue. According to (3), the LST of the stationary distribution of W is given by

φ F (s) = 1 -ρ 1 -φ I (s) = 1 -ρ 1 -ρ μ μ+s which is the LST of F (x) = 1 -λ μ e -(μ-λ)
x . Hence it follows that the probabilities η and π are given by

η = μ -λ μe (μ-λ) -λ = πe -(μ-λ) and π = μ -λ μ -λe -(μ-λ) = ηe (μ-λ) .
To find the transform φ H of Z, the residual lifetime from Theorem 1, let be the number of finite renewals in a renewal process X 1 , X 2 , . . . with defective inter-arrival distribution F I and let S = i=1 X i . The random variable has a geometric distribution with P( = 0) = 1ρ, so that S is the sum of exponential random variables with rate μ and thus exponentially distributed with rate μ(1ρ) = μλ. It follows that

P(Z ≤ x) = P(S > 1, S -1 ≤ x, > 0) = P(S > 1, X ≤ x, > 0)
where X is exponential with mean 1/μ, independent of S. Hence

P(Z ≤ x) = ρP S > 1| > 0 1 -e -μx = ρe -(μ-λ) 1 -e -μx with transform φ H (s) = e -(μ-λ) λ μ+s . It follows from (9) that φ D (s) = πφ H (s) 1 -φ I (s) + πφ H (s) = λη s + μπ .
Equations ( 8) and ( 7) yield

φ L (s) = 1 - s λη φ D (s) = μπ μπ + s and then φ I (s) = φ L (s) π + (1 -π)φ L (s) = μ s + μ .
Hence I is exponentially distributed, as expected by the lack of memory property.

Idle period with set-up time

We now consider the idle period of the finite G/M/1 queue, with set-up time a ∈ (0, 1], i.e. after each busy period the workload starts with an initial value being equal in distribution to Z 1 + a. The upper diagram in Fig. 4 shows the workload process V t , together with the first cycle C 1 of that queue.

From V t we construct a new process R t , representing the time elapsed since the arrival of the customer being served. R t is obtained from the risk reserve process R t by removing the time intervals where R t = 1 (a similar construction can be found in [START_REF] Brill | Single-server queues with delay-dependent arrival streams[END_REF]).

Next we define the process W t = 1 -R a+t for 0 ≤ t < C 1a, yielding a first cycle for W t . In other words, we remove the set-up times and flip the process R t . At time t = C 1a we let the process W t restart at 1a and define the second cycle of W t via W t = 1 -R 2a+t for C 1a ≤ t < C 2 -2a. Continuing in this way, cycle by cycle, we define a regenerative process W t for t ≥ 0.

By construction, the idle periods I 1 , I 2 , . . . of the finite G/M/1 queue are identical to the overshoots of the process W t over level 1. Since the overflows occur only once in a cycle, it follows that the I i are i.i.d. with common distribution function denoted by F I .

Being a regenerative process with finite cycle mean, the process W t is stable in the sense that W t converges in distribution to some W as t → ∞. Let f denote the density of the distribution of W . Proof Let D x (C) and U x (C) denote the number of down-and up-crossings of level x by W t during the first cycle C. By level crossing theory [START_REF] Brill | Level crossings in point processes applied to queues: Single-server case[END_REF][START_REF] Cohen | On up-an down-crossings[END_REF][START_REF] Perry | A controlled M/G/1 workload process with an application to perishable inventory systems[END_REF] the longrun average number of down-crossings is given by E(D x (C))/E(C) = f (x), for all x ≥ 0. By similar reasoning as in [START_REF] Perry | A controlled M/G/1 workload process with an application to perishable inventory systems[END_REF] we find that for x < 1a the average number of up-crossings is given by

Lemma 3 The density function f fulfills the equation

f (x) = ⎧ ⎪ ⎨ ⎪ ⎩ ch(x) + ρh * f (x), 0 ≤ x < 1 -a, ch(x) -d + ρh * f (x), 1 -a ≤ x < 1, ch(x) + ρ 1 0 h(x -y)f (y) dy, x ≥ 1, (13)
E(U x (C)) E(C) = f (0) 1 -F S (x) + λ x 0 1 -F S (x -u) f (u) du = ch(x) + ρh * f (x),
and equating the two averages leads to the first line in [START_REF] Cohen | Single-server queue with uniformly bounded virtual waiting time[END_REF]. For x ≥ 1 the number of up-and down-crossings is again equal and the average number of up-crossings is f (0)(1 -F S (x)) + λ 1 0 (1 -F S (xu))f (u) du, since there are no jumps from above level 1.

If x ∈ [1a, 1) then D x (C) = U x (C) -1, since after crossing level 1 the process never returns to [0, x) during the cycle. Hence The constants can then be determined from the normalizing condition ∞ 0 f (u) du = 1.

E(D x (C)) E(C) = E(U x (C)) E(C) - 1 E(C) = f (0) 1 -F S (x) + λ x 0 1 -F S (x -u) f (u) du - 1 E(C) .
We now prove the main result of this section, relating the distribution of the idle period I of the finite G/M/1 queue to the stationary density f of W t .

Theorem 4 The distribution function of I is given by

F I (x) = 1 - f (1 + x) f (1) , ( 15
)
where f and F denote the density and distribution of W = lim t→∞ W t (the latter limit is defined in terms of weak convergence).

Proof The conditional density f c W of W -1, given that W > 1, is

f c W (x) = f (1 + x) 1 -F (1)
.

By looking at the renewal process I 1 , I 2 , . . . we conclude that f c W is also the density of the (equilibrium) forward recurrence times of that process. Hence f c W (x) = (1 -F I (x))/E(I ), and 1))/E (I ) .

F I (x) = 1 -E(I ) • f c W (x) = 1 - f (1 + x) (1 -F (
By Lemma 3, f (1) = 1/E(C) and by renewal theory, 1 -F (1) = E(I )/E(C). Consequently, 1-F (1) E(I ) = f (1) and thus [START_REF] Cohen | The Single Server Queue[END_REF] follows.

Fig. 2

 2 Fig. 2 Visualization of the proof of Theorem 1

Fig. 4

 4 Fig. 4 Construction of the processes R and W

1 0

 1 Now, since W t crosses level 1 exactly once every cycle, we have f (1+) = 1/E(C). But f (1+) = d and hence the second equation in[START_REF] Cohen | Single-server queue with uniformly bounded virtual waiting time[END_REF] follows.Solving for the renewal equationf (x) = ch(x) + ρh(x) * f in [0, 1a) (i.e.finding a solution f that fulfills the equation only on [0, 1a)) we obtainf (x) = ch * (x), x ∈ [0, 1a), where (x) = ∞ n=0 ρ n h * n (x)and h * n is the n-fold convolution of h with itself. Similarly for x ∈ [1a, 1), we have the renewal equationf (x) = (ch(x)d) + ρh * f (x). Letting d(x) = d1 {x≥1} we get f (x) = ch -d * (x) = ch * (x)d x 0 1 {x-u≥1} d (u) = ch * (x)d (x -1)-. (14)Finally for the interval [1, ∞) we havef (x) = ch(x) + ρ h(xy)f (y) dy,where f (y) is already known for y ∈ [0, 1). To find the two constants c and d, note that from d = ch(1) + ρh * f (1-) and (13) it follows that f (1-) = 0. By using[START_REF] Cohen | On up-an down-crossings[END_REF] we obtain ch * (1) = d (0) and since (0) = 1,d c = h * (1).

Acknowledgements

The authors gratefully acknowledge the comments of the anonymous referees, whose careful reading helped to improve the presentation of the paper considerably.

Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.


M G Brik
email: hbrik@fi.tartu.eeh

First-principles calculations of electronic, optical and elastic properties of ZnAl 2 S 4 and ZnGa 2 O 4

Keywords: A. inorganic compounds, C. Ab initio calculations, D. electronic structure

The optimized crystal structures, band structures, partial and total densities of states (DOS), dielectric functions, refractive indexes and elastic constants for ZnAl 2 S 4 and ZnGa 2 O 4 were calculated using the CASTEP module of Materials Studio package. Pressure effects were modeled by performing these calculations for different values of external hydrostatic pressure up to 50 GPa. Obtained dependencies of the unit cell volume on pressure were fitted by the Murnaghan equation of state, and the relative changes of different chemical bond lengths were approximated by quadratic functions of pressure. Variations of applied pressure were shown to produce considerable re-distribution of the electron densities around ions in both crystals, which is evidenced in different trends for the effective Mulliken charges of the constituting ions and changes of contour plots of the charge densities. The longitudinal and transverse sound velocities and Debye temperatures for both compounds were also estimated using the calculated elastic constants.

Introduction

A large group of compounds with spinel structure are generally described by the AB 2 X 4 formula (where A and B are di-and trivalent cations, respectively, and X is a divalent anion). There are two types of cation sites in the spinel structure: tetrahedral for A species and octahedral for B species.

Compounds from this family usually have wide band gaps, which can be attractive for various optical applications [1]. Wide band gap also makes these compounds suitable for doping with different impurities to get luminescence in particular spectral region. Optical properties of several impurity ions in these compounds were studied previously, e.g. Cr 3+ in ZnAl 2 S 4 [2][3][4][5][6][7], Fe in ZnGa 2 O 4 [8]; Mn in ZnGa 2 O 4 [9]; optical properties of ZnGa 2 O 4 nanoparticles were reported in Ref. [10]. Structural and elastic properties of ZnX 2 O 4 (X=Al, Ga, In) were calculated from first principles in Ref. [11].

In spite of considerable efforts on both experimental and theoretical studies of these compounds, previously reported results for the calculated properties of the bulk materials are somewhat different (comparison between the previous and present calculations will be given below). Besides, it turned out that ab initio calculations of the optical properties and careful studies of the pressure effects on the interionic distances and electron density distribution around ions in these compounds (especially ZnAl 2 S 4 ) are lacking, to the best of the author's knowledge. So, this was one of motivations for the present study, which reports the results of first-principles calculations of a wide range of physical properties for ZnAl 2 S 4 and ZnGa 2 O 4 single crystals. The calculated results include the optimized crystal structures (for various pressures from 0 to 50 GPa), band structures and density of states, dielectric functions and refractive indexes, elastic constants, dependencies of the chemical bond lengths and unit cell volumes on pressure, sound velocities and Debye temperatures. In addition, influence of the applied pressure on the charge density distributions, Mulliken charges, and chemical bond orders was considered.

The paper is structured as follows: in the next section the crystal structure of the considered crystals and the computational details will be described briefly, then the calculated results will be presented and compared with available experimental data and results of other calculations. The paper will be concluded with a short summary.

Crystal structure and computational details

Both considered crystals have a cubic crystal structure, space group Fd-3m with lattice constants (in Å) 8.3342 for ZnGa 2 O 4 [START_REF] Wendschuh-Josties | Neues Jahrbuch fuer Mineralogie[END_REF] and 10.009 for ZnAl 2 S 4 [13]. There are eight formula units in one unit cell. Each Al (Ga) atom is six-fold coordinated by the S (O) ions, respectively; each Zn ion is surrounded by 4 anions (O or S). The experimental crystal structural data were taken as an initial input for optimizing the crystal structure and calculations of the electronic, optical and elastic properties, as described below.

All calculations were performed in the density functional theory (DFT) framework, using the CASTEP (Cambridge Serial Total Energy Package) module [14] of Materials Studio 4.0. The total plane-wave pseudopotential method forms a basis of the CASTEP calculations. The exchangecorrelation effects were treated within the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof functional [15]. The Monkhorst-Pack scheme k-points grid sampling was set at 5×5×5 for the Brillouin zone. The plane-wave basis set energy cutoff was set at 310 eV for ZnAl 2 S 4 and 340 eV for ZnGa 2 O 4 ; ultrasoft pseudopotentials were used for all chemical elements. The convergence parameters were as follows: total energy tolerance -1×10 -5 eV/atom, maximum force tolerance 0.03 eV/nm, maximal stress component 0.05 GPa and maximal displacement 0.001 Å.

Results of calculations

Optimized crystal structures, electronic and optical properties

The optimized lattice constants are shown in 1 and2, respectively.

Coordinates of the special points of the Brillouin zone are as follows (in terms of unit vectors of the reciprocal lattice): W(0.5, 0.25, 0.75), L(0.5, 0.5, 0.5), G(0, 0, 0), X(0.5, 0, 0.5), K(0.375, 0.375, 0.750). The calculated band gaps E g were 2.699 eV (ZnAl 2 S 4 ) and 2.571 eV (ZnGa 2 O 4 ), whereas the corresponding experimental values are 3. 42 eV (ZnAl 2 S 4 , [18]) and 4.0 eV (ZnGa 2 O 4 , [19]). Such an underestimation of the calculated band gaps is related to well-known DFT limitations, namely not taking into account the discontinuity in the exchange-correlation potential [20], and is a common feature of all DFT calculations. In order to overcome such a discrepancy, the so called scissor operator [21] can be used. This operator produces a simple rigid shift of the unoccupied conduction band with respect to the valence band, and thus effectively eliminates the difference between the theoretical and experimental gap values. In this case the value of the scissor operator was 0.7 eV for ZnAl 2 S 4 and 1.5 eV for ZnGa 2 O 4 . The calculated band structures in Figs. 1 and2 are corrected with inclusion of the corresponding scissor operators. As seen from these figures, both compounds have indirect band gap, which is in agreement with previous literature data [8,18]. However, there exists some difference in how this "indirect character" of the band gap is formed in both compounds. The top of the valence band for ZnAl 2 S 4 is flat near the G point, whereas in ZnGa 2 O 4 the top of the valence band and the bottom of the conduction band are simply realized at different points of the Brillouin zone. The "zoomed" views of the band gaps for both crystals clearly illustrate this difference. Since both energy gaps are indirect, the phonon contribution to the absorption processes should be important.

Composition of the calculated energy bands can be resolved with the help of partial density of states (PDOS) and total density of states (DOS) diagrams. Figs. 3 and4 show the PDOS and total DOS for ZnAl 2 S 4 and ZnGa 2 O 4 , respectively. These diagrams allow concluding that the conduction band in ZnAl 2 S 4 is about 5 eV wide and is formed by the Al 3s and 3p states, which are hybridized with the S 3p states and Zn 4s, 4p states. The valence band is wide -about 7 eV -and consists of two sub-bands, clearly seen in the band structure as well: the upper one (between -6 and 0 eV) is a mixture of the S 3p states and Al 3s, 3p states. The lower one is narrow (between -7 and -6 eV) and is due to the completely filled Zn 3d states. Another band between -14 and -11 eV is created by the S 3s states with a minor contribution coming from the Al 3s, 3p states as well. One of the main optical characteristics of a solid is its complex dielectric function 2 1

) (

ε ε ω ε i + = . The imaginary part ) Im( 2 ε ε = of dielectric function is calculated in CASTEP
numerically by a direct evaluation of the matrix elements between the occupied and unoccupied electronic states. The real part ) Re(

1 ε ε = of ) (ω ε
is calculated then using the Kramers-Kronig transform. Such an approach, of course, has certain limitations. First of all, the local field effects (related to the fact that the electric field experienced at a particular site of a crystal lattice is screened by the polarizability of the atom at this site) are neglected. Also phonon contributions to the optical spectra, which are especially important for the crystals with indirect band gap, are also not taken into account, when calculating optical properties in CASTEP. However, even with these limitations the calculated spectra give reasonable agreement with experimental results (see, for example, calculations of absorption spectra for PbWO 4 [22,23], TiO 2 [24] etc).

The dielectric functions for ZnAl 2 S 4 and ZnGa 2 O 4 are shown in Figs. 5 and6, respectively. The instrumental smearing of 0.2 eV was used to model the broadening effects. The regions, in which the imaginary part 2 ε is different from zero, can be related to the absorption spectrum (which, in this case, will be due to the electronic transitions from the valence band to the conduction band). So, from these A further insight into the electronic properties of these compounds and chemical bond formation can be gained by considering the charge densities distributions (Figs. 7 and8; one primitive cell is shown). As seen from these figures, in ZnAl 2 S 4 there is strong hybridization between the S and Zn atoms, whereas interaction between the Al and S atoms is not that pronounced. On the contrary, in ZnGa 2 O 4 oxygen atoms interact strongly with both Zn and Ga atoms as well.

Different character of charge density distributions can be also seen in different values of the

Mulliken charges [25] for ions in both compounds. They are as follows: Zn (+0.96), Ga (+1.

C C C = = .
Calculations of these constants were performed for zero-th external pressure and optimized crystal structures; the results are shown in Table 2. For both compounds the values of the 11 C constant are approximately 2-2.2 times greater than each of two remaining constants, which is in line with a wellknown observation that for cubic crystals elastic constants decrease with increase of the lattice parameter [START_REF] Kittel | Introduction to Solid State Physics[END_REF]. It is also easy to notice that all elastic constants for ZnAl 2 S 4 are about twice as small as those for ZnGa 2 O 4 , suggesting that the same external pressure should produce greater deformations in the former host.

Variations of the interionic distances with external hydrostatic pressure were studied by optimizing crystal structures for both considered crystal in the pressure range from 0 to 50 GPa. Fig. 9 shows the relative 0 y y variation of the Ga-O, Zn-O distances in ZnGa 2 O 4 and Al-S, Zn-S distances in ZnAl 2 S 4 crystals (y stands for the interionic distance at some pressure P, whereas y 0 is the same distance at ambient pressure). The Ga-O chemical bond appears to be the stiffest, whereas the Zn-S bond is the weakest among the four considered ones.

Relative variations of these interionic distances were approximated by the quadratic leastsquares fits 3. They can be readily used for estimations of the interionic distances for any arbitrary pressure up to 50 GPa, which is especially important for doped crystals, when the pressure dependent luminescence of dopants is studied.

Dependence of volume V on pressure P for solids can be modeled by the Murnaghan equation [START_REF] Murnaghan | Proc. Natl. Ac. Sci[END_REF] '

1 0 0 0 0 ' 1 B B P B V V - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = (1)
where V 0 is the volume at ambient pressure, 0 B and ' 0 B are the bulk modulus and its pressure derivative, respectively. The 0 V V ratios for unit cells of both crystals were calculated for pressures 0, 10, 20, 30, 40, and 50 GPa (Fig. 10, symbols). The least-squares fits to Eq. ( 1) are shown in Fig. 10

C C =

[29] is practically fulfilled. S. Haussühl [START_REF] Haussühl | [END_REF] formulated that exactly this relation holds true for atoms, which interact only by central forces. It has been also noted [START_REF] Haussühl | [END_REF][31][32] Dependence of the calculated band gaps for both compounds on pressure is presented in Fig. 11.

Linear fits (also shown in the Figure ) describe well the behavior of the calculated band gaps. However, the value of E g increases with pressure for ZnGa 2 O 4 and exhibits opposite behavior for ZnAl 2 S 4 .

Usually, band gaps increase with pressure, which also can be interpreted in terms of a blue shift of charge transfer transitions (which are, in other words, band-to-band transitions) with applied hydrostatic pressure [34]. Nevertheless, examples of decreased band gaps with applied pressure have been reported for a number of compounds: PbS [35], Li 3 N [36], BN nanotubes [37], MnIn 2 S 4 [38],

Li 4 CaB 2 O 6 [39]. The reason for such a behavior can be either a high covalency of chemical bond, or significant re-distribution of electron density in the space between atoms [39]. Speaking about the crystals considered in the present paper, it can be mentioned that ZnAl and8; in all figures one primitive cell is depicted for the sake of simplicity).

Enhanced interactions in the Zn-O, Ga-O, Zn-S pairs are clearly seen form the shown contour plots.

It is also interesting to follow dependence of the Mulliken charges on pressure, which is represented in Table 4. In ZnGa 2 O 4 the absolute values of the oxygen and gallium Mulliken charges are increasing with pressure, whereas the zinc Mulliken charge is decreasing. Such a result can be interpreted as re-distribution of the electron densities between these ions, namely, a shift of the electron density (negative charge) from gallium ions towards oxygen ions and from oxygen ions towards zinc ions, with the former shift dominating over the latter one.

In ZnAl 2 S 4 the absolute values of the Mulliken charges for all three ions are decreasing, which can be thought of as a shift of electron density (negative charge) from sulfur ions towards zinc and aluminum ions. Additional information about the nature of the chemical bond between atoms can be extracted by using the concept of the Mulliken bond order. Although the absolute values of the bond order depend on the choice of the basis set, and thus their validity should not be overestimated [40],

their relative values, calculated consistently for similar or even identical systems under various external conditions, like pressure, for example, can reveal some trends in behavior of particular chemical bonds and ions involved into their formation. A greater value of this parameter indicates the covalent character of the considered bond, whereas lower values of the bond order stress the ionic nature of the bond [40]. Analysis of the calculated Mulliken bond orders (Table 5) indicates that the Ga -O bond is the most ionic, whereas the Zn -S bond is the most covalent. On average, chemical bonds in ZnAl 2 S 4 are more covalent than in ZnGa 2 O 4 -this conclusion is also in line with the results of elastic constants estimations. Increase of pressure leads to a monotonic increase of the bond order (increased covalency) for each individual bond.

The results of calculation of the elastic constants can be used further to estimate one of the most important thermal characteristics of the studied compounds, the Debye temperature D θ . This parameter is very important for an analysis of the specific heat, thermal conductivity and melting temperature. It can be estimated in the framework of the Debye-Grüneisen model by means of the following simple equation [41]

m A D v M N n k h 3 / 1 4 3 ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ρ π θ , (2)
where h and k are the Planck's and Boltzmann's constants, respectively, N A is the Avogadro's number, ρ is the crystal's density, M is the molecular weight, n denotes the number of atoms per unit cell (56 for both studied compounds), m v is the mean sound velocity expressed in terms of the longitudinal l v and transversal t v sound velocities as follows:

3 / 1 3 3 1 2 3 1 - ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = l t m v v v . (3)
The l v and t v sound velocities are calculated as [START_REF] Schreiber | Elastic Constants and Their Measurements McGraw-Hill[END_REF][START_REF] Kanoun | [END_REF][44]

ρ 3 4 3 G B v l + = , ρ G v t = , (4)
with B being the bulk modulus and 2

R V G G G + =
the isotropic shear modulus, which is calculated as the average value of the Voigt's shear modulus G V (an upper limit for G values) and the Reuss's shear modulus G R (a lower limit for G values). Analytical expressions for G V and G R are as follows:

3 C C C G C C C G R V + - = + - = , (5)
and the values of the elastic constants

44 12 11 , , C C C
are given in Table 2.

The density of ZnAl 2 S 2 is 3.29 g/cm 3 [45] and the corresponding value for ZnGa 2 O 4 is 6.17 g/cm 3 [46].

A summary of calculated values obtained from Eqs. ( 2) -( 5) is presented in Table 6.

The ZnGa 2 O 4 shear moduli, like other elastic constants, are approximately twice as large as their ZnAl 2 S 4 counterparts. Also the sound velocities for ZnGa 2 O 4 are somewhat greater than for ZnAl 2 S 4 .

The values of Debye temperature decrease with increasing the lattice constants and anion's atomic number and weight; therefore, the highest acoustic phonons frequency should be greater in ZnGa 2 O 4 .

Conclusion

The first-principles calculations of the structural, electronic, optical and elastic properties for 3) are shown by solid lines.

  The conduction band in ZnGa 2 O 4 is considerably wider (about 7 eV) than in ZnAl 2 S 4 and is composed of Zn 4s, 4p states, Ga 4s, 4p states and O 2p states. The valence band is somewhat narrower in ZnGa 2 O 4 (about 6 eV in comparison with 7 eV in ZnAl 2 S 4 ) and does not exhibit any subband structure. 4s states of Ga occupy the bottom of the valence band, whereas the rest of the valence band consists of the Zn 3d states and O 2p states. Two remaining narrow bands at low energies are due to the Ga 3d states (at about -12.5 eV) and O 2s states (between -18 eV and -16 eV). It is also noticeable that the Zn 3d states in ZnGa 2 O 4 are spread over a wider region (almost 6 eV), whereas in ZnAl 2 S 4 the Zn 3d states are localized in vicinity of -7 eV, at the very bottom of the valence band.

figures absorption starts at about 3

 3 figures absorption starts at about 3 eV for ZnAl 2 S 4 and 4 eV for ZnGa 2 O 4 , in agreement with the band gap estimations. The real part of ) (ω ε in the limit of zero energy (or infinite wavelength) is equal to the square of the refractive index n. From Figs. 5 and 6 it is possible to estimate the value of n to be 2.24 for ZnAl 2 S 4 and 1.65 for ZnGa 2 O 4 . Unfortunately, no experimental values of the refraction index were found in the literature, so these estimations remain to be purely theoretical. Difference in electronic and band structure of these two compounds is also manifested in different shapes of the imaginary part of the dielectric function. Different band gaps are obviously related to different absorption edges. Besides, the widths of non-zero part of imaginary component of dielectric functions (in other words, the total width of the absorption spectra) are also different. It is directly linked with the widths of the conduction bands: a wider conduction band in ZnGa 2 O 4 corresponds to a more extended absorption spectrum. In addition, a greater number of absorption peaks in ZnGa 2 O 4 can be attributed to more significant contribution of the Zn 3d states to the total absorption, since they are distributed almost over the whole valence band (in ZnAl 2 S 4 they are localized at the bottom of the valence band).

  are plotted by solid lines in Fig. 9; they excellently follow the calculated values shown by different symbols for each chemical bond. The values of the approximation coefficients are listed in Table

  2 S 4 is a more covalent compound than ZnGa 2 O 4 , with a quite considerable contribution of the sulfur ions into the conduction bands (hybridization with Zn and Al ions), which falls in line with some of the above given examples. Finally, variation of the chemical bonds lengths with pressure should affect the charge densities distributions as well. Figs. 12 and 13 show how the charge density distributions are modified at the highest considered pressure of 50 GPa for ZnAl 2 S 4 and ZnGa 2 O 4 , respectively (these figures should be compared with Figs. 7

ZnAl 2 S 4 and

 24 ZnGa 2 O 4 have been performed using the CASTEP module of Materials Studio package. Both crystals are characterized as indirect band compounds. Structure of the calculated energetic bands was analyzed using the partial density of states diagrams. Although the structure of both crystals is the same, their electronic, elastic and bonding properties are quite different. It was shown that the conduction band is considerably wider in ZnGa 2 O 4 , whereas the valence band is spread over a wider energy interval in ZnAl 2 S 4 . The Zn 3d states are localized at the very bottom of the valence band in ZnAl 2 S 4 , whereas the same Zn 3d states occupy practically the whole valence band in ZnGa 2 O 4 . Special attention has been paid to pressure effects and their influence on the electronic and structural properties. Variations of the Ga -O, Zn -O, Al -S, Zn -S chemical bonds lengths with pressure were approximated by the quadratic functions of applied pressure. It was shown that the Ga -O bond is the stiffest, and the Zn -S is the weakest one. Dependencies of the unit cell volumes on pressure were fitted by the Murnaghan equation of state; determined from these fits the bulk moduli for both compounds are very close to those obtained from available in CASTEP direct calculations of elastic constants.Analysis of the chemical bonds in both crystals (with calculations of the Mulliken bond orders)showed ZnAl 2 S 4 to be considerably more covalent than ZnGa 2 O 4 . In addition, chemical bonds in the tetrahedral complexes formed by Zn ions and its neighbors (oxygen or sulfur) are more covalent than the chemical bonds in the octahedral complexes formed by the (Ga,O) and (Al, S) ions. Detailed analysis of the applied hydrostatic pressure effects has shown that it not only decreases the interionic distances, but also causes considerable re-distribution of the electron density, which was confirmed by calculating effective Mulliken charges for different pressures and following different trends of their variation. Finally, the calculated elastic constants were used to estimate the sound velocities and Debye temperatures for both compounds.

Fig. 3 .

 3 Fig. 3. Partial and total densities of states for ZnAl 2 S 4 . From top to bottom: sulfur, aluminum, zinc PDOS and total DOS.

Fig. 4 .

 4 Fig. 4. Partial and total densities of states for ZnGa 2 O 4 . From top to bottom: oxygen, gallium, zinc PDOS and total DOS.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Calculated dielectric function of ZnAl 2 S 4 .

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. A contour plot of the charge-density distribution in ZnAl 2 S 4 (a normal to the shown slice is [0.707, -0.707, 0.000])

4 Fig. 10 .

 410 Fig. 10. Dependence of relative change of volume on pressure for ZnGa 2 O 4 and ZnAl 2 S 4 . The calculated values are shown by symbols, and the fits to the Murnaghan equation (1) by solid lines.

Fig. 11 .

 11 Fig. 11. Pressure dependence of the calculated band gaps E g for ZnGa 2 O 4 (squares) and ZnAl 2 S 4 (circles). Linear fits and their equations are also shown.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. A contour plot of the charge-density distribution in ZnAl 2 S 4 (a normal to the shown slice is [0.707, -0.707, 0.000]; pressure 50 GPa)

Table

  ZnGa 2 O 4 were calculated by other authors as well (with rather wide range of the calculated lattice constants), whereas similar calculations for ZnAl 2 S 4 are scarce. Calculated band structures of ZnAl 2 S 4 and ZnGa 2 O 4 are shown in Figures

, in comparison with experimental findings and results of other available calculations for ZnGa 2 O 4 . As seen from this Table, agreement between the results of the present calculations and experimental data is very good. Structural properties of

  22), O (-0.85) in ZnGa 2 O 4 and Zn (+0.79), Al (+0.72), S (-0.56) in ZnAl 2 S 4 . Increased positive value of the Mulliken charge for Zn ions in ZnGa 2 O 4 can be attributed to the stronger mixture with surrounding oxygen ions, when the negative charge is shifted towards oxygen ions, thus resulting in a greater absolute value of the anion's (oxygen) charge as well, when compared with that of sulfur in ZnAl 2 S 4 .

	This is also in line with a shorter calculated Zn -O distance (1.99713 Å) in comparison with Zn -S
	separation of 2.35347 Å. Additional confirmation of the obtained result can be based on the optical
	electronegativities of oxygen and sulfur ions, which would enhance covalency when going from O 2-to
	S 2-ligands [26].
	3.2. Elastic properties and pressure effects
			Elastic properties of cubic crystals are described by three independent elastic moduli
	11 C	,	12 C	,	C	44	[27], with the following identities	11 C	=	C	22	=	33 C	,	12 C	=	C	23	=	31 C	,	44	55	66

Table 2 )

 2 by solid lines. . Estimated pressure derivatives of the bulk moduli are also within the typical range for solids. It is easy to note from Table2that the Cauchy condition for cubic crystals

	From these approximations, the values of 0 B and ' 0 B for ZnAl 2 S 4 are 78.91 ± 0.47 GPa and 3.97 ± 0.04,
	respectively and for ZnGa 2 O 4 163.46 ± 1.31 GPa, 4.40 ± 0.09, accordingly. The extracted from the
	plots values of the bulk moduli 0 B are very close to those obtained as the results of the elastic constants
	calculations (44 12

Table 1 .

 1 Crystal lattice constants a and unit cell volume V for ZnGa 2 O 4 and ZnAl 2 S 4 crystals This work; c Ref.[11]; d Ref.[16]; e Ref.[19]; f Ref.[17]; g Ref.[13]

					ZnGa 2 O 4				ZnAl 2 S 4
		Exp. a	Calc. b	Calc. c	Calc. d	Calc. e	Calc. f	Exp. g	Calc. b
	a, Å	8.3342	8.4412	8.4063	8.38	8.174	7.977	10.009	10.0340
	V, Å 3 578.884 601.468 594.039	588.480	546.140	507.597	1002.702	1010.235
	a Ref. [12]; b

Table 2 .

 2 Calculated elastic constants (all in GPa, except for non-dimensional Poisson ratios). The values in parenthesis for ZnGa 2 O 4 are the results of calculations from Ref.[11].

	Compound	11 C	12 C	C	44	Bulk modulus B 0	Young Modulus	Poisson ratio
	ZnAl 2 S 4	112.42 ± 2.05	50.84 ± 3.08	59.63 ± 6.13	71.37 ± 2.17	80.75	0.3114
	ZnGa 2 O 4	249.45 ± 3.54 (228)	121.81 ± 0.95 (120)	115.50±2.45 (107)	164.36 ±1.34 (156)	169.52	0.3281

Table 3 .

 3 Coefficients of quadratic approximations O 4 and ZnAl 2 S 4 .

	y	0 y	=	a	+	1 b	P	+	2 b	P	2	(P is pressure in GPa) of

Table 4 .

 4 Mulliken charges of individual ions in ZnGa 2 O 4 and ZnAl 2 S 4 as a function of external

	pressure
	Pressure,		ZnGa 2 O 4			ZnAl 2 S 4
	GPa	Zn	Ga	O	Zn	Al	S
	0	0.96	1.22	-0.85	0.79	0.72	-0.56
	10	0.93	1.25	-0.86	0.79	0.72	-0.56
	20	0.91	1.27	-0.86	0.78	0.71	-0.55
	30	0.90	1.29	-0.87	0.77	0.70	-0.54
	40	0.90	1.30	-0.88	0.76	0.70	-0.54
	50	0.89	1.32	-0.88	0.76	0.69	-0.53

Table 5 .

 5 Mulliken bond orders for different chemical bonds in ZnGa 2 O 4 and ZnAl 2 S 4 as a function of

	external pressure
	Pressure,	ZnGa 2 O 4		ZnAl 2 S 4
	GPa	Zn -O	Ga -O	Zn -S	Al -S
	0	0.43	0.27	0.50	0.45
	10	0.45	0.29	0.58	0.47
	20	0.47	0.30	0.61	0.49
	30	0.48	0.31	0.64	0.50
	40	0.49	0.33	0.67	0.52
	50	0.50	0.34	0.70	0.53

Table 6 .

 6 Calculated shear moduli, sound velocities and Debye temperatures for ZnAl 2 S 4 and ZnGa 2 O 4

	Compound	G V , GPa	G R ,	G, GPa
	ZnAl 2 S 4	48.09	43.37	45.73	6342.4	3728.4	4132.5	943.8
	ZnGa 2 O 4	94.83	87.24	91.04	6805.2	3841.2	4272.6	1170.8

l v , m/s t v , m/s m v , m/s D θ , K

Acknowledgment:

The author thanks Prof. Ü. Lille (Tallinn University of Technology) for allowing to use the computational facilities and Materials Studio package. Financial support from Estonian Science Foundation (Grants No. 7456, JD69, 6999 and 6660) is also appreciated.