Marko P J Punkkinen

Se Qing-Miao Hu

Kyun Kwon

Borje Johansson

János Kollár

Levente Vitos
email: levente@kth.se

Q.-M Hu

S K Kwon

Surface properties of 3d transition metals

Keywords: surface properties, density-functional methods, magnetism user supplied): surface energy, surface stress, magnetic transition metals surface relaxation, surface energy, surface stress, magnetic transition metals, density functional theory

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Surface energy (γ) and surface stress (τ ij ) are two fundamental physical parameters of surfaces. For solids, they differ by the excess or residual surface stress τ ij -γδ ij = ∂γ/∂ ij , where ij stands for the strain within the surface plane. For high symmetry surfaces, often the scalar surface stress is used (τ ), which corresponds to the average stress along in-plane directions. For liquids, the surface free energy and the surface stress are equal each other due to the fact that in this case the surface energy does not change when the surface is strained, i.e. ∂γ/∂ = 0.

Knowledge of the geometry of solid surfaces is amongst the basic questions of surface science. This is because in general the surface parameters exhibit strong structure dependence [1]. Transition metal surfaces are of particular interest, since they act as catalysts in various chemical reactions. It is well known that real surfaces of pure metals can adopt different structures from those of ideally truncated crystals. Upon layer relaxation only rigid inward or outward displacement of the atomic layers occurs, while in case of reconstruction, the displacement of atoms may alter the two-dimensional symmetry of the surface. Experimental studies have † Deceased. demonstrated that the surface layer of clean transition metal surfaces relaxes inward [START_REF] Martienssen | Springer Handbook of Condensed Matter and Materials Data[END_REF], i.e. the inter-layer distance between the topmost two atomic layers is smaller than that of the bulk. Outward expansion of the top layer has also been found for some surfaces of noble metals. The top layer relaxation is often accompanied by relaxation of the subsurface atomic layers, resulting in an oscillatory behavior of multi-layer relaxation. Several theoretical models have been proposed to explain the surface relaxation of transition metals. Generally, the magnitude of relaxation is larger for rough surfaces than for smooth ones. The model, proposed by Finnis and Heine [START_REF] Finnis | [END_REF] based on Smoluchowski smoothing [4], states that when cutting a perfect crystal, charge redistribution gives rise to an inward electrostatic force on the top layer nuclei. Since this effect increases with surface roughness, the above model confirms the general relationship between relaxation and surface roughness. It also explains the contraction found for most transition metal surfaces, but it fails to describe the noble metal surface relaxation. Pettifor [5] suggested that the crystal geometry is determined by the balance of the negative pressure of the localized d bonds and the homogeneous positive pressure of the sp electrons. According to Heine and Marks [6], at metal surfaces the mobile sp electrons can flow into the vacuum, while the d bonds between the first and second layers remain practically unchanged. Therefore, the d electrons give rise to an inward force on the surface atoms, which is proportional to the strength of the d bonds. Because the strength of the d bonds shows a well known parabolic behavior across the transition metal series, the magnitude of the top layer relaxation is also expected to follow a parabolic trend. This model predicts positive (outward) relaxation for noble metal surfaces. A third approach, introduced by Kádas et al. [7] using ab initio results based on the exact muffin-tin orbitals method [8,9], gives a quantitative foundation for the model by Heine and Marks. According to that, the layer relaxation (δ) is proportional to (V /B)∆n sp , where ∆n sp is the depletion of sp electrons in the surface layer relative to bulk, V and B are the atomic volume and bulk modulus. Since ∆n sp follows a parabolic trend across the nonmagnetic transition series [7] the above expression would indeed dictate a parabolic shape for δ as a function of d-occupation number. However, due to the large equilibrium volume and the rather small bulk modulus of the early transition metals, the multiplier V /B is relatively large at the beginning of the series resulting in large layer relaxations for those elements.

One possible way to approach the problem of surface relaxation is to treat the surface as a two-dimensional defect with a potential which is screened towards the bulk metal. It is well known that free surfaces cause Friedel oscillations in the charge density towards the bulk metal [10]. The periodicity of the Friedel oscillation is twice the Fermi wave vector k F and it is generally incommensurate with the relevant lattice constant (inter-layer distance) of the metal. This causes the layer relaxation to fall off quickly towards the bulk. However, if the lattice constant and 2k F are commensurate, then the layer relaxation and the Friedel oscillation strengthen each other, which results in a slowly decaying layer relaxation towards the bulk. Striking examples are the hexagonal close-packed Tc and Re [1,11,12], where the magnitude of the deep layer relaxations were found to be similar to those of the top-layers.

Both the surface energy and surface stress are sensitive to the actual surface geometry. In particular, the layer relaxation dependence of the surface stress is closely linear whereas that of the surface energy quadratic [1]. Due to the linear term, any accurate surface stress calculation must be preceded by a precise surface geometry study. Kwon et al. [1] connected the layer relaxation dependence of γ and τ to the surface elastic constants. Nichols et al. [13] According to a simple model described by Ibach [14], upon cleaving a metal surface the electronic charge density of the broken bonds is redistributed between surface atoms and their backbonds. Hence, surfaces should possess inward layer relaxation as a result of the strengthening of the backbonds of surface atoms. On the other hand, since lattice constraint by sub-surface layers hinders in-plane relaxation of surface atoms, the increased charge density between surface atoms is expected to result in tensile (positive) surface stress on metal surfaces. First-principles calculations performed for non-magnetic transition metals [1,11,15] indeed confirmed that the surface stress on clean metal surfaces is tensile. However, a recent density functional study on 3d metals [16] demonstrated that magnetism can overwrite this normal behavior in magnetic transition metals resulting in compressive surface stress for some of the close-packed surfaces. The atomic-scale mechanism behind this anomalous surface stress has been shown to be the enhanced magnetism near the free surfaces, which favors a larger lattice constant as compared to the bulk. In this paper, we provide further details of our study performed for 3d metals [16] by presenting ab initio data on the top layer relaxation, surface energy and surface stress. In Section 2 we outline the method of the surface energy and stress calculation, while in Section 3 we summarize and discuss the results.

Computational method

Surface energy and stress

The surface stress is defined to be the reversible work per unit area to stretch the surface elastically. It can be expressed by the Shuttleworth equation [17]

τ ij = 1 A ∂Aγ ∂ε ij = γδ ij + ∂γ ∂ε ij , (1)
where τ ij denotes the components of the surface stress tensor, A is the surface area, ε ij stands for a deformation tensor element, γ denotes the surface energy, and δ ij is the Kronecker delta. The surface energy is defined as the reversible work per surface area to create a surface, and may be determined as

γ = E s -E b A , (2)
where E s is the total energy of the surface region and E b is the total energy of a bulk region, both of them referring to the same number of atoms. Throughout the paper, upper s index refers to surface and upper b index to bulk quantities. In equations ( 1) and ( 2), all quantities are evaluated for the unstrained lattice, i.e. for ε ij = 0. Since around the equilibrium, the variation of the bulk energy E b (ε) as a function of lattice strain vanishes in linear order, in fact it is the linear strain dependence of E s (ε) that leads to nonzero stress. A general expression of the layer relaxation dependence of the above surface parameters may be formulated as follows. Here for the sake of simplicity we assume that only the top layer is relaxed, viz. there is only one inter-layer distance λ different from the ideal bulk value. Using notation ε for the biaxial in-plane strain [15,18] and δ = (λλ 0 )/λ 0 for the deviation relative to the equilibrium (surface or bulk) inter-layer separation λ 0 , the fully relaxed slab energy relative to the bulk energy can formally be expanded as

E s (ε, δ) -E b (ε, δ) ∼ E s (0, 0) -E b (0, 0) + + A 0 τ 0 ε + 1 2 γ ε ε 2 + 1 2 γ δ δ 2 + τ , 0 εδ + ... , (3)
where τ 0 , γ ε , γ δ and τ , 0 are constants and A 0 is the area of the undistorted surface. Note that in expansion (3) the pure linear term in δ vanishes because the slab energy corresponds to the layer relaxed structure (viz. ∂E s (ε = 0, δ)/∂δ = 0 for δ = 0) and the bulk energy E b (ε, δ) contains only second and higher order terms in δ and ε. Using the above expression in connection with definitions ( 1) and (2), in the limit ε → 0 we find

γ(δ) = γ(0) + 1 2 γ δ δ 2 + ... (4)
and

τ (δ) = τ (0) + τ , 0 δ + ..., (5)
where γ(0) = (E s (0, 0) -E b (0, 0))/A 0 and τ (0) = τ 0 . Equations ( 4) and ( 5) are equivalent to those derived by Kwon et al. [1]. According to these expressions, the surface energy depends on the layer relaxation only in the second order but the surface stress due to the mixed term in (3) exhibits a much stronger δ dependence.

Numerical details

In order to calculate the surface stress, we consider a slab formed by N atomic layers parallel to the surface. For surface calculations, the slab is embedded in vacuum, while for bulk calculations the slab is periodically repeated along the direction perpendicular to the atomic layers. The bulk is relaxed to its equilibrium with equilibrium inter-layer distance λ 0 . Next, for the slab with in-plane lattice parameters fixed to those of the bulk, we perform multi-layer relaxation (apart from a few central layers, where the inter-layer distance is kept at λ 0 ) and find the λ inter-layer distances which minimize the total slab energy. In other words, the only allowed relaxation is the one which is perpendicular to the surface plane, i.e. we do not consider in-plane relaxation. Starting from the so relaxed slab geometry, we elongate the lattice vectors lying in the surface plane by ε but keep the interlayer distances fixed. For this distortion, the deformation tensor has the form

ε ij =   ε 0 0 0 ε 0 0 0 0   . (6)
The total energies of the surface and bulk slabs, i.e. E s and E b , are computed for five different deformations (ε = 0, ±0.01 and ±0.02) and fitted by a quadratic polynomial, i.e.

F

E s/b (ε) ≈ E s/b (0) + c s/b 1 ε + c s/b 2 ε 2 . (7)
Here E b (0) is the total energy of the undistorted bulk and E s (0) is the total energy of the undistorted slab (corresponding to δ = 0). The surface stress is determined from the linear coefficients of the slab and bulk energies at ε = 0 (unstrained lattice), viz.

τ = c s 1 -c b 1 4A , (8)
where the factor 4 arises from the two surfaces of the slab and the two deformed in-plane lattice vectors. In the case of low symmetry surfaces, equation (8) gives the average of the two main stress tensor components. Obviously, at equilibrium the bulk stress should vanish (c b 1 = 0) and thus it is expected that the surface stress reduces to c s 1 /(4A). However, due to different numerical parameters in bulk and slab calculations, this is not always the case and thus equation ( 8) is recommended in stress calculations based on the total energy.

In the present application, the free surfaces were modeled by periodically repeated slabs separated by vacuum layers. At ambient conditions, V and Fe have the body centered cubic (bcc) structure; Sc, Ti, Co and Zn have the hexagonal close-packed (hcp) structure, and Ni and Cu are face centered cubic (fcc) metals. Chromium has an incommensurable antiferromagnetic state, which however can be approximated by a commensurable antiferromagnetic state with B2 structure [19]. In the following, when referring to bcc Cr we mean in fact B2 Cr. The ferromagnetic state was used for Fe, Co and Ni. For each system, the surface properties were computed for the most stable facet and for some systems also for the second most stable facet. The slabs were formed by 8 atomic layers for the fcc (111) and (100) surfaces, 12 atomic layers for the bcc (110) and hcp (0001) surfaces, and 16 atomic layers for the bcc (100) surface. For these systems, the thickness of the vacuum layer was set to 22 -23 Å. Manganese adopts a complex antiferromagnetic structure (α-Mn) with 58 atoms in the unit cell. The (110) facet of α-Mn was reported to be the most stable surface [20]. This surface was modeled by a slab containing 130 Mn atoms separated by a vacuum layer of thickness 12.2 Å.

All numerical calculations were performed using density functional theory [21] in combination with the generalized gradient approximation [22]. The surface geometry, surface energy, and surface stress were computed using the projector augmented wave (PAW) method as implemented in the Vienna Ab initio Simulation Package (VASP) [23,24]. In addition to the PAW calculations, we also carried out a series of tests for the unrelaxed bcc (001) surfaces of V, Cr, Mn and Fe. In these tests we employed two alternative ab initio total energy approaches. The first tool is based on the exact muffin-tin orbitals (EMTO) method [8,9,[START_REF] Vitos | Computational Quantum Mechanics for Materials Engineers[END_REF][START_REF] Andersen | Lectures on Methods of Electronic Structure Calculations[END_REF]] and the second one on the full-potential augmented plane wave with local orbitals (FP APW+lo) method as implemented in WIEN2k package [START_REF] Schwarz | [END_REF]28].

The EMTO method [START_REF] Andersen | Lectures on Methods of Electronic Structure Calculations[END_REF] is a screened Korringa-Kohn-Rostoker type of method [29][30][START_REF] Zabloudil | Electron Scattering in Solid Matter A Theoretical and Computational Treatise[END_REF]. In the present application the one-electron Kohn-Sham equation is solved within the scalar-relativistic approximation (omitting the spin-orbit term) and using the Green's function technique. The one-electron potential is represented by optimized overlapping muffin-tin potential spheres. By using overlapping spheres, one describes more accurately the exact crystal potential, compared to the conventional muffin-tin or non-overlapping approach [START_REF] Vitos | Computational Quantum Mechanics for Materials Engineers[END_REF]32]. In combination with the Table 1. Calculated bulk lattice parameters (in Å) and optimized c/a ratios for the 3d metals. For comparison, the experimental data are given in parentheses [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Lawson | [END_REF]. full charge density technique [START_REF] Kollár | From ASA Towards the Full Potential[END_REF][34][35] for the total energy calculation [9,[START_REF] Vitos | Computational Quantum Mechanics for Materials Engineers[END_REF], the EMTO method is also suitable to describe accurately the total energy with respect to anisotropic lattice distortions. The EMTO method has been applied successfully in the theoretical study of the thermophysical properties of metallic alloys [9,[START_REF] Vitos | Computational Quantum Mechanics for Materials Engineers[END_REF][36][37][38][39][40][41] and complex oxides [42][43][44][45]. The FP APW+lo method treats the potential and charge density without any shape-approximations inside the muffin-tin spheres and in the interstitial region [46]. The wave-function is expressed as the augmented plane wave (APW) plus local orbitals for the semicore states [28]. The APW with each radial wave function is computed at a fixed linearization energy. The local orbitals are evaluated at the same energy as the APW's and vanish at the muffin-tin boundary. The APW+lo scheme leads to significantly smaller basis sets and therefore reduces the computation time as compared to the linearized augmented plane wave method without losing accuracy.

Metal

In PAW calculations, for all slabs, the in-plane lattice constant was fixed to the theoretical bulk equilibrium value, and the inter-layer distances, except for the central layers, were allowed to relax to their equilibrium values. The optimization of the atomic structure was performed using conjugate-gradient minimization of the total energy with respect to the atomic coordinates. The atoms were relaxed until the remaining forces were less than 20 meV/ Å. The Brillouin zone sampling was performed by the Monkhorst-Pack scheme [47]. The chosen plane wave cutoff energy (450-500 eV) in PAW calculations and the k-mesh in EMTO and FP APW+lo calculations assured ∼ 1 meV per atom convergence in the total energy. For each surface, the mean surface stress was determined from the variation of the surface energy upon biaxial in-plane strain [15,18].

Results and discussion

Layer relaxation

The calculated bulk equilibrium lattice parameters and c/a-ratios for the 3d metals including Zn are listed in Table 1. They are in good agreement with former theoretical results obtained using the same exchange-correlation functional [50] as well as with the quoted experimental data [START_REF] Kittel | Introduction to Solid State Physics[END_REF][START_REF] Lawson | [END_REF].

In Table 2, we give our results for the multi-layer relaxation obtained for some low-index surfaces of 3d transition metals. Examining the layer relaxations, we see that for most 3d metals the relaxation decays relatively fast with the distance from the surface. For most of the systems, especially where the expected crystal structure is not altered by magnetic effects, the top layer relaxation follows similar trend as those observed for the 4d and 5d metals [11]. Exceptions are the bcc (110) surfaces of Cr and Fe, where the calculated relaxations of -1.95 % and -0.05 % are significantly smaller than those obtained for the other metals from the same groups. In particular, for the top layer relaxation of the close-packed surfaces of W and Os, Zolyomi et al. [11] predicted -3.76 % and -3.72 %, respectively, whereas for Mo and Ru Kwon et al. [1] obtained -4.74 % and -3.96 %, respectively. Note that the same is true for the bcc (100) surface of Cr when compared to that of Mo or W [1,11]. Manganese is a somewhat special case. If atomic layers are defined for the (110) surface of α-Mn, there will be a less dense layer between the first two relatively close-packed layers. The inter-layer distance between this loosely packed second layer and the third layer is small, which explains the large δ 23 from Table 2. Omitting this less dense layer, the first two inter-layer separations for α-Mn relax by δ 13 = 0.99% and δ 34 = -1.83%, respectively. These figures are close to those obtained for Cr and Fe. Thus, according to the present ab initio calculations, the close-packed surfaces of the magnetic 3d metals show unusually small layer relaxation. We suggest that the relatively small surface relaxation obtained for some of the magnetic metals is attributed to the positive magnetic pressure around the free surfaces (see Section 3.3).

Surface energy

The calculated surface energy, stress, and excess surface stress are listed in Table 3 for some low index surfaces of 3d transition metals. For the thermodynamically stable surfaces, the present surface energies (γ 3d ) are compared to those obtained for the 4d (γ 4d ) and 5d (γ 5d ) metals in Figure 1. We observe that the three transition metal series possess rather similar surface energies both in magnitude and in trend as a function of d-occupation. This similarity is even more interesting if we realize that there are differences in the morphology of the most stable surfaces.

For instance, Cr has as the most stable surface the bcc (100) facet compared to the bcc (110) facet of Mo and W. All three bcc metals from the Cr group possess the largest surface energy across the transition series.

The local maximum (minimum) seen for Ru and Os (Tc and Re) is missing from γ 3d . However, this deviation may be ascribed to the complex crystallographic structure of α-Mn. Indeed, the surface energies calculated for the hypothetical bcc (100) facet show a shallow minimum around Mn. These results were obtained for the unrelaxed bcc (100) surface and are shown in Figure 2. We find that all three ab initio methods (PAW, EMTO and FP APW+lo) give the same surface energy trends as going from V to Fe.

For illustration of the magnetic contribution in the surface energy, Figure 2 displays also the bcc (100) surface energies obtained from nonmagnetic PAW and EMTO calculations. Clearly, there is a large drop in the surface energy of bcc Mn due to the magnetism: 1.200 J/m 2 calculated in PAW and 1.102 J/m 2 in EMTO. The magnetic effect is somewhat smaller for Cr and Fe: 0.750 J/m 2 and 0.490 J/m 2 in PAW (0.529 J/m 2 and 0.746 J/m 2 in EMTO).

Another important feature of the theoretical data from Figure 2 is that in all three calculations and for both magnetic states, the surface energy of Cr turns out to be larger than that of bcc V. This might be in apparent contradiction with previous theoretical calculations based on the linear muffin-tin orbitals method [51,52]. However, if we consider the same bcc (100) or bcc (110) crystallographic facets for V and Cr, then the four sets of theoretical results are consistent with each other. Interestingly, none of the theoretical set of results confirms the experimental estimate by Tyson and Miller [53], as it was inspired in Ref. [51]. The estimated surface energies (γ est 3d ) for V, Cr, Mn, and Fe are 2.62, 2.35, 1.54, and 2.42 J/m 2 , respectively. The deep local minimum in these "experimental" surface energies is absent from γ 3d from Figures 1 and2. Although theory might predict a weak local minimum in γ 3d near Mn when calculating for the same crystallographic facet (bcc (100), Figure 2), this theoretical minimum is much shallower than that in γ est 3d . Using continuum model approach [54], Punkkinen et al. [16] the deep minimum in the estimated surface energies for bcc 3d metals is, to a large extent, a consequence of the compressive surface stress of Cr and Mn (see Section 3.3).

Surface stress

For the thermodynamically stable surfaces, the present surface stresses (τ 3d ) are compared to those obtained for the 4d (τ 4d ) and 5d (τ 5d ) metals in Figure 3. While τ 4d and τ 5d are relatively similar, we observe that τ 3d exhibits a rather different behavior as a function of d-occupation (N d ). For N d ≤ 3 and N d ≥ 7 the present τ 3d and (τ 3d -γ 3d ) follow similar trends as those calculated for the 4d and 5d metals. However, for Cr, Mn and Fe we can see a deep minimum in τ 3d . Most surprisingly, for Cr and Mn the surface stress of the thermodynamically most stable surfaces becomes compressive: -0.32 J/m 2 for Cr and -0.22 J/m 2 for Mn. That is, the thermodynamically stable surfaces of Cr and Mn prefer a larger lattice constant than their bulk counterparts. Apparently, these two systems are the only transition metal surfaces where the Ibach's model fails.

Recently, we demonstrated that the marked deviation between τ 3d and its 4d and 5d counterparts is due to the enhanced surface magnetism in Cr, Mn, Fe, Co and Ni as compared to their bulk [16]. We define the magnetic contribution to the surface stress (τ m ) as the difference between the nonmagnetic (NM) and magnetic (Mag) values. It turns out that the magnetic surface stress for Cr, Mn, and Fe is of order of the corresponding NM surface stress. Namely, τ m Cr = 5.41 Jm -2 (102% of NM τ Cr ), τ m Mn = 5.00 Jm -2 (99%) and τ m Fe = 3.88 Jm -2 (71%). That is, surface magnetism drastically reduces the surface stress of Cr, Mn and Fe. For Co and Ni we have somewhat smaller magnetic effects, τ m Co = 1.11 Jm -2 (30%) and τ m Ni = 0.44 Jm -2 (15%).

In order to visualize the impact of magnetism on the surface stress, we realize that the decrease of the surface stress is associated with the excess magnetic pressure (P mag ) around the free surface. For 3d itinerant ferromagnetic systems, using the Stoner model and the atomic sphere approximation, we can estimate the magnetic pressure contribution relative to the non-polarized case as [55] 3P

mag V = -µ 2 I∆ d S 2 + α∆E(µ) (9)
where I, ∆ d , S stand for the Stoner exchange parameter, the width of the 3d band, and the atomic radius, respectively. The band energy change with respect to the non-polarized value, belonging to the actual value of the magnetic moment µ is denoted by ∆E(µ). The parameter α characterizes the decrease of the d-band width ∆ d with increasing volume

α = - ∂ln∆ d ∂lnS . ( 10
)
Here the symbol ∂ indicates a restricted variation during volume change (the potential is kept frozen [55]). Similarly, using these quantities we can express the energy change during spin polarization as In leading order in µ 2 we have

E mag = - 1 4 µ 2 I + ∆E(µ). (11)
∆E(µ) ∼ µ 2 4N ( F ) , (12)
so the magnetic energy becomes

E mag ∼ µ 2 4N ( F ) [1 -IN ( F )] . (13)
When the Stoner condition is satisfied, i.e. IN ( F ) > 1, E mag is always negative, as it should be. Finally, for the leading term in the pressure we obtain

3P mag V ∼ µ 2 ∆ d S 2 N ( F ) α 4∆ d S 2 -IN ( F ) . ( 14
)
For the 3d elements where magnetism occurs from Cr to Ni the values of the α parameter lies between 4 and 5, and thus the first term in the expression of 3P mag V is typically 10 (usually lies between 8 an 11). This should be compared with the values of IN ( F ), which is the largest for bcc Fe when equals to 2 and usually lies between 1 and 2. Thus we can generally say that the first term in equation ( 14) is roughly an order of magnitude larger than the second, which means that the magnetic pressure contribution for these elements is always positive. Due to the reduced coordination, the surface magnetization is generally enhanced compared to the bulk value, resulting in larger surface magnetic moment (µ s ) than the bulk counterpart (µ b ). Hence, according to equation ( 14), an excess surface magnetic pressure should appear near free surfaces. Assuming that this excess surface pressure is responsible for the magnetic surface stress, we arrive at τ m ∼ ∆µ 2 , where ∆µ 2 ≡ µ 2 s -µ 2 b represents the surface induced enhancement of the square of the magnetic moment. Punkkinen et al. [16] confirmed the almost perfectly linear relationship between τ m and ∆µ 2 . Therefore, we can conclude that the enhanced surface magnetism is the primary factor responsible for the unusually small or negative surface stress of the magnetic 3d metals. Another consequence of the positive surface magnetic pressure is the unusually small surface relaxation obtained for Cr and Fe (and to some extent also Mn).

Conclusions

We have presented a systematic theoretical study of the layer relaxation, surface energy, and surface stress of 3d transition metals, and examined the atomic number dependence of the above surface parameters. In addition, we have pointed out the unusually small surface relaxation for some of the magnetic metals, and the exceptionally low excess surface stress in the middle of the series. Finally, we have demonstrated that the revealed anomalous trends are the consequence of the enhanced surface magnetism experienced first of all in Cr, Mn and Fe.

F

. Theoretical surface energies for 3d metals. For comparison, the theoretical surface energies obtained for the 4d (squares) and 5d (triangles) series are also shown. Each data corresponds to the thermodynamically stable facets (see Table 3) of the low temperature crystallographic structures: hcp for Sc, Ti, Zn, Co, Y, Zr, Tc, Ru, La, Hf, Re and Os; bcc for V, Fe, Nb, Mo, Ta and W; B2 for Cr; α-Mn; and fcc for Ni, Cu, Rh, Pd, Ag, Ir, Pt and Au.

  .J. Punkkinen, Q.-M. Hu, S. K. Kwon, B. Johansson, J. Kollár, and L. Vitos

Figure 2 .

 2 Figure 2. Comparison between three sets of surface energies calculated for the non-relaxed bcc (001) surface of V, Cr, Mn and Fe using the VASP, EMTO and WIEN2k computer programs. Results are shown for non-magnetic (NM) and magnetic (FM) phases.

Figure 3 .

 3 Figure 3. Theoretical surface stress for 3d metals. For comparison, the theoretical surface stresses obtained for the 4d (squares) and 5d (triangles) series are also shown. Each data corresponds to the thermodynamically stable facets of the low temperature crystallographic structures (see caption for Figure2).

Table 2 .

 2 Calculated layer relaxations (in percent) for some low-index surfaces of 3d transition metals. In the case of Mn, the numbers in parenthesis represent the relaxations for the first two close-packed layers, viz. δ13 and δ34.

	May 6, 2011	4:19	Philosophical Magazine	manuscript
				Philosophical Magazine Surface Relaxation
					δ 12	δ 23	δ 34
			Sc	hcp(0001) -2.58	1.19	0.29
			Ti	hcp(0001) -7.14	2.87	-0.49
			V	bcc(100)	-12.41 0.24	2.87
				bcc(110)	-5.24	0.58	0.01
			Cr bcc(100)	-5.06	3.62	-1.04
				bcc(110)	-1.95	0.44	0.20
			Mn α(110)	-6.95 12.30	-1.83
					(0.99)	(-1.83)
			Fe	bcc(100)	-0.18	2.88	0.78
			bcc(110) F Co hcp(0001) -2.89 -0.05 Ni fcc(100) -3.74 o fcc(111) -1.31 -0.09 0.71 1.85 0.55 r Cu fcc(100) -1.61 1.48	-0.15 -0.72 0.04 0.11 1.25
			fcc(111) Zn hcp(0001)	-0.31 P 0.69 -0.65 0.04	0.31 -0.10
					e e r
						R e v i e w
						O n
						l y

Table 3 .

 3 Calculated surface energy (γ), total (τ ) and excess (τγ) surface stress (in J/m 2 ) for 3d transition metals. The results correspond to the fully layer-relaxed geometry (except bcc Mn, where the ideal unrelaxed geometry was used).

		Surface	γ	τ	(τ -γ)
	Sc	hcp(0001) 1.26 0.73	-0.53
	Ti	hcp(0001) 1.97 0.65	-1.32
	V	bcc(100)	2.40 2.12	-0.28
		bcc(110)	2.41 2.12	-0.29
	Cr bcc(100)	3.06 -0.32	-3.38
		bcc(110)	3.10 0.79	-2.31
	Mn bcc(100)	2.14 -2.24	-4.38
		α(110)	2.59 -0.22	-2.81
	Fe	bcc(100)	2.50 1.39	-1.11
		bcc(110)	2.45 1.56	-0.89
	Co hcp(0001) 2.11 2.20	0.09
	Ni	fcc(100)	2.22 1.73	-0.49
		fcc(111)	1.92 2.16	0.24
	Cu fcc(100)	1.44 2.55	1.11
		fcc(111)	1.30 1.87	0.57
	Zn hcp(0001) 0.32 1.44	1.12

  The Swedish Research Council, the Swedish Foundation for International Cooperation in Research and Higher Education (STINT), the Swedish Steel Producer's Association, the Carl Tryggers Foundation, the National Research Foundation of Korea (KSRCP-2008-08) and the Hungarian Scientific Research Fund (research projects OTKA 84078 and K77771-IN83114) are acknowledged for financial support. The calculations were performed at the Finnish IT Center for Science.

	May 6, 2011	4:19	Philosophical Magazine	manuscript
			Philosophical Magazine & Philosophical Magazine Letters
			5. Acknowledgments
			o F o F o
			r r r
			P P P
			e e r e e r e e r
				R e v i e w R e v i e w R e v i e w
				O n O n O n
				l y l y l y

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

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


Isabel Brandariz

Emilia Iglesias

Potentiometric characterization of cyclodextrin inclusion complexes of local anesthetics

Keywords: potentiometric technique, local anesthetics, inclusion complexes, cyclodextrin potentiometric technique, local anesthetics, inclusion complexes, cyclodextrin

The complexation of several local anesthetics, Xc, with  and cyclodextrins, CD, was studied, tetracaine and dibucaine complexation constants were determined, while it was found that prilocaine and lidocaine do not form inclusion complexes.

Introduction

Cyclodextrins (CD) form inclusion complexes with a great variety of guests and act as drugdelivery systems. This ability depends on the compatibility of the guest with the size of the conical cylinder cavity of the cyclodextrin and the affinity in terms of polarity with the hydrophobic wall of the cavity (1)(2)[START_REF] Dodziuk | Cyclodextrins and their Complexes[END_REF]. An interesting group of compounds that form this kind of complexes are local anesthetics that have been studied previously in our laboratory [START_REF] Iglesias | [END_REF](5)(6)(7)(8). Stability constants of procaine were obtained using a wide variety of techniques (fluorescence or UV-vis spectroscopies, electrical conductivity, potentiometry and kinetic experiments) and conditions [START_REF] Iglesias | [END_REF](5)(6)(7); in the same manner, tetracaine was studied by fluorescence spectroscopy (8) and the corresponding stability constants of the complexes were determined.

The present paper aims to compare those data with the same constants obtained by potentiometry for the case of tetracaine, with which compound both techniques can be applied, and to determine potentiometrically the stability constants of the anesthetics, dibucaine, lidocaine and prilocaine, whose constants can not be obtained by fluorescence spectroscopy, as it will be discussed below. monocation and the neutral molecule of tetracaine are able to enter into the CD cavity with the unprotonated secondary amine group end first. As the aromatic ring of the drug, the hydrophobic part of the molecule, is inside the cavity, the fluorescence emission intensity is enhanced compared to that of the molecule in water, this fact was used to determine the stability constants for the complexes of Tc and HTc + with CD. The same approach can not be used for dibucaine, because no appreciable variation in the fluorescence emission was observed in water and in the presence of CD, or with lidocaine and prilocaine, because no fluorescence was observed in the last two cases.

Taking this problem into account, a potentiometric method was considered to study the possible complexation of these anesthetics with cyclodextrins. In order to apply this technique, a small amount of sodium chloride (0.1 M) was added to all solutions to keep the ionic strength constant during the titration.

If the acid-base equilibria of an anesthetic, Xc, are represented by,

the complexation equilibrium of anesthetic cation is, then, given by,

and the complexation equilibrium of neutral anesthetic may be represented by,

The constants C K 2 and C K 3 are used in accordance with previous notation [START_REF] Iglesias | [END_REF]8) where C K 1 will represent the stability equilibrium constant of the hypothetical complex between CD and the anesthetic dication. It was also concluded that the fully protonated anesthetic molecule do not form inclusion complexes. These equations are valid for tetracaine and dibucaine. In the cases of lidocaine and prilocaine the first equilibrium, H K 1 , does not exist because there is only one amine group available to be protonated. in basic media to avoid precipitation. Under these conditions 0.01 M NaOH was used as titrant and was freshly prepared before use with water bubbled with nitrogen.

Experimental

On the other hand, the necessary amount of NaCl (Merck, p.a.) to adjust the ionic strength, I, to the desired value was added to all solutions (I=0.1 M).

Titrations were carried out in a dual-wall cell, kept at constant temperature (25.0 °C) by a thermostat-controlled system. Purified nitrogen was bubbled through the solutions in order to ensure thorough homogenization and CO 2 removal. A Crison microBu 2030 automatic buret furnished with a 2.5 mL syringe for dispensing the titrant was used. The buret was controlled with a computer that was used to read the emf values from a Crison micropH 2000 pH meter connected to two electrodes: a glass electrode (Radiometer pHG211) and a reference electrode (Radiometer REF201).

Equilibrium constants were calculated from the potentiometric titration data by means of the Hyperquad program (9). The electrode response at constant ionic strength is given by E = E 0′ + s log [H + ], where E is the emf, and the electrode parameters, the formal potential, E 0′ , and the slope, s, are required by the program to determine the constants. A separate experiment was carried out to obtain the electrode parameters: emf was measured in solutions of known proton concentration, at the desired ionic strength, s and E 0′ were obtained from a linear regression; see reference (10) Measurements in proton concentration, [H + ], allows the determination of stoichiometric or concentration constants. The values of s and E 0′ for the corresponding electrode has been used to

plot p[H]=-log[H + ].

Results

Tetracaine complexes.

Acid-base constants of TCA were previously obtained by potentiometry, while stability constants of tetracaine complexes with β-CD were determined by spectrofluorimetry; are listed in Table 1, together with the values found in this work, for purposes of comparison. When tetracaine hydrochloride was dissolved in water and titrated with strong acid, an equilibrium is established between dication and cation, Equation (1a). When β-CD is added to this system, the titration curve is displaced toward more acidic pH values due to the cation being taken by the CD, Equation (2), shifting equilibrium (1a) to the right. This effect is clearly observed in Figure 1a, where solutions of equal concentration of tetracaine hydrochloride were titrated with and without CD, with the same electrode, the same day, and with the same titrant solution. The difference between both curves of the figure is proof of the complexation with the CD, which can be used to determine the complexation constant when H pK 1 is known. Several titrations of tetracaine with cyclodextrin in acid media were introduced together in Hyperquad to yield the value listed in Table 1. It may be argued that these data are not of practical interest because too low a pH is needed to found them, but this is not true because TcH.CD is present at higher pH values, too. In fact, at slightly basic pH, Equation (1b) governs the system, and the equilibrium is established between the neutral molecule, Tc, and the cation, TcH + ; when β-CD is added to this system, it will remove both species partially to form the complexes Tc.CD and TcH.CD + , and (Table 1).

Another choice considers both constants as unknown parameters in Hyperquad and determines both values with basic titrations. Of course, the error was greater when two constants were determined instead of only one with the same set of titrations, but the results of these two procedures agree quite well within the experimental error, Table 1. In other systems (when some problems difficult the experimental determination of the constants) it may be necessary to use both kind of titrations, for example in dibucaine, as it will be discussed later.

Complexation constants of tetracaine with γ-cyclodextrin are also included in Table 1.

Titrations in acid and basic media were used to obtain these values. It is interesting to note that when tetracaine was titrated with acid in the presence of γ-CD (four titrations

HCO

, was added in a concentration of 0.033 M to adjust the pH, then it was in a much greater quantity than TCA and it could be able to compete with the anesthetic to form complexes with γ-CD. This result emphasizes the interest of using different methods to study the complex formation.

Dibucaine complexes

It was mentioned in the experimental section that the neutral molecule of dibucaine is less soluble than the other anesthetics used in this study; solutions as low as 4.5 10 -5 M were used to avoid precipitation in basic media (when CD was present 10

. So, it is possible to calculate the complexation constants of these species, since when CD is added to the solution, it removes part of these molecules to form the complexes; but if the constants are quite similar, no shift in the equilibrium is expected: the system does not detect the presence of the CD. To illustrate this point, a simulation has been performed and the result can be seen in Figure 2. The HySS program (12) was used to reproduce the expected experimental titration curves for dibucaine in the absence (circles) and in the presence of β-CD (lines from a to g) for different values of among curves are bigger (notice that pH-scale is the same in both plots). It is a good way to show the difficulties that appear when low concentrations are used, as it was stated before. Besides, this is the reason why an effect that is almost imperceptible in Figure 2a, becomes clearly observed in Figure 2b i.e., the separation of the curves, at high pH, when β-CD (curves a to g) is present and when the anesthetic is titrated alone (circles). The explanation is quite simple, the acid-base equilibrium of β-CD was taken into account in the simulations, and its influence can be noticed at the end of the titrations. It should be noted that β-CD concentration is higher than that of the anesthetic, not only in the simulation, but in the experimental titrations too, where β-CD concentration used was equal or several times higher than that of the anesthetic. As this effect was observed in the experimental data, it was considered interesting the determination of the acid-base constants of the cyclodextrins used in this work and it will be discussed later. Although there is a separation of the curves at high pH, because of the reason explained above, at lower pH, when both constants are equal, the plots coincide very well, curves a and circles Figure 2a-2b, while the other curves b-g increase their separation respect to the circles as long as the difference between , with β-CD, using drug-sensitive electrodes (14).

Lidocaine and prilocaine complexes

Lidocaine and prilocaine have only one amino group and only one acid-base equilibria, Equation (1b). Titration with base of the prilocaine hydrochloride alone, with β-CD and γ-CD, is presented in Figure 3. The concentration of the anesthetic was the same in the three experiments, and were performed with the same electrode, the same day, and with the same titrant solution. As it can be appreciated all the curves coincide quite well except at high pH, where the influence of the acid-base equilibrium of the CD is appreciable; the same behaviour can be observed in Figure 2b. At the end of the titration, prilocaine has been converted from the protonated to the neutral form, while CD begins to ionize causing the separation of the curves. Similar behaviour was observed for lidocaine. Several titrations in the presence of CD (β-CD or γ-CD) were fitted with Hyperquad by fixing the pK of CD (Table 2), the pK H of the anesthetic was obtained and the results are listed in Table 2, and it agrees with literature reports (13, 15). Titration curves and pK H values are equal in the presence and in the absence of CD, therefore it is impossible to determine the stability constant of the inclusion complexes because no separation in the curves appears. There are various possible explanations. One of them is evident and may be attributed to the absence of complex formation; on the other hand, a more subtle explanation is that complexes were formed, but the stability constants are approximately equal and no variation is expected in the curves (see Figure 2) or even the anesthetic does not penetrate very deeply in the cavity of the CD and no disturbance in the acidbase equilibrium is observed. Previously ( 16), inclusion complexes of the neutral molecules of some anesthetics, including prilocaine and lidocaine, with β-CD were determined from solubility measurements. Quite low values for the stability constants were obtained such as

1 3 26 - = M K C
for lidocaine and

1 3 96 - = M K C
for prilocain. To illustrate the expected shifts in the titration curves that would be originated by such low complexation constants, a new simulation with the HySS program (12) was performed. The results appear in Figure 3b, when the concentrations used were the same that those appearing in the experimental titrations of Figure 3a when the constants are too low, which it seems to be the case here, and, of course, this precludes the potentiometric determination of the constants.

Acid-base constants of β β β β-cyclodextrin and γ γ γ γ-cyclodextrin

As it has been discussed before, the last points of the titration curves in basic media are influenced by the presence of CD. This fact encouraged us to estimate the pK of βand γ- pK of the CD, the obtained results were 13.77 for β-CD and 13.79 for γ-CD, which agree very well with the values found using a hydrogen electrode (17). On the other hand, slightly higher values for the pK of both cyclodextrins (0.3 units more) can be found in the literature (18). This difference can be ascribed to the experimental error at high pH and besides to the calibration procedure that is quite different from ours. Previously (18) calibration of the electrode was performed with standard buffers of known pH, after that, pH was measured and the proton concentration was obtained, using the proton activity coefficient, calculated from Debye-Hückel equation with a certain value of the size parameter in this equation. In fact, it is no necessary, to have any supposition about activity or activity coefficients of individual ions to perform a calibration in concentration (19). Solutions of known concentration can be used to obtain a calibration line without any kind of supposition, that could affect the obtained values in an indeterminate way.

Discussion

Cyclodextrins are truncated cone-shaped molecules with a hollow, tapered cavity of 7.9-Å depth. The top and bottom diameters of the cavity are 6.8 and 7.8 Å for β-cyclodextrin, and 8.3 and 9.5 Å for γ-cyclodextrin, respectively (

As it was stated before, the ability to form inclusion complexes depends on the compatibility of the guest molecule with the size of the cavity and the affinity in terms of polarity with the hydrophobic wall of the cavity. To explore the possible inclusion of the anesthetic molecules into the CD, the size of these compounds was optimized using PM6 method by means of MOPAC (20), (identical results were found with PM3). It is expected that the molecule will enter with its more hydrophobic part oriented towards the cavity, that is, the aromatic ring. Some relevant distances are depicted in scheme 1. It can be shown that tetracaine has a size perfectly compatible with the cavity of CDs and give rise to the highest constants. On the other hand, lidocaine is the molecule with the biggest diameter because of the methyl groups in the benzene ring, which fact difficults the complexation with β-CD, but not with γ-CD that has the biggest cavity. Therefore, the absence or the small amount of complex formation may be ascribed to a lack of affinity with the hydrophobic cavity of the CD. The same seems to be true for prilocaine, that has a smaller size than lidocaine.

Both compounds exhibit a much greater solubility in water than dibucaine or tetracaine, and it seems they prefer the polar environment of water instead of the hydrophobic cavity of the CDs. On the contrary, dibucaine with a large aromatic part, even relatively big, can form complexes in significant amount. A detailed study of the geometry of these molecules, neutral, protonated and in different solvents can be found in ref (21). The use of the potentiometric technique to determine inclusion constants of several local anesthetic with CDs has been discussed and the complexation constants determined. It was found that tetracaine, neutral or protonated, forms complexes in significant amount with β and γcyclodextrins. This system had been studied previously with similar results, in this laboratory, using fluorescence. Although this is a very sensible technique, there are other anesthetics that can not be studied with this method, i.e. no fluorescence or change of fluorescence is observed in the presence of CD, so potentiometry was employed. This is the case of dibucaine that forms complexes with both CDs; also prilocaine and lidocaine, which do not form complexes or do it in a small amount that can not be detected with the potentiometric technique.

Conclusions

  Tetracaine [TCA, 2-(Dimethylamino)ethyl 4-(butylamino)benzoate] has two amino groups, therefore, depending on the acidity of the medium it may exit as a neutral molecule, a monocation, or a dication. The same is true for dibucaine [DCA, 2-Butoxy-N-[2-(diethylamino)ethyl]N-(2-Methylphenyl)-2-(propylamino)propanamide] only one amino moiety is present and there is just one acid-base equilibrium between neutral and monoprotonated molecules, see scheme 1. It was that (8) the

  High-purity TCA.HCl, DCA.HCl, LDA.HCl, PDA.HCl and β-cyclodextrin were purchased from Sigma-Aldrich (St Lous, MO, USA), while γ-cyclodextrin were purchased from Cyclolab R & D Laboratory Ltd., Budapest, Hungary.

2 .

 2 An aqueous aliquot containing the anesthetic hydrochloride (TCA, DCA) was titrated with a standard solution of hydrochloric acid (Merck p.a.) to determine H pK 1 , and with a standard solution of sodium hydroxide (MERCK p.a.) in order to determine H pK When cyclodextrin was added to the solution of the anesthetic hydrochloride, the titration with HCl yields case of LDA.HCl and PDA.HCl only titrations with base (with and without CD) were performed because there is only one acid-base equilibrium. The anesthetic concentration was made varied from M was higher in acid medium than in basic medium, because in the later medium the anesthetic is in the form of a neutral molecule that easily precipitates. This is especially true for dibucaine, whose concentration was kept as low as M

  explain the system behaviour. In Figure1b, TCA.HCl is titrated with strong base, with and without β-CD, the difference between the curves is caused by the formation of the complexes, and this kind of titration can be used to determine their stability constants. At this point two options are possible, acid titrations, this value can be fixed in the Hyperquad program to obtain the other complexation constant,

  Figure 2a were: [DCA]=1.5 10 -4 M and [β-CD]=1.5 10 -3 M titrated with [NaOH]=0.01 M. In Figure 2b these values were increased 10 times to illustrate the influence of concentration, i.e.

  for the case of prilocaine, i.e. 40 ml of [PCA]=1.0 10 -3 M and [β-CD]=4.9 10 -3 M titrated with [NaOH]=0.02 M, while the constants the following values (a)1.98, (b)1.88, …..(g)1.58. The new fact in this plot, compared to the other simulations, is that the curves are very close one to another, showing the difficulties that appear

Figure 1 .Figure 2a .

 12a Figure 1. Tetracaine hydrochloride titration without (circles) and with β-CD (squares) a) in acid media, [TCA]=3 mM, [β-CD]=3 mM, [HCl]=0.1, v 0 =40 mL b) in basic media, [TCA]=1 mM, [β-CD]=1 mM, [NaOH]=0.010 M, v 0 =40 mL (v 0 is the initial volume of titrand)

  Scheme 1 161x229mm (300 x 300 DPI)

  The comparison of the data obtained in this paper with those determined previously from fluorescence measurements(8), shows that the stability constants with β-CD agree very well, while C K 3 with γ-CD, obtained in this work, is significantly greater. The explanation may reside in the fact that concentration of TCA in fluorescence experiments was quite low ~10 µM (it was selected taking into account that the absorbance should be ~0.2).

									But buffer,	3 -	3 /CO -2
			F o r
								P
									e
									e r
									R
									e
									v i e
									w
									O n l
									y
									), Hyperquad fitting
	yields	log 2 = C K	40 . 2	and	pK	1 = H	26 . 2	, therefore
	value of	1 =	. 1	96

H pK 1 can be determined simultaneously with C K 2 in the presence of γ-CD; as it can be seen this result agrees perfectly with 24 .

2 1 = H pK

on Table

1

obtained without γ-CD. If the complexation equilibrium is omitted from the system, an incorrect H pK was determined, which clearly shows that the presence of γ-CD affects the acid-base equilibrium (of course, similar behaviour is observed with β-CD).

Table 1 .

 1 On the other hand, although the sizes of the cavities of the CDs are adequate to include the molecules of the anesthetics used in this study (with the exception of lidocaine in β-CD); the absence of complexation can be ascribed to different affinity of guest the molecule and the cavity. Acid-base and complexation constants of tetracaine and dibucaine at 25ºC and 0.1 M ionic strength of NaCl

	Aknowledgments
	Financial support from the Direccion General de Investigacion (Ministerio de Educacion y Ciencia)
	of Spain and FEDER (Project CTQ2005-07428/BQU) and from Direccion General de Programas y
	Transferencia de Conocimiento (Ministerio de Ciencia e Innovacion) of Spain (Project CTQ2008-
	04429/BQU)is gratefully acknowledged.

*

ref

(8)

Table 2 .

 2 Acid-base constants of lidocaine, prilocaine, β-cyclodextrin and γ-cyclodextrin at 25ºC and 0.1 M ionic strength of NaCl Structures of the anesthetics tetracaine, dibucaine, prilocaine and lidocaine.

		pK	H
	Lidocaine	7.92±0.01
	Lidocaine with β β β β-cyclodextrin	7.89±0.01
	Lidocaine with γ γ γ γ-cyclodextrin	7.90±0.01
	Prilocaine	8.01±0.01
	Prilocaine with β β β β-cyclodextrin	8.03±0.01
	Prilocaine with γ γ γ γ-cyclodextrin	7.99±0.01
	β β β β-cyclodextrin	11.95±0.01
	γ γ γ γ-cyclodextrin	11.88±0.02

URL: http:/mc.manuscriptcentral.com/tandf/gsch Email: suprachem@mail.cm.utexas.edu

URL: http:/mc.manuscriptcentral.com/tandf/gsch Email: suprachem@mail.cm.utexas.edu Supramolecular Chemistry