Shadyar Farhangfar
email: shadyar.farhangfar@gmail.com

Size-Dependent Thermoelectricity in Nanowires

Keywords: numbers: 73.50.Lw, 73.23.-b, 73.63.Nm

By employing a semiclassical Boltzmann transport equation and using an energy and size dependent relaxation time (τ ) approximation (RTA), we evaluate the thermoelectric figure-of-merit ZT of a quantum wire with rectangular cross-section. The inferred ZT shows abrupt enhancement in comparison to its counterparts in bulk systems. Still, the estimated ZT for the representative Bi 2 Te 3 nanowires and its dependence on wire parameters deviate considerably from those predicted by the existing RTA models with a constant τ . In addition, we address contribution of the higher energy subbands to the transport phenomena, the effect of chemical potential tuning on ZT , and correlation of ZT with quantum size effects (QSEs). Particularly, we show that within each subband, ZT has an optimum value which depends on wire dimensions and chemical potential. The obtained results are of general validity for a wide class of systems and may prove useful in the ongoing development of the modern thermoelectric applications.

Introduction

Thermoelectricity deals with conversion of heat to electricity and vice versa. Last few years, mainly due to the recent advances in material science and fabrication techniques on the one hand and a quest for alternative sources of energy generation on the other hand, have witnessed a rejuvenation of interest in the thermoelectric phenomena [1,2,3,4,5,6,7,8,9,10,11,12,13]. Performance of a thermoelectric device is defined through its dimensionless figure-of-merit ZT ≡ S 2 σT /κ. Here, T is the absolute temperature, S stands for the thermopower (the Seebeck coefficient), and σ and κ represent the electrical and thermal conductivities, respectively. At first glance, a large ZT might be achievable by maximizing σ and S, and by minimizing κ. The problem, however, lies in the fact that these quantities are interdependent and each of them is rather sensitive to the material properties, to the temperature, and to the size and geometry of the underlying building blocks in a miniaturized thermoelectric module. Consequently, often one of the two not-completely-independent strategies is exploited: increasing of the power factor P ≡ S 2 σ which is largely an electronic property, or decreasing of κ which consists of electronic κ e and lattice κ ph. contributions by manipulation of the latter. A large P can be achieved by tuning of the chemical potential level through doping or external electric fields (gating) and by having control over the scattering mechanisms of the charge carriers in the device. Reduction in κ ph. , instead, can be obtained by exploitation of material systems built up from heavy elements and by intentional enhancement of the phonon scattering events through introduction of lattice imperfections (impurity atoms, superlattice structures, etc.) or by intensifying of the surface scattering of the heat carriers from the specimen boundaries in the systems with reduced dimensionality.

In what follows, assuming a constant κ ph. , we focus on the charge carrier contributions to the thermoelectric properties of a one-dimensional (1D) wire with discrete energy levels. To demonstrate the applicability of the present approach, we evaluate ZT in bismuth telluride nanowires and discuss the obtained results with an eye on the predictions of the earlier models.

Theoretical model

Current generic models (as compared to the ab initio calculations) for evaluation of the thermoelectric figure-of-merit are mostly based on the solution of Boltzmann transport equation (BTE) with various degrees of sophistications. An efficient approximation for solving BTE can be attained by its linearization and through introduction of a relaxation time τ , the time period within which the system gains its equilibrium after removal of the external stimulus. Consequently, the generalized transport coefficients of the system L (l) can be obtained through the relation [START_REF] See | Electrons and Phonons[END_REF]

L (l) ≡ νq 2-l 2π - ∂f 0 ∂E τ (k)(E(k) -ζ) l v(k)v(k)dk, (1)
where, q is the electrical carrier charge, E(k) is the energy-wavevector dispersion relation, hv(k) ≡ ∇ k E(k) defines the velocity operator (h is the reduced Planck constant), and ν is the valley degeneracy. The equilibrium Fermi-Dirac distribution function is given by f

0 (E)≡[1+exp β(E -ζ)] -1 , where β ≡ 1/k B T (k B is the Boltzmann constant)
and ζ denotes the chemical potential. Correspondingly, the measurable transport quantities are defined by the following expressions:

σ ≡L (0) , S ≡ 1 T L (1) L (0) -1 , κ e ≡ 1 T L (2) -L (1) L (0) -1 L (1) . (2)
Here, σ designates the electrical conductivity tensor, S is the thermopower (Seebeck coefficient), and κ e is the electronic contribution to the thermal conductivity [START_REF] See | Electrons and Phonons[END_REF].

The precise form of the dispersion relation above can only be derived through detailed band structure calculations and depends strongly on the material properties and on the boundary conditions in the system of interest. Similarly, the relaxation time for any specific process depends on the density of energy states (DOS) and on the types of the scattering mechanisms taken into account. In most calculations, however, one expresses DOS via the single-particle eigenenergies and assumes a power-law dependence of the relaxation time on energy, τ (E) ∝ E α . It was only under such conditions that A. F. Ioffe arrived at his compact expressions in terms of the Fermi-Dirac integrals for the evaluation of the transport coefficients in Eq. ( 1) above [START_REF] Ioffe | Semiconductor Thermoelements and Thermoelectric Cooling[END_REF][START_REF] Nolas | Thermoelectrics, Basic Principles and New Materials Developments[END_REF]. Despite this, such basic assumptions have been widely overlooked in the latter evaluations of the ZT values in low-dimensional systems, leading ultimately to the estimations of unrealistically high, monotonously size-dependent, ZT 's in atomically thin nanowires [17,18,19,[START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF]. It is worth reminding that in the limit of ultra narrow nanowires, the whole concept of a size-independent DOS and any argument subsequent to it will collapse. For example, predicting an extremely high value of ZT ∼ 40 (based on evaluation of transport properties with a size-independent relaxation time and DOS) for nanowires as thin as ∼1 Å, where there is only a very sparse, discontinuous and strongly size-dependent density of states, will lead to misleading conclusions. (We are not going here, by any means, to underestimate the value of the pioneering work of Dresselhaus and coworkers [17,18,19,[START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF], rather to point out the shortcomings of their approach explicitly, and to take their widely cited model one step further.) In such cases, one has to resort to the other proper techniques like, e.g., tight-binding models and check for the consistency of the obtained results. Furthermore, as shown elsewhere [START_REF] Mahan | [END_REF], even minor variations in the form of DOS [and subsequent changes in τ (E)] will affect ZT values dramatically. It is also elucidating to remind that, based on a purely mathematical analysis of the functional in Eq. (1), Mahan and Sofo estimate in [START_REF] Mahan | [END_REF] a universal (that is, independent of size, temperature, and specific material properties) finite upper limit for ZT in semiconductor thermoelectric materials. Based on these facts and to circumvent the above-mentioned crucial shortcomings in the existing theoretical models for low-dimensional thermoelectric systems, starting with the basic principles, below we obtain the transport coefficients and ZT values specifically for a 1D nanowire with quantized energy levels. The key issue in such a treatment is the derivation of a proper, i.e. size and energy dependent, expression for the relaxation time.

To this end, neglecting the many many-body electronic correlations and localization effects, we first assume that the single-particle charge carriers are confined in a wire aligned along the z-axis. The wire has dimensions w, t, and L (width, thickness, and length, respectively) and it's volume is Ω≡wtL. The wave functions and eigenenergies of such a particle can be derived easily by solving of the corresponding Schrödinger equation and are given by

Ψ ij (k) = 2 √ Ω sin iπx w sin jπy t exp (ikz) (3)
and

E ij (k)≡E ij + h2 k 2 2m z , (4)
where

E ij ≡ (h 2 π 2 /2)(i 2 /w 2 m x + j 2 /t 2 m y
) are the subband energies, i, j ∈ N, and k is the wavenumber of the traveling particle wave along the wire axis. Above, m x , m y , and m z are the components of the carrier effective mass along the wire axes.

Let us now assume that N scatterers each with a scattering strength V 0 are randomly distributed at positions R j along the wire,

V (r) = N j=1 V 0 δ (r -R j ) . (5)
Above, δ(r) is the three-dimensional Dirac delta function and r ≡ (x, y, z) indicates the position vector. The corresponding scattering rates can now be obtained through the Fermi's golden rule

τ -1 mnk→m ′ n ′ k ′ = 2π h | Ψ mn (k) | V (r) | Ψ m ′ n ′ (k ′ ) | 2 . (6)
Here, | Ψ mn (k) and | Ψ m ′ n ′ (k ′ ) stand for the initial and final states of the single-particle wave functions, respectively. Assuming now an elastic scattering with conservation of momentum, the transition rate between the states (m, n) and (m ′ , n ′ ) can be evaluated by averaging over the configuration of the scattering centers [22]; the result is:

τ -1 mn→m ′ n ′ = π 2 ̺ hV 2 0 Λ mn m ′ n ′ g(E)| E=E m ′ n ′ (k ′ ) , (7)
where

Λ mn m ′ n ′ ≡ (2 + δ mm ′ )(2 + δ nn ′ ), ̺ ≡ N/Ω
is the volume density of scatterers, and the energy density of states g(E) is given by

g(E) = s 2πh Ω wt m z 2 i,j Θ (E -E ij ) E -E ij , ( 8
)
where s is the spin degeneracy, and Θ designates the Heaviside step function. In Eq. ( 7) above, as in most cases the density and strength of scatterers are unknown and may vary from one individual nanowire to another, it will be helpful to make use of the macroscopic definition of carrier mobility µ ≡ |q| τ /m * (|q| is the magnitude of carrier charge and m * is the effective mass; for a 1D wire aligned along the z-axis, m * = m z ) and define analogously a mean value µ ij for the carrier mobility within each subband (i, j) through the relation µ ij ≡ |q| τ ij /m * . Substituting for the relaxation time in Eq. ( 7) above, the quantity ̺V 2 0 for the subband (i, j) can now be accounted for as 1

̺V 2 0 | ij = sν 4 Ω wt µ ij m z |q| h2 m z 2k B T m,n 1 τ mn . (9)
Here, the lifetime of the state (m, n) is defined through the Matthiessen's rule τ mn

-1 ≡ m ′ n ′ τ mn→m ′ n ′ -1
, while the corresponding expectation values are given by

τ mn→m ′ n ′ ≡ τ mn→m ′ n ′ (ε)εg(ε)∂ ε f dε εg(ε)∂ ε f dε , (10)
where, ε ≡ βE stands for the reduced energy. At this point, it should be emphasized that substitution for ̺V 2 0 in Eq. ( 9) above is not a principal restriction on the present model. If the density of scatterers ̺ and their strength V 0 (or only the product ̺V 2 0 ) are known for any single nanowire, one can evaluate all the kinetic coefficients accordingly. The use of a constant mobility here is only to facilitate the comparison to the experimental data and to the results of previous models with a constant relaxation time and a fixed mobility. Correspondingly, ̺V 2 0 can be used as a fitting parameter to experimental data, and will thus provide useful information on the presence and distribution of scattering centers in a fabricated specimen.

Next, substituting for the corresponding expressions and making a coordinate transformation where ε and the reduced chemical potential ζ * ≡ βζ (for electrons) are measured from the bottom of the conduction band [START_REF] See | Electrons and Phonons[END_REF], the transport coefficients in Eq. ( 1) can be expressed as

L (l) = sν wt 4πµ β l+0.5 m z 2 q (1-l) h m,n 1 τ mn i,j F l ij ≡ i,j L (l) ij , (11)
where

F l ij ≡ F l ij (ζ * ) and
F l ij ≡ ∞ 0 cosh -2 ε + ζ * 2 τ ij (ε)(ε + ζ * ) l ε -ε ij dε. ( 12
)
In Eq. ( 11) above, L

ij can be interpreted as the contribution of a single subband (i, j) to the generalized transport coefficient L (l) . (The contribution of the holes can be accounted for analogously.) It is also noticeable that the expression in Eq. ( 12) is reminiscent of the Fermi integral (see Eq. ( 14) in the Appendix) appearing in the corresponding models with a size-independent relaxation time introduced originally in [START_REF] Ioffe | Semiconductor Thermoelements and Thermoelectric Cooling[END_REF] and exploited later by Dresselhaus and coworkers in their primary [18] and subsequent studies [19,[START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF].

Results and discussion

Next, we evaluate ZT for a bismuth telluride nanowire at T = 70 K and at 300 K. To make the comparison between the predictions of the present model and those of the former studies more transparent, we solve the equations under exactly the same assumptions and use the same material parameters, m x = 0.32, m y = 0.08, and m * = m z = 0.02 for the electronic effective masses (all in units of the free electron mass), and take the same value for the subband mobilities of electrons along the wire axis as adopted in previous calculations [18,19,[START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF], µ ij ≈ µ = 1200 cm 2 V -1 s -1 . In both the models addressed here (the present one and that with a constant-τ in Ref. [18,19,[START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF]), only electronic contributions to the transport properties are taken into account. The lattice thermal conductivity of bulk Bi 2 Te 3 is κ bulk ph. ≈ 1.5 W/Km. Indeed, generally, thermal conductivity comprises of contribution of phonons with a wide range of wavelengths. In a more precise picture, particularly at lower dimensions, one has to deal with such a spectrum of phonon wavelenghts rather than with the simplified concept of a single phonon mean free path (and a single value for κ). For example, in the case of silicon, phonon mean free paths can span several orders of magnitude at room temperature and can have significant contribution at wavelenghts as large as 1 µm [23]. However, here, as in the corresponding studies with a constant relaxation time [18,19,[START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF], it is assumed that the phonon confinement effects become considerable only if the lateral dimensions of the wire, w and t, are comparable in size to the phonon mean free path λ ph. . As for bismuth telluride nanowires λ ph. ∼ 1 nm [18], for wires with w, t > ∼ λ ph. , one can safely assume κ ph. ≈ κ bulk ph. [START_REF]This is a conservative estimation. If one takes the phonon confinement effects into account[END_REF]. Figure 1 shows dependence of the electronic contribution to the figure-of-merit of a wire with square cross-section on its size w and on ζ * . Here, in contrast to the constant-τ model of Hicks and Dresselhaus (HD) (see, the Appendix), dependence of ZT on wire size is non-monotonous and after reaching its optimum at a certain width and chemical potential, falls starkly to zero. This is a direct consequence of disappearance of electronic density of states for ultra thin wires. The most distinct differences between the two models are as follow: the peak values of the σ, κ e , and ZT in the present model show a shift toward larger wire thicknesses, ∼ 10 nm, as compared to the plain w -2 -dependence of the corresponding quantities predicted by the HD model; the thermopower S given here is temperature dependent and reaches the overall high values of the HD model only in very narrow wires; in contrast to the HD model, here the position of optimum ZT in ζw-plane has a visible dependence on temperature and, as temperature falls down, shifts toward higher thicknesses; the asymptotic values of the transport coefficients for a model square wire with w = 50 nm given by the HD model are unrealistically smaller than those obtained here (see Table 1 in Appendix).

To illustrate the effect of subband contributions and that of QSEs, the calculations were repeated for the first three subbands, up to (i, j) = (1, 3) (see, Fig. 2). Here, especially at T = 70 K, one can distinguish two clear maxima for σ, κ e , and S which arise from QSEs [22,[START_REF] Shadyar | [END_REF] and which are smeared out in ZT . This observation is rather general for other subbands too and is a consequence of the fact that σ and κ e have similar functional dependencies on the wire size and the QSEs distinguishable independently in each of them (and to a lesser extent in S) are compensated in ZT by one another. Another interesting observation is the abrupt enhancement of κ e at thinner wire diameters and in the region where ζ * approaches its optimum. Here, noticeably, κ e ≫ κ ph. . Alike, the HD model predicts an even stronger enhancement of κ e in thinner wires (Appendix). Dependence of the thermal conductivity on dopant concentration (cf. tuning of ζ), leading ultimately to the takeover of κ ph. by the electronic contribution κ e , was observed in experiments with bismuth-antimony alloys some years ago [26].

Summary

The approach described here, which is based on the formulation of the scattering rates of the charge carriers in terms of an energy and size dependent expression for the density of states, can be readily extended to address transport and thermoelectric properties of a wide range of quantum systems subject to more sophisticated forms of the relevant DOS's. The model also holds potential to account for the non-diagonal contributions of the transport tensor elements to the kinetic properties. It can equally be exploited in the metallic or semiconducting regimes and it considers inherently the important issue of quantum size effects in low-dimensional structures.

Appendix

Table 1 outlines contributions of the higher energy subbands to the transport properties and to ZT in a 50-nm-thick Bi 2 Te 3 nanowire with effective mass components m x = 0.32, m y = 0.08, and m z = 0.02 along the wire axis (all in units of the free electron mass). The calculations were done for a relatively thick wire for two main reasons. First, at around 50 nm, in both the models (the HD model and the one presented here) the transport coefficients approach their asymptotic values and, secondly, fabrication of such wires is experimentally easier. Below, we also summarize the main results of the constant-τ model and depict its predicted values of the transport properties as a function of wire thickness in Fig. 3. A comparison of this figure with Figs. 1 and2 in the text, should further illustrate the enhanced differences between predictions of the two models in thinner nanowires.

Constant relaxation-time model

Derivation of the transport properties in bulk systems, based on a constant relaxation-time approximation, was originally given in [START_REF] Ioffe | Semiconductor Thermoelements and Thermoelectric Cooling[END_REF] and later discussed in detail by Nolas et al. in [START_REF] Nolas | Thermoelectrics, Basic Principles and New Materials Developments[END_REF]. Later, extension of these results to the case of one-dimensional nanowires was obtained by Dresselhaus and coworkers [18,19]. A more recent account is given in [START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF]. Below, we summarize the results. correspond to the subband numbers up to which the electronic (vs. holes) contributions are taken into account (κ ph. = 1.5 W/Km). HD represents the one-band constantτ model [18,19,[START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF]. Note that with valley degeneracy ν = 6, σ and κ e become simply six times larger while S remains intact. Dependence of ZT on ν is of the form ZT ∝ (κ e + κ ph. /ν) -1 and its values for ν = 1 and ν = 6 are given explicitly.

The electronic (vs. holes) contribution to the transport coefficients σ, S, and κ e are defined through the generalized transport coefficients in Eq. ( 1). Substituting for a constant relaxation time in terms of a fixed mobility as assumed in [17,18,19,[START_REF] Dresselhaus | Thermoelectrics Handbook: Micro to Nano[END_REF], µ = |q| τ /m * , one can obtain the following expressions:

L (0) = 1 2 D e F -1/2 , L (1) = (k B T ) D e q 3 2 F 1/2 - 1 2 ζ * F -1/2 , L (2) = (k B T ) 2 D e q 2 5 2 F 3/2 -3ζ * F 1/2 + 1 2 ζ * 2 F -1/2 ,
where

D e ≡ ν 2 |q| πw 2 2k B T h2 1/2 √ m * µ, (13)
and

F i ≡F i (ζ * ) ≡ ∞ 0 ε i dε e (ε-ζ * ) + 1 . (14)
Using the definitions in Eq. ( 2), the above expressions can be further simplified to:

σ = ν|q| πw 2 2k B T h2 1/2 √ m * µF -1/2 , S = k B q 3F 1/2 F -1/2 -ζ * , κ e = 1 |q| 2ν πw 2 2k B T h2 1/2 k 2 B T √ m * µ 5 2 F 3/2 - 9F 2 1/2 2F -1/2 .
The dimensionless figure-of-merit, ZT = S 2 σT /(κ e + κ ph. ), can now be written as:

ZT = 1 2 3F 1/2 F -1/2 -ζ * 2 F -1/2 1 B + 5 2 F 3/2 - 9F 2 1/2 2F -1/2 (15)
with

B ≡ 2ν πw 2 2k B T h2 1/2 k 2 B T √ m * µ |q| κ ph. . (16)
The most distinct difference of the above expressions with those presented in the Text is lack of oscillatory dependence of the transport coefficients on the wire size and on the quantization of the energy levels, i.e., absence of the quantum size effects. Also, it is noticeable that here the Seebeck coefficient S is wholly independent of size and the temperature, and σ and κ e have only a bare monotonic dependence of the form ∝ 1/w 2 on the wire width. Such dependencies lead ultimately to unrealistically large values for σ, κ e and, eventually, ZT in ultra thin wires. Figure 3 presents the obtained results for a Bi 2 Te 3 nanowire with the same parameters as given in the Text. This work is dedicated to the memory of Amirkhan Qezelli, my uncle and childhood friend.

Figure 1 .

 1 Figure 1. (a) Dependence of the transport properties and the thermoelectric figureof-merit ZT (only the electronic contribution) on wire width w and on the chemical potential ζ (measured in units of thermal energy k B T ) for a square Bi 2 Te 3 nanowire at T = 300 K. The valley degeneracy is ν = 6. Only the first subband contribution is taken into account. (b) T = 70 K. Notice the many fold decreases of the peak values and their shift toward larger wire thicknesses.

Figure 2 .

 2 Figure 2. (a) Dependence of the transport coefficients on the wire size and on the chemical potential at T = 300 K. Same parameters as in Fig. 1, except that here (m, n) = (1, 3). The clearly visible bumps in σ and κ e are manifestations of the quantum size effects. (b) T = 70 K.

T

  (K) (m, n) σ(10 6 /Ωm) κ e (W/Km) S(µV/K) ZT (ν = 1) ZT (ν = 6) 70 HD 1.1 × 10 -5 8.

2

 2 Te 3 nanowire (aligned along the z-axis) with w = 50 nm at a fixed chemical potential ζ * = -3 (measured from the bottom of the conduction band). The indices (m, n)

Figure 3 .

 3 Figure 3. (a) Thermoelectric figure-of-merit ZT for a square Bi 2 Te 3 nanowire at T = 300 K obtained by the constant-τ model. The chemical potential ζ is measured in units of thermal energy k B T . The effective mass values are the same as in the Text (m x = 0.32, m y = 0.08, and m z = 0.02 along the wire axis). (b) T = 70 K. Notice that the predicted enhancement of the transport properties here emerge only in ultra thin wires (w < ∼ 1 nm). This has to be compared to the predictions of our model presented in the Text.

Table 1 .

 1 Electrical conductivity σ, electronic thermal conductivity κ e , the Seebeck coefficient S, and the thermoelectric figure-of-merit ZT for a quadratic Bi

	300 70	HD (1,1)	1 × 10 -6 2.2 × 10 -5 7.2 × 10 -5 0.041 0.052	-401.3 -401.3 -92.9	8.1 × 10 -5 4.9 × 10 -4 7.2 × 10 -4 4.3 × 10 -3 0.016 0.082
	300	(1,1)	0.089	0.494	-86.0	0.099	0.265
	70	(1,2)	0.144	0.182	-103.9	0.065	0.252
	300	(1,2)	0.348	1.92	-88.5	0.240	0.377
	70	(1,3)	0.242	0.350	-105.5	0.102	0.314
	300	(1,3)	0.758	4.13	-92.3	0.344	0.442
	70	(2,1)	0.160	0.206	-95.9	0.061	0.228
	300	(2,1)	0.354	1.96	-86.6	0.230	0.360
	70	(2,2)	0.552	0.691	-107.8	0.205	0.477
	300	(2,2)	1.39	7.62	-89.1	0.362	0.420
	70	(2,3)	0.908	1.31	-109.7	0.273	0.492
	300	(2,3)	3.02	16.41	-93.0	0.438	0.471
	70	(3,1)	0.345	0.428	-100.4	0.126	0.359
	300	(3,1)	0.791	4.37	-87.5	0.309	0.393
	70	(3,2)	1.17	1.46	-112.9	0.353	0.611
	300	(3,2)	3.10	16.97	-90.0	0.408	0.438
	70	(3,3)	1.86	2.67	-114.2	0.406	0.578
	300	(3,3)	6.74	36.48	-94.1	0.471	0.487


Árpád Barna

László Kotis

János Lábár

Attila Sulyok

Attila L Tóth

Miklós Menyhárd

Peter Panjan

Growing imbedded Ni 3 C rich layer with sharp interfaces by means of ion beam mixing of C/Ni layers

C/Ni bilayers of various layer thicknesses (20-40 nm) were ion bombarded using Ga + and Ni + projectiles with energies in the range of 20-30 keV. Ion bombardment resulted in the growth of a Ni 3 C rich layer with the following features: a. / sharp carbon / Ni 3 C rich layer interface b. / the amount of Ni 3 C produced by the irradiation proportional to the square root of the fluence and depended on the type of projectile c. / good correlation between the distribution of vacancies produced by the ion bombardment and the distribution of the Ni 3 C. The formation of the metastable Ni 3 C compound was explained by a vacancy assisted process. The sharp interface is the consequence of a relaxation process removing the intermixed Ni from the carbon layer. The square root of fluence dependence of the thickness of the Ni 3 C rich layer can be explained by a usual diffusion equation considering moving boundaries.

INTRODUCTION

There is a constant need for materials with desired physical and chemical behavior. This need can be, as least partially, fulfilled by non-equilibrium phases, which might exhibit unique, tailored features. Similarly, nano structures also display exceptional physical and chemical features. There are various methods which can be applied to produce non-equilibrium structures; one of them is ion bombardment, which can also be used for producing such structures in nano regions.

Ion irradiation of solid prompts a series of processes, like collisional events, large temperature increase in small volumes, radiation enhanced diffusion, chemically driven transport, etc, resulting in serious changes of the irradiated solid. One of the processes occurring during ion irradiation of a solid is ion beam mixing (IBM), which might cause the broadening of the originally abrupt interface between two pure materials by prompting the transport of the atoms from one pure region to the other [1][2][3]. IBM has been applied for producing tailored materials from the eighties [START_REF]Beam-Solid Interactions; Physical Phenomena[END_REF]. Roughly speaking, the IBM process can be divided into two parts as a./ the mixing taking place due to collisional cascade, thermal spike etc initiated directly by the projectile b./ the relaxation of the highly non-equilibrium damaged region. The relaxation in many cases is terminated by a fast cooling period, consequently the damaged region generally does not reach equilibrium condition, rather it remains in a metastable state.

A huge variety of metastable states can be observed ranging from miscible amorphous state [START_REF] Milosavljevic | [END_REF] to the immiscible amorphous states [6,7]. The observed kinetics vary greatly too, the mass transport in limiting cases might be symmetric or asymmetric [8,9]. Applying focused ion beam (FIB) nano patterning of ion beam affected material is possible. Thus using ion bombardment, we can produce nano regions of non equilibrium material of the desired properties.

that until reaching a critical fluence, the carbon/mixed layer (containing Ni 3 C) interface remains sharp (despite the large carbon transport occurring through the interface). The amount of produced Ni 3 C depends on the square root of the fluence and the type of projectile besides the ion energy. We will also show that the alterations produced by the ion bombardment can only be explained by considering the relaxation processes taking place after primary mixing is terminated.

EXPERIMENTAL

C/Ni multilayered specimens have been prepared by sequential sputter deposition of pure carbon and nickel on silicon (111) substrate. The nominal structures of the specimens were: specimen A; 3x(C 28 nm /Ni 40 nm)/Si substrate, specimen B; 3x(C 21 nm /Ni 37 nm)/Si substrate.

The actual structures of the specimens before and after irradiation were determined by cross sectional transmission electron microscopy (XTEM). The samples for XTEM studies were prepared by applying the usual procedure [START_REF] Barna | Sample Preparation Techniques for TEM[END_REF]; first mechanical polishing was applied, which was followed by Ar + -ion milling. XTEM images were taken by a Philips CM20 200 kV analytical microscope.

The samples were irradiated by applying various energies and fluences of Ga + and Ni + ions.

A 200x200 µm 2 area of the specimens were irradiated by Ga + in a LEO 1540XB (FEG SEM -FIB) cross beam system, using the Canion FIB optics. The area to be irradiated had been located using scanning electron microscopy (SEM) secondary electron image (SEI), and then checked on a single pass FIB SEI image. The angle of incidence of the irradiation was 5 o (with respect the surface normal). The energies of the Ga + projectiles were 20 and 30 keV, and the ion current density varied in the range of 3-20 nA/mm 2 . The applied fluences, 0, 25, 50, 100, 200, 400 and 800 Ga + /nm 2 , were determined by the time of the single pass irradiation.

10x10 mm sample was irradiated by Ni + in Rossendorf [START_REF]AIM -Center for Application of Ion Beams to Materials Research[END_REF]. The angle of incidence of the irradiation was 0 o , and energies of 20 and 30 keV were applied.

Auger electron spectroscopy (AES) depth profiling was carried out in our dedicated systems by a STAIB DESA 100 pre-retarded CMA in direct current mode. The primary electron current was 30 nA with a diameter of about 60 µm.1 keV Ar + ions were used for AES depth profiling with angle of incidence of 74 o (with respect to the surface normal). The ion current was kept constant during sputtering. The sample was rotated (4 rev/ min) during ion bombardment. These parameters were chosen to minimize the ion bombardment induced surface and interface morphology development.

Deriving elemental concentrations from the Auger peak intensities

The shapes of the C Auger peaks in graphite and carbide strongly differ, thus the two phases can easily be separated. This way we could distinguish five constituents in the samples, those of C, Ni, Si, Ga and Ni 3 C. We used the usual relative sensitivity factor method for calculation of the concentrations [START_REF] Childs | Handbook of Auger Electron Spectroscopy[END_REF]. To determine the relative sensitivity factors for our spectrometer and measurement conditions for the species above we followed various ways. There were always regions in the samples, which contained only pure C, Ni and Si. Using the Auger intensities measured in these regions we could determine directly the relative sensitivity factors of these elements. In the cases of Ni 3 C and Ga, there were no pure regions, thus other solutions were applied. There were regions in the samples where only Ni 3 C and Ni were found. In these regions the decrease of the Ni signal (the shape of which is not affected by the chemical interaction) with respect to the pure Ni signal can be attributed to the presence of C; this allows us to determine the amount of carbide present and thus the determination of the relative sensitivity factor of the Ni 3 C as well. The relative sensitivity factor of Ga was determined by fitting the measured numbers of Ga atoms to the implanted one in a sample where we could assume that all Ga projectiles remained in the sample. The same relative sensitivity factors were applied for all measured spectra.

In this system the volume density strongly varies (the volume densities of Ni and graphite are 91 and 113 atoms/nm 3 , respectively). Thus, the AES depth profiles will not be given in the usual concentration vs. depth representation; rather the representation of volume densities (atoms/nm 3 ) vs. depth (nm) will be used. Having this representation, one can integrate the depth profile to result in the total numbers of atoms / unit area, which allows the comparison of the calculated number of atoms with those determined by XTEM.

Calculation of the depth scale

The procedure to be applied here has been published recently [18] and will only be shortly summarized as follows. In AES depth profiling, the number of atoms, i t N ∆ , of element i removed from an area A during time t of ion bombardment is

J t Y X A N T i i t * * * * ∆ = ∆ ,
where X i is the atomic concentration of element i, Y T is the total sputtering yield of sputtered material, and J is the ion current density. It is assumed that the intermixed region is an ideal solution and thus the average atomic density of the mixture is

( ) = i i i a X ρ ρ / / 1
, where i is the atomic density of pure element i. The total sputtering yield of the mixture was assumed to be ( )

= i i i T Y X Y
, where the Y i is the partial sputtering yield of element i. Though it is well known that this equation is generally not valid [START_REF] Sigmund | Sputtering by Particle Bombardment I Top[END_REF], we used it since Transport Range of Ions in Matter (SRIM) [START_REF]SRIM Transport Range of Ions in Matter[END_REF] simulation predicted an error of less then 20 %. The number of atoms in the unaltered layers, assuming nominal densities, is known from the XTEM studies, which provides the thickness of layers. The calculated (from AES depth profile) and measured (by XTEM) number of atoms can be equaled for a given layer providing the ion current density.

Using this value, we have the depth scale. The samples have four layers, which have not been affected by the implantation. We always calculated the number of atoms in the unaltered layer as well. Using the above procedure the deviation of the calculated and measured (by XTEM) numbers of atoms in the four unaltered layers were less than 6%, proving the validity of the method.

Results

XTEM results

XTEM imaging was used to characterize the initial and final (irradiated) structures of the samples. Figure 1 shows the XTEM image of the as received sample B. The image shows that the Ni and C layers are polycrystalline and amorphous, respectively. The interfaces are flat.

The intrinsic rms roughness of the interfaces was estimated, based on the XTEM images, to be less than 1 nm. Similar results were obtained in case of the as received sample A as well.

In figure 2 interface, but it is still less than 1 nm. We have carried out many experiments varying the projectiles (Ga + and Ni + ), their energies and the samples. In the following, as examples, some typical results will be summarized.

AES depth profiling results

DISCUSSION

In the discussion we will address two questions a./ how can a sharp C/Ni 3 C rich layer interface form, b./ how can we describe the growth of the Ni 3 C rich layer. In both cases we will deal with a one dimensional problem. These explanations, however, are not applicable for our case, where C with much higher melting point and cohesive energy than those of Ni, seems to be the more mobile species.

Formation of the sharp interface between the carbon and the

Stepanova developed a phenomenological model to describe the asymmetric mixing considering the kinematical factors as well as the changed sputtering yield due to compound formation [START_REF] Stepanova | [END_REF]22]. However, the asymmetry found in our case is so high, that it is unlikely that one could find parameters using of which the given calculation would describe our case. To understand our results we should first emphasize that AES depth profiling as well as the XTEM show the condition of the material after the final stable or metastable state is reached;

thus, we study the condition after ion mixing and possible relaxation. To explain the asymmetrical material transport found in our case let us first consider the measured penetration of the projectiles.

The distributions of the implanted Ga measured on sample B using various fluences at 20 keV irradiation energy are summarized in figure 9. The actually measured Ga densities are multiplied to help the easier comparison of the curves (see later). The shape of the Ga distribution is characterized by a maximum at the interface. In the figure, we also show (by a line) the shape of the Ga distribution according to SRIM simulation (applying the same bombardment conditions on the same sample structure). The simulated Ga penetration also

shows the maximum at the interface. The maximum based on the simulation can be explained by the lower projected range of Ga + in the Ni layer (7.2 nm) than that in the C layer (17.7nm).

Because of the difference of the projected ranges, Ga accumulation at the interface is expected.

The measurement of the implanted Ni distribution is not so straightforward since we cannot distinguish between the implanted Ni and those originally in the Ni layer. Thus, we are restricted to the measurement of the Ni distribution in the C layer. The atomic numbers of Ni [START_REF] Desimoni | [END_REF] and Ga (31) are only slightly different, thus Ni penetration and the mixing it causes should be similar to that of Ga. This is really the case if we consider the ballistic model. In figure 7, we show the Ni distribution (by a line) according to the SRIM simulation. The shape is similar to that of Ga distribution shown in figure 9. The maximum at the interface is explained similarly as before; the projected ranges of Ni + /20 keV are 18.8 and 7.7, in C and Ni, respectively. Thus, we expect that the distribution of the implanted Ni would be similar to that of Ga. However, according to figure 7 Similar explanation was applied for the IBM of Ag/Fe system. In that system the high energy ion irradiation resulted in much less interface broadening than that expected based on the ballistic model; to explain this result, it was supposed that the ballistic mixing is followed by a demixing process [START_REF] Crespo-Sosa | [END_REF]25]. density, increases with increasing fluence. This process is similar to any heterogeneous compound formation on a surface e.g. silicon oxidation [26] or when melting or solidification occurs. In all these cases, the diffusion driven material transport results in a moving boundary [START_REF] Gupta | The classical Stefan problem Noth-Holland series in Applied mathematics and Mechanics[END_REF]. The corresponding differential equation (Stefan problem [START_REF] Gupta | The classical Stefan problem Noth-Holland series in Applied mathematics and Mechanics[END_REF]) is Fick's second equation

The growth of Ni 3 C layer

(in one dimension): 2 2 x c D = t c ∂ ∂ ∂ ∂ ( 1
)
where the boundary conditions at the two interfaces are c is the concentration of the reactant, s(t) is the coordinate of the moving boundary H and k are parameters describing the input and reaction of the reactant at the interfaces (the initial

) ( 0 1 t s = x = x dt ds k = x c ) H(c = x c ∂ ∂ - ∂ ∂
surface is fixed at x=0). According to the equations there is no compound formation in the reacted zone (0<x<s(t)), which is saturated, thus the current of reactant entering at x=0, is completely consumed at x=s(t). The solution of this equation provides s(t), that is, the growth of the reacted zone as a function of time. In steady state condition (∂ 2 c/∂ x 2 =0) the solution of equation (1) gives s(t)~t and s(t)~ t dependencies for low and high t, respectively. In materials science, the most famous application of this equation is most likely the Deal-Grove equation describing the oxidation of silicon [26].

Desimoni et al. [START_REF] Desimoni | [END_REF] applied the same equation for the description of compound formation during ion bombardment. In their approach, adsorption was replaced by the bombardment induced current. They considered asymmetric and symmetric mixing as well. On the other hand, they did not consider the actual concentration distribution, rather only the amount of reacted material has been calculated.

Considering our present experiment, for the first glance, this equation seems to work properly.

A closer look reveals some problems, however. Equation (1) assumes strong compound formation (phase change) at the s(t) interface, and no compound formation in (0<x<s(t)) region, which is saturated. However, in our case a./ the Ni 3 C is a metastable compound with positive heat of formation, b./ in the 0<x<s(t)region there is Ni + Ni 3 C, that is, the reacted region is not saturated by Ni 3 C, c./ the compound-(undamaged) matrix interface is not sharp.

In the following, we show that despite of these, seemingly contradictory facts, our system meets with the requirements above and the layer growth is similar to any heterogeneous reaction.

Formation of Ni 3 C.

The heat of formation of Ni 3 C is in the range of +0.02-+0.05 eV determined by experiment [29], and theory [START_REF] Miedema | Cohesion in Metals: Transition Metal Alloys[END_REF]; thus it is a metastable compound which cannot be formed in equilibrium conditions. Ion beam mixing is far from being an equilibrium process and thus it can overcome the problem above. Really, IBM is frequently applied to produce various nonequilibrium phases like amorphous alloys with negative or positive heat of mixing [1][2][3][START_REF]Beam-Solid Interactions; Physical Phenomena[END_REF][START_REF] Milosavljevic | [END_REF][6][7][8][9]. To explain these phenomena it is frequently argued that the fast quenching is responsible for the process. There are, however, explanations that are more refined. It was observed that multilayer systems of immiscible elements could be converted by IBM to amorphous nonequilibrium alloys [6,7]. In these cases the total energy of the system decreases despite of the positive heat of formation of the alloy, since the mixing removes the high energy interfaces.

We can also find similar explanation for our case. We assume that the metastable Ni Thus, the energy balance is negative, since a high energy vacancy is annihilated in the reaction, and thus the formation of the Ni 3 C compound is energetically favorable.

Layer growth

According to eq. 2 b, compound formation is only possible in regions where vacancies are available. The Ga + and Ni + projectiles evidently produce vacancies. The amount of vacancies is calculated by the ballistic models. The actual vacancy concentration is somewhat uncertain, however, since the ballistic models do not consider the relaxation processes. Experimentally we cannot measure the vacancy distribution. The projectile distribution can be measured, however, and it is accepted that the vacancy distribution correlates with it; thus based on the projectile distribution we can estimate the vacancy distribution. The Ga distributions were measured, and are shown in figure 9. The measured densities are multiplied by factors given in the figure to make the comparison easier.

The depth distributions of Ga produced by various fluences in the range of 100 -400 Ga + /nm 2 (a factor of 4 change in fluence) are rather similar. In case of fluence 600 Ga + /nm 2 a new feature also appears, while at low fluence (50 Ga + /nm 2 ) the signal is too much noisy for comparison.

The similarities of the Ga distributions can be explained by considering that at these bombardment conditions the carbon sputtering yield is low, consequently the thickness of the carbon layer only slightly decreases, which decrease is balanced by the presence of the implanted Ga atoms. Therefore, the penetration of the Ga does not change strongly as a function of fluence, thus the shape of the Ga distribution also does not change as a function of fluence if we are in the given fluence range. That is, the only change with increasing fluence is the density increase of the Ga, maintaining the same shape. At higher fluence ranges, where the removal of the carbon cannot be neglected, the Ga distribution begins to change.

From the point of view of vacancy production we suppose that a threshold concentration of vacancies necessary for compound formation. Since the vacancy concentration correlates with the Ga concentration, the Ga concentration should exceed also a threshold concentration for the compound production. It also follows that the thickness of the compound should increase with increasing fluence, since the necessary threshold density of Ga moves deeper with increasing fluence; consequently there should be a correlation between the Ga penetration and carbide formation. This is really the case. In figure 10 Though the two projectiles are rather similar, SRIM calculation predicts more vacancies in the Ni for Ni + irradiation than that for Ga + irradiation. This also supports the assumption that the number of created vacancies determines the compound formation.

We can also account for the problem of saturation. Obviously if there are no more vacancies available the reaction stops. Thus, the saturation in the 0<x<s(t)region is connected to the lack of vacancies and not to the availability of Ni atoms; the coexistence of the Ni 3 C and Ni does not mean that the region is not saturated.

The fluence dependence of the Ni 3 C can also be explained. Figure 8 shows the total amount of The usual solutions of the Stephan problem predict sharp transitions on both interfaces. These can be accepted if we deal with melting or compound production of high energy of formation.

Even in these cases the discontinuous concentration distribution poses mathematical problems. It has been shown, however, that supposing a homogeneous reaction, a broad interface on the growing side [32] is possible, as in our case.

Summarizing; the observed compound formation was found to be a vacancy assisted solid state reaction, and the kinetics of the growth of the compound layer is described by equation ( 1). This process is expected to occur for other material pairs as well, where the above conditions apply. E.g. we have observed that the irradiation of Si/Cr layer pair by 30 keV Ga + ions also resulted in a broad intermixed layer with sharp interface [33].

In this experiment applying broad or/and overlapping ion beams, we have produced an imbedded layer containing a non-equilibrium compound, which allowed an easy structure and chemical characterization. Obviously, using the same arrangement but applying finely focused ion beam and irradiating not overlapping regions, one can produce nano particles of desired size (limited by the beam size). The lateral arrangement of the imbedded nano particles can follow any tailored pattern because of the versatile control of the FIB.

Conclusions

The interface of C/Ni bilayer was bombarded using Ga

  the XTEM image of sample B irradiated by 200 Ni + (20 keV)/nm 2 is shown. Only the top C and following Ni layers were affected by the irradiation. In the first C layer a broad distribution of heavier atoms (darker contrast would be identified as Ni, see later) appeared. In the first Ni layer we could identify a new layer with lower contrast. The detailed chemical analysis of the sample was carried out by applying AES depth profiling. The rms roughness of the C / new layer interface seemed to be somewhat higher than that of the intrinsic C/Ni

Figure 3 (

 3 Figure 3 (a) shows the whole as measured AES depth profile recorded on sample B irradiated by 200 Ni + (20 keV)/nm 2 . Carbon being in carbide and graphite are shown separately; Ni 3 C

Figure 3 b

 3 Figure 3 b shows, in good agreement with the XTEM image, that the Ni irradiation caused the

First in figure 5

 5 (a) and (b) we show the Ni 3 C density distributions for sample B applying Ni + and Ga + irradiations, respectively. The applied fluences (numbers of ions/nm 2 ) / ion energies (keV) are shown in the figures. The depth scales are shifted to have the C/Ni 3 C transitions at the depth where the original (non-irradiated ) C/Ni interface is. Figures 6 (a) and (b) show the Ni (all Ni being in carbide and metallic form) and C density distributions, for various fluences applying Ga + irradiations. In figure 7, we show the Ni (all Ni in carbide and metallic form) density distributions after Ni + irradiation with various fluences; the fluences (units number of Ni + /nm 2 ) applied are given in the figure. Last, in figure 8, we show the total amount of Ni 3 C /nm 2 formed applying Ni + /20keV, Ga + /20 keV and Ga + /30 keV irradiations of sample B as a function of the square root of the applied fluence. Similar results were obtained if sample A was irradiated.

  Ni 3 C rich layer.

Figure 6 a

 6 Figure 6 a and b shows that there is a sharp interface between the C and Ni (containing Ni 3 C).

  Ballistic model based simulations generally predict that a sharp interface (two pure materials in contact) broadens due to the ion bombardment resulting in an error function concentration distribution. The parameter of the error function increases with fluence. This is exactly the picture suggested by the SRIM[START_REF]SRIM Transport Range of Ions in Matter[END_REF] simulation (based on ballistic model) if we apply it for our case. It predicts that C enters to the Ni layer and Ni appears in the C layer as well. The broadening of the interface is more or less symmetrical. The SRIM simulation cannot account for the accumulation of the defects during the ion bombardment. The more refined TRIDYN (TRIDYN_FZR, FZR-317[START_REF] Tridyn_Fzr | Forschungzentrum Rossendorf[END_REF]) simulation (also based on ballistic model) gives a more detailed description. According to the TRIDYN simulation, Ni penetrates to the carbon layer and its penetration distance and concentration increase with fluence. Thus, the results of both simulations disagree with the experimental finding, which shows a strongly asymmetric material transport.

  the Ni distribution is strongly different from that of Ga. It shows, at any fluence, a maximum within the C layer. Leaving the maximum toward the interface the Ni concentration strongly decreases. Close to the C/Ni 3 C rich layer interface the Ni concentration is low; at lower Ni + fluences it even vanishes. Such behavior is completely absent in the Ga distribution. This type of distribution can neither be explained by the primary penetration of the Ni, nor with ion mixing, which supposes the transport of Ni to the C layer. We can explain this distribution, however, assuming that after the primary processes (primary ion penetration and ion mixing), during the relaxation process, desegregation occurs; all Ni transported into C either as primary ions or mixed atoms in some vicinity of the C/Ni interface, return either to the Ni layer or the Ni enriched region of the C layer. Exactly the same desegregation process should occur in the case of Ga irradiation. Based on the SRIM simulation it is clear that the IBM drives the Ni atoms to the C layer into a 2-3 nm thick layer. If the mixed Ni atoms remain there then if the next cascade hits that region (containing mixed Ni) it drives the Ni atoms even deeper. The TRIDYN simulation predicts such a process, which is resulted in Ni penetration into the C layer. The missing Ni penetration (figures 2,4,6) can only be explained by the Ni desegregation, which should occur before the next cascade.

Figures 5 a

 5 Figures 5 a and b show that the density of the Ni 3 C does not depend on the fluence (except the lowest fluences) applying either Ni + or Ga + projectiles. The (constant) density of Ni 3 C

  3 C forms according to the following chain of reactions: Ni = Ni(int) + Ni(vac) (due to Ga + or/and Ni + bombardment) 2.a 3Ni + Ni(vac) + C (int) = Ni 3 C 2.b adding the two equation gives 4Ni + C = Ni 3 C + Ni(int) 2.c where int and vac stand for interstitial and vacancy, respectively. The experimentally reported formation energy of the self vacancy of Ni is in the range of 1.3-1.85 eV [31]. The formation energy of Ni 3 C is +0.02-0.05 eV [29,30]. We have not found data for the energy of the C interstitial in Ni, but it should be positive. The energy of Ni(int) is supplied by the ion beam.

  we show the Ga and Ni 3 C density distributions as a function of depth for fluences ranging from 50 Ga + /nm 2 to 600 Ga + /nm 2 . To show the relative position of the Ga and Ni 3 C distributions, the distributions are shifted to have the same depth for the Ni 3 C/pure Ni transitions. It is clear that relative position of the two distributions is unchanged while the fluence increases by more than a factor of 10. This correlation supports the assumption that compound formation occurs if the vacancy concentration exceeds a threshold value. The same vacancy-mediated Ni 3 C formation process is responsible for the projectile dependent Ni 3 C density; figures 5 a and b show that applying the same bombarding energy Ni irradiation results in higher concentration of the Ni 3 C than that produced by Ga irradiation.

Ni 3 C

 3 figure 8 might be also explained by equation (1).

  + and Ni + projectiles of energies in the range of 20-30 keV. The ion bombardment resulted in a asymmetric material transport, C had been transported to Ni, but not vice versa. Because of the asymmetric material transport, the carbon/Ni 3 C rich layer interface was found to be sharp. The carbon transported to the Ni layer formed metastable Ni 3 C compound. The concentration of the Ni 3 C compound depended on the type of projectile. With increasing fluence only the thickness of the Ni 3 C rich layer increased but not the concentration of the Ni 3 C. The total amount of Ni 3 C produced by the irradiation was found to be proportional to the square root of the fluence at lager fluence values. The sharp carbon/Ni 3 C rich layer interface could be explained by relaxation processes; the Ni transported into the C layer by the cascade mixing, returns to the Ni layer by means of desegregation. Metastable Ni 3 C formation was explained by a vacancy mediated process, while the fluence dependence of the layer growth can be described by diffusion equation considering a moving boundary.

Figure captions 1 . 4 .FiguresFigure 1 .

 141 Figure captions

Figure 2 .

 2 Figure 2. XTEM image of irradiated sample B. Irradiation: Ni + , 20 keV, 200 ions/nm 2 .

Figure 3 .

 3 Figure 3. a The whole as measured AES depth profile (Auger signal vs. sputtering time) of sample B irradiated by Ni + (20 keV, 200 ions/nm 2 ).

Figure 3b .

 3b Figure 3b. The beginning part of the profile shown in figure 3 a.

Figure 4 .

 4 Figure 4. The as measured AES depth profile of sample B irradiated by Ga + (20 keV, 200 ions/nm 2 ).

Figure 5 .

 5 Figure 5. a Ni 3 C density profile measured on sample B; projectile Ni + (20 keV). The parameters in the figure are the fluence (unit, numbers of ions/nm 2 )/ ion energy (keV).

Figure 5 .

 5 Figure 5. b. Ni 3 C density profile measured on sample B; projectile Ga + . The parameters in the figure are the fluence (unit, numbers of ions/nm 2 )/ ion energy (keV).

Figure 6 .

 6 Figure 6. The Ni (a) and C (b) density profiles in case of Ga + (20 keV) irradiation. The parameters in the figure are the fluences (unit, numbers of Ga + ions/nm 2 ) applied.

Figure 7 .

 7 Figure 7. The Ni density profiles for Ni + (20 keV) irradiation. The parameters in the figure are the fluences (unit, numbers of Ni + ions/nm 2 ) applied. The line shows the shape of the Ni distribution according to SRIM simulation.

Figure 8 .

 8 Figure 8. The total amount of Ni 3 C / nm 2 produced by various projectiles and energies as a function of the square root of fluence.

Figure 9 .

 9 Figure 9. The Ga density profile in sample B at 20 keV Ga + bombardment. The measured density values are multiplied by the factors shown in the figures to make the comparison easier. The parameters in the figure are the fluences (unit, numbers of Ga + ions/nm 2 ) applied. The line shows the shape of the Ga distribution given by SRIM simulation.

Figure 10 .

 10 Figure 10. The Ga (line) and carbide (symbol) density profiles. Parameters in the figure are the fluences (unit, numbers of Ga + ions/nm 2 ) applied. Sample B, Ga + 20 keV. The distributions are shifted to have the same depth for the Ni 3 C/pure Ni transitions.

Acknowledgements

Ni irradiation was made possible by the AIM -Center for Application of Ion Beams to Materials Research, RITA Contract No. 025646.


Anton Krynkin
email: a.krynkin@salford.ac.uk

Olga Umnova

Alvin Y B Chong

Shahram Taherzadeh

Keith Attenborough

Scattering by coupled resonating elements in air

Keywords: Sonic crystals, Scattering, Helmholtz resonance, Elastic shell, Band gaps PACS numbers: 43.20.+g, 43.50.+y, 46.70.De, 43.20.Ks

Scattering by (a) a single composite scatterer consisting of a concentric arrangement of an outer N-slit rigid cylinder and an inner cylinder which is either rigid or in the form of a thin elastic shell and (b) by a finite periodic array of these scatterers in air has been investigated analytically and through laboratory experiments. The composite scatterer forms a system of coupled resonators and gives rise to multiple low frequency resonances. The corresponding analytical model employs polar angle dependent boundary conditions on the surface of the N-slit cylinder. The solution inside the slits assumes plane waves. It is shown also that in the low-frequency range the N-slit rigid cylinder can be replaced by an equivalent fluid layer. Further approximations suggest a simple square root dependence of the resonant frequencies on the number of slits and this is confirmed by data. The observed resonant phenomena are associated with Helmholtz-like behaviour of the resonator for which the radius and width of the openings are much smaller than the wavelength. The problem of scattering by a finite periodic array of such coupled resonators in air is solved using multiple scattering techniques. The resulting model predicts band-gap effects resulting from the resonances of the individual composite scatterers below the first Bragg frequency . Predictions and data confirm that use of coupled resonators results in substantial insertion loss peaks related to the resonances within the concentric configuration. In addition, for both scattering problems experimental data, predictions of the analytical approach and predictions of the equivalent fluid layer approximations are compared in the low-frequency interval.

Introduction

Infinite periodic structures do not support wave propagation in certain frequency intervals known as stop bands (band gaps) [1,2]. For finite arrays of periodically arranged scatterers, referred to as Sonic Crystals, the stop bands correspond to the frequency intervals of very low transmission [3,4,5]. The position of the stop bands can be tuned to selected frequency intervals by changing the spacing between the scatterers. This makes Sonic Crystals attractive for applications as noise barriers particularly for narrow band sources.

The performance of the Sonic Crystals can be improved by increasing the filling fraction [2,6], varying the arrangement of the scatterers [7] and also by replacing the scatterers with the resonant elements [8,9,10,11,12,13,14,15,16]. The latter allows increasing insertion loss in the lowfrequency range if the resonant frequencies of the scatterers lie below the first stop band associated with the array periodicity. The most obvious choice for such resonant elements are split ring resonators [9,13,15] which essentially are 2D Helmholtz resonators. Alternatively an array of thin elastic shells in air has been shown to possess strong multiple resonances in the low-frequency range [11,16].

Arrays of composite scatterers such as concentric split ring resonators [9,17] are shown to support band gaps generated by the multiple resonances of each scatterer in the array. In this paper we investigate the concentric arrangement of two types of the resonant elements i.e. thin elastic shells and split ring resonators with multiple slits, referred to as N-slit cylinders. These composite scatterers, referred to as composites, serve a double purpose. First, they support multiple resonances including Helmholtz and annular cavity resonances and also axisymmetric resonance of the elastic shells. The latter is defined by the geometry of the shell and not by the area of the cavity that therefore gives advantage in tuning of the composite multiple resonant behaviour compared to that for the concentric split ring resonators. It is also shown that the coupling between two resonators leads to the shift of the axisymmetric shell resonance to the lower frequency range. The analogous effect for coupled mechanical oscillator (multiple degree of freedom system) is well known [START_REF] Morse | Theoretical Acoustics[END_REF]. Second, for the purpose of creating a noise barrier an array of elastic shells exposed to the environment would not be practical. A simple solution is to protect the shells by using outer cylinders containing slits surrounding inner elastic shells arranged concentrically.

The scattering problem for the hollow cylinder with non-uniform boundaries has been solved semi-analytically in [START_REF] Tsalamengas | [END_REF]20] and references therein. This involves solving the linear algebraic system of equations derived from the analytical evaluation of the integral equations. An alternative approach has been developed in grating theory [21,22] where solution of the studied partial differential equations is subject to the boundary conditions dependent on the coordinates.

Three configurations involving the N-slit cylinder are studied in this paper (a) alone, (b) concentric with an inner rigid cylinder and (c) concentric with an inner elastic shell. To solve the scattering problem for N-slit cylinders we use boundary conditions dependent on polar angle [21]. In section 2 this approach is employed to solve a single composite scatterer problem that is similar to that solved for an electromagnetic case [22]. However in the present model the finite thickness of the wall of the N-slit cylinder is taken into account. This avoids the use of adjustable parameters introduced for numerical stability [21,22] and provides a more accurate description of the real structure used in the experiments reported later. The solution inside slits is replaced by jump conditions [23] that describe the slit interface as a moving piston. It should be noted that use of jump conditions to replace the slits makes it possible to accurately predict the total wave field in the low-and mid-frequency range. It is also shown that in the low-frequency range the slit cylinder can be replaced by an equivalent fluid layer. The approximations are based on the results derived for the perforated plate [24,25,[START_REF] Shenderov | Wave Problems of Underwater Acoustics[END_REF]. The developed approximations are of particular interest for studying the low-frequency resonant behaviour of the single composite arrangement and a finite periodic array of them. This behaviour is attributed to the Helmholtz-like resonances with wavelength much bigger than the characteristic sizes of the cylinder (i.e. radius and width of the slits).

Experimental validation of the analytical and numerical predictions has been carried out. In particular the results obtained with the grating theory are compared with the low-frequency approximations for the cylinders with different number of slits and various inner structure (N-slit cylinder alone, rigid core and elastic shell core). However the analytical models could be used also to study the acoustical properties of arrays incorporating resonators of various thickness and with different slit width.

In section 3 the model is generalised for finite arrays of N-slit cylinders and composite elements. In section 4 the experimental setup is described and the model predictions are compared with the data.

Single scatterers

Formulation

Consider the two-dimensional problem of acoustic wave scattering by a single N-slit rigid cylindrical shell of thickness h and external radius r o . The sound is generated by the cylindrical point source placed at the origin of the coordinate system which is defined by either the Cartesian (x, y) or polar (r, θ) coordinates. Throughout the paper the time-harmonic dependence is taken as exp(-iωt). The widths of consecutive slits in Oxy plane are denoted by d n , n = 1..N and they are infinitely long in the direction of the cylinder main axis Oz, see Figure 1(a). The external length of the rigid arc is given by D n . It is also assumed that the thickness of the rigid shell is much smaller than its radius so that the following geometric simplifications can be applied: Without loss of generality it is assumed that radius of the vector passing through the middle of the first slit makes either zero or π angle with the Ox axis.

For simplicity, the acoustic environments outside, inside of the N-slit cylinder and in the slits are assumed to be identical and are described by density ρ and sound speed c. We consider three configurations: (a) an N-slit cylinder alone, (b) with a rigid cylindrical core inside and (c) with an elastic shell inside.

The displacement potential p(r) in the acoustic medium satisfies the Helmholtz equation

∆p α + k 2 p α = 0, (1)
where k = ω/c is the wavenumber, index α relates solution p to one of the regions (i.e. 'o' outside the configuration, 'i' inside it and 's' within the slits) and Laplacian ∆ is given by either 1 r

∂ ∂r r ∂ ∂r + 1 r 2 ∂ 2 ∂θ 2 or ∂ 2 ∂x 2 +
∂ 2 ∂y 2 . The outer solution p o (r) also satisfies Somerfeld's radiation conditions

∂p o ∂r -ikp o = o r -1/2 , as r → ∞, (2)
where r = x 2 + y 2 .

To proceed with the boundary conditions imposed on the surface of the N-slit cylinder we first solve an auxiliary problem for the slits. The geometry of this problem is illustrated in Figure 1(b). Waves propagating inside the slit are described by plane wave solution as

p s = C 1 e iky + C 2 e -iky , (3)
where C 1 and C 2 are unknown coefficients. It is assumed that contribution of the x-dependent components to p s is negligible when slits are small compared to the wavelength (that is kd n 1 and kh 1). Solution (3) is subject to the continuity conditions imposed on the slit faces y = ±h/2 as

p α = p s , q α = ∂p s ∂y , (4)
where q α = ∂p α /∂n is the normal derivative on the slit faces and index α = o, i. Using solution (3) and its unknown coefficients C 1 and C 2 , equations (4) can be transformed to the jump conditions [23]

p i = p o -hq o , q i = k 2 hp o + q o , (5)
that relate the wave field inside the N-slit cylinder to that of the outer region at the slit interface. Note that due to assumption kh 1 the trigonometric functions in ( 5) are replaced by their leading orders.

The solution of equation ( 1) is subject to the jump conditions (5) as well as to the Neumann conditions imposed on the rigid surface of the N-slit cylinder. The former and the latter can be combined into the set of two boundary conditions [21] that are

∂p i ∂r = ∂p o ∂r + f (θ)k 2 hp o , ∂p o ∂r - f (θ) h (p o -p i ) = 0, (6)
where stepwise function f (θ) of the polar angle θ introduces the distribution of slits along the N-slit cylinder surface over the interval θ ∈ [0, 2π] as

f (θ) = 1 if θ ∈ [-φ 1 /2, φ 1 /2] ∪ .. ∪ [2π -φ 1 /2 -φ N -ψ N , 2π -φ 1 /2 -ψ N ], 0 otherwise. ( 7
)
For the non-zero value of this function equations ( 6) are transformed to the jump conditions ( 5) whereas the Neumann conditions can be obtained from ( 6) by setting f (θ) to zero.

The periodic angle distribution of N slits can be introduced through the alternative form of f (θ) that is

f (θ) = N n=2 H θ - 2π(n -1) N + φ n /2 -H θ - 2π(n -1) N -φ n /2 + 1 j=0 {H (θ -2πj + φ 1 /2) -H (θ -2πj -φ 1 /2)} , (8)
where H(θ) is the Heaviside function.

The total wave field outside the slit cylinder is represented by p o which can be found as

p o = H (1) 0 (kr) + +∞ n=-∞ A n H (1) n (kr)e in θ , (9)
where r(r, θ) and θ(r, θ) are the polar coordinates of receiver with origin placed at the centre of the scatterer and A n are unknown coefficients. The solution inside the slit cylinder is given by

p i = +∞ n=-∞ [B n J n (kr) + C n Y n (kr)] e in θ, (10)
where coefficients B n and C n have to be found. The coefficient C n can be derived in terms of B n or set to zero according to the type of the core that results in the following

C n = B n C n , n ∈ Z, (11)
where

C n = 0, N-slit cylinder alone, ( 12
)
C n = - J n (ka 1 ) Y n (ka 1 )
, rigid core, ( 13)

C n = - [J n (kR)] 2 [1 -(k s R) 2 + n 2 ] J n (kR)Y n (kR)[1 -(k s R) 2 + n 2 ] + [n 2 -(k s R) 2 ]ρ(ρ s πRh) -1 , elastic shell core, ( 14
)
where primes (') denote the derivative with respect to polar coordinate r. In equation ( 14) the elastic shell is described by frequency parameter k s = ω/c s , density ρ s , dilatational wave speed c s , radius of the elastic shell mid-surface R = (a 1 + a 2 )/2 and its half thickness h s = (a 1 -a 2 )/2. Here the dilatational wave speed c s for a thin elastic plate is given by

c s = E ρ s (1 -ν 2 ) , (15)
with Young's modulus E and the Poison's ratio ν. In this paper we consider thin viscoelastic shells made of latex [16]. For the numerical predictions the size of each elastic shell is always defined by the outer radius a 1 = 0.02 m and thickness 2h s = 0.00025 m.

In order to find the unknown coefficients A n we apply Graf's addition theorem [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] to the outer solution (9). This enables us to express (r, θ) in terms of (r, θ) [16]. Then, substituting ( 9) and (10) into the boundary conditions (6) and taking the inner product A n δ m,n 2πhH (1) n (kr o )

-F n-m H (1) n (kr o ) -H (1) n (kr o )I n + k 2 h 2π H (1) n (kr o ) ∞ j=-∞ F j-m F n-j I j = (16) ∞ n=-∞ H (1) n (kQ)e -in(π+α) -δ m,n 2πhJ n (kr o ) + F n-m [J n (kr o ) -J n (kr o )I n ] - k 2 h 2π J n (kr o ) ∞ j=-∞ F j-m F n-j I j , m ∈ Z,
where δ m,n is Kronecker delta, vector Q = Q(cos α, sin α) is the radius vector to the centre of the slit cylinder,

I n = J n (kr i ) + C n Y n (kr i ) J n (kr i ) + C n Y n (kr i ) , (17)
and F n is the Fourier transform of function f (θ) given by

F n =            N l=1 φ l for n = 0, 2 n N -1 l=0 sin nφ l 2 e -2inπl/N for n = 0. (18)
Note that factor C n and the geometrical parameters of the slits only appear in I n and F n , respectively. This makes form of the system (16) invariant with respect to scatterer core type and the arrangement and size of the slits.

Taking φ l = 0, l ∈ Z, the solutions of system ( 16) is reduced to the case of scattering by rigid cylinder that is

A m = - J m (kr o ) H (1)
m (kr o ) To find the numerical solution the infinite system (16) has to be truncated at some number m, n, j = -M..M . In general the convergence of the numerical solution is dependent on the radius r o of the N-slit cylinder, number of slits and their angles φ l (i.e. their length) as well as on the frequency. Reducing the angle φ n of the slit results in faster convergence so that numerical solution approaches the value defined in [START_REF] Tsalamengas | [END_REF]. It is found that for kr o < 5 and φ l < 0.1, l ∈ Z, truncation number M between 30 and 40 gives accurate results to less than three significant figures. Throughout this paper we use M = 35.

H (1) m (kQ)e -im(π+α) , m ∈ Z. ( 19
)

Results

Knowing the coefficients A n , n ∈ Z, enables us to compute scatterer insertion loss given by IL = 20 log 10 H

(1)

0 (k o r) p o , (20)
with acoustic potential p o found by formula (9). The insertion loss computed for a single N-slit rigid scatterer is compared with data in Figure 2 (the experimental setup is described in section 4). The results illustrated in Figures 2(a),(b) and(c) are obtained for N-slit cylinders alone and with a rigid core. For the given scatterer the prediction exhibits Helmholtz resonance in the low-frequency range (i.e. kr o < 1). In both semi-analytical results and data, this resonant behavior is observed around 600 Hz for a single slit rigid cylinder, see Figure 2(a), and around 1500 Hz for a 4-slit rigid cylinder, see Figure 2(b). It is observed that increasing the number of slits results in significant shift of this Helmholtz-type resonance toward higher frequencies. The resonance can also be shifted by reducing the cavity area. This effect is seen in Figure 2(c) where the presence of a concentric internal rigid cylinder results in a shift of the resonance to a higher frequency f ≈ 1700 Hz.

The second type of the resonant behaviour observed in Figure 2 is related to the cut off frequency of the acoustic mode in the circular/annular cavity. For the single slit rigid cylinder this resonance is found around 4200 Hz (Figure 2(a)). The 4-slit rigid cylinder (Figure 2(b)) exhibits this type of resonance at a slightly higher frequency f ≈ 4400 Hz. For the annular cavity formed by the concentric 4-slit rigid cylinder and rigid internal cylinder one of the cut off frequencies can be found around 3500 Hz (Figure 2(c)).

Figure 2(d) illustrates predicted and measured insertion losses for a 4-slit rigid cylinder with an inner elastic shell. The presence of the inner elastic shell with radius a 0 = 0.02 m causes resonant behavior to appear around 1000 Hz, 2700 Hz and 3400 Hz. The first two resonances can be related to the interaction of the axisymmetric resonance of the elastic shell and the Helmholtz resonance. The frequency of the axisymmetric resonance of the elastic shell, which is originally at about 1200 Hz, is decreased to about 1000 Hz by the concentric arrangement. The second resonance at f ≈ 2700 Hz is related to the Helmholtz resonance. This frequency is higher than that for the empty slit cylinder due to the effective reduction in the size of the cavity.

Low-frequency approximations

In the low-frequency range a simpler approximation can be used to describe scattering by Nslit cylinders. The approximation is based on the results derived by Horace Lamb [24] for the electromagnetic waves transmitted through a metallic grating. Here we only consider a case when N identical slits of width d n = d, n = 1..N, are distributed periodically along the scatterer surface. In this case N-slit cylinder can be replaced by a homogeneous fluid shell. The latter is characterised by two effective parameters that are density ρ l and speed of sound c l .

Similar approximations have been used for a periodic array of circular scatterers [START_REF] Hu | [END_REF]. The fluid layer is defined through the thickness h l = hF and density ρ l = ρ/F 2 with filling fraction F = Nd/(2πr o ) satisfying the following conditions

F < 1 4 and λ = c f > F. (21)
The parameters h l and ρ l were obtained by matching transmission coefficients derived for the infinite homogeneous and perforated plates. Having thickness as an effective parameter for the infinite plate gives quite accurate results in the low-frequency regime. However the fact that h l is not equal to the actual thickness of the element leads to unnecessary complications in the model for the cylindrical shell.

In contrast here the thickness of the equivalent fluid layer is made equal to the actual thickness of the N-slit cylinder and, therefore, the inner and outer radii of the corresponding fluid shell are uniquely defined through those of the N-slit cylinder. The obtained approximation gives first resonance position within 5% of that for the N-slit cylinder provided that length of d satisfies conditions (21). 9)-( 16) (---) for single N-slit scatterer compared with its approximations (----) defined by ( 29), ( 30 To find the unknown density ρ l and speed of sound c l we first write down the plane wave transmission coefficient of a fluid layer of thickness h as

T = 4Z -1 e -ikh (1 + Z -1 ) 2 e -ik l h -(1 -Z -1 ) 2 e ik l h , Z = ρ l c l ρc , ( 22
)
where k l is the wave number in the fluid layer and Z is its relative impedance. Note that coefficient T is derived for the plane waves propagating perpendicular to the surface of the layer. In the low frequency regime when kh 1 and k l h 1 the transmission coefficient can be approximated by

T ≈ Z -1 Z -1 -0.5ik l h(1 + Z -2 ) . ( 23
)
For a perforated plate such that kd 1 the transmission coefficient T is described by oscillating piston-type waves propagating in slits. Using existing analytical results [START_REF] Shenderov | Wave Problems of Underwater Acoustics[END_REF] we arrive at

T = 4Fe -ikh (1 + Ẑ-1 ) 2 e -ikh -(1 -Ẑ-1 ) 2 e ikh , ( 24
) with Ẑ = F - ikd F 2 π 3 ∞ n=1 sin(Fπn) n 3 = F -ik∆, ( 25
)
where ∆ is referred to as end correction. The lower limit of the end correction can be estimated from [24,25]

∆ = d π log 1 sin(πF/2) . ( 26
)
By assuming low frequency and small filling fraction (i.e. F 1) the coefficient T is approximated as

T ≈ F F -0.5ikh (1 + F 2 + 2∆/h) . ( 27
)
In order to be able to replace the perforated plate by an effective fluid layer we equate T and T . Comparison of these two coefficients gives relative impedance Z and the relationship between wave numbers k l and k in the following forms

Z = 1 F k l = k h h + 2∆ 1 + F 2 ≈ k h (h + 2∆) . ( 28
)
It is then possible to derive the density and sound speed of the effective fluid layer as

c l = h h + 2∆ c ( 29
)
ρ l = h + 2∆ hF ρ. (30)
To derive the solution of the appropriate scattering problem for the layered cylinder continuity boundary conditions have to be imposed at the faces r = r o , r i of the fluid layer in the following form

p l = - ρ ρ l p α ∂p l ∂r = ∂p α ∂r , (31)
where solutions p α are given by ( 9) and ( 10) whereas solution inside the fluid layer is defined through the new set of constants D n and E n , n ∈ Z, as

p l = +∞ n=-∞ [D n J n (k l r) + E n Y n (k l r)] e in θ . (32)
Note that conditions (31) and solution (32) replace the boundary conditions ( 6) imposed at the faces of the N-slit cylinder.

Using boundary conditions (31) the unknown coefficients A n , D n ,E n and C n can be found. The coefficient A n takes the following form

A n = -Z n H (1) n (k o Q)e -in(π+α) , (33)
where

Z n = ρ l W 1 W 3 J n (kr o ) -ρJ n (kr o ) [W 1 J n (k l r o ) + W 2 Y n (k l r o )] ρ l W 1 W 3 H n (r o ) -ρH n (r o ) [W 1 J n (k l r o ) + W 2 Y n (k l r o )] , (34)
within which

H n (r o ) = H (1) n (kr o ), W 1 = [-ρY n (k l r i )Y n (kr i ) + ρ l Y n (k l r i )Y n (kr i )] C n + ρ l J n (kr i )Y n (k l r i ) -ρJ n (kr i )Y n (k l r i ), W 2 = [ρJ n (k l r i )Y n (kr i ) -ρ l J n (k l r i )Y n (kr i )] C n -ρ l J n (kr i )J n (k l r i ) + ρJ n (kr i )J n (k l r i ), W 3 = W 1 J n (k l r o ) + W 2 Y n (k l r o ) W 1 , ( 35
)
with C n defined in (12).

To validate the low-frequency approximation its predictions are compared with those of the model presented in section 2.2. In Figure 3 the results are shown for the same types of scatterers. The accuracy of the approximation deteriorates with increase in frequency. However in the frequency range that is inside of the limit ( 21) the approximation reproduces the shift of the resonances observed for the concentric arrangement in Figures 3(c) and(d).

The low-frequency resonance displayed in Figure 3 below 2000 Hz can be found by constructing approximation of the appropriate eigenvalue solution. In the eigenvalue problem the outer acoustic potential ( 9) is defined through Bessel functions of second kind only. The solution of this problem corresponds to zeros of the determinant of coefficient A n in (33) with H n (r o ) = Y n (kr o ). To obtain the approximations of the resonant frequencies we need to consider two different cases defined by the core type factor C n in ( 12)-( 14). For N-slit cylinders that are either empty or concentric with an inner rigid cylinder the leading order approximations can be derived by assuming that

kr o = O(β), k l r o = O(β), ρ ρ l = O(β 2 ), (36)
where β 1 is a small parameter. These four relations define the frequency regime and contrast between fluid layer and outer/inner environment. Note that according to equation (30) and the assumption that F 1, density ρ l is always much greater than ρ. Expanding Bessel functions and collecting similar orders of smallness the first resonant frequency i.e. n = 0 can be approximated by

f = c 2π 2ρ ρ l (r 2 i -a 2 1 ) log(r o /r i ) . (37)
The resulting values of resonant frequency can be further approximated by using the fact that ratios h/r i and h/r o are small. This follows from the assumption made in the beginning of section 2.1. The expansion log(r o /r i ) = h/r i + O(h 2 /r 2 i ) and substitution of the equation (30) give

f ≈ c 2π Nd π(r 2 i -a 2 1 )(h + 2∆) . ( 38
)
For a single slit N = 1 equation (38) reduces to the well-known value of Helmholtz resonance. From equation (38) it is also seen that increasing the number of slits leads to the shift of the Helmholtz resonance toward higher frequencies by approximately √ N . This is in accordance with the results observed in Figures 3 (a Next is considered the case of the concentric N-slit cylinder and elastic shell which is described by factor C n in equation (14). The interaction of the axisymmetric and Helmholtz resonances leads to the considerable shift of the low-frequency resonance as well as the high-frequency resonances. The new position of the low-frequency resonance can be approximated by introducing the following order of smallness of the physical and geometrical parameters

kr o = O(β), k l r o = O(β), ρ ρ s = O(β 3 ), ρ l ρ s = O(β), c s c = O(β), c s c l = O(β), h s R = O(β), (39)
For n = 0 the leading order approximation of the resonant frequency becomes

f = c 2π 1 2h s R(r 2 i -R 2 ) log(r o /r i ) ρ ρ s r 2 i + c 2 s c 2 h s R (r 2 i -R 2 ) log r o r i + 2 ρ ρ l h s R ± ρ ρ s r 2 i + c 2 s c 2 h s R (r 2 i -R 2 ) log r o r i -2 ρ ρ l h s R 2 + 8 ρ ρ l ρ ρ s h s R 3 log r o r i 1/2 . (40)
In this equation the presence of two solutions explains the interaction of the first two low-frequency resonances. Two limiting cases can be constructed which show the degeneration of the two resonances into single one. It can be seen that if radii r i and r o tend to infinity the equation (40) transforms to the resonance of the elastic shell under membrane compression loading [16] i.e.

f = c 2π 1 R c 2 s c 2 + R h s ρ ρ s . (41)
Comparison of the resonant frequency predicted by (41) for an elastic shell made of latex with the smallest value in (40) obtained for the concentric 4S rigid cylinder and latex elastic shell enables us to estimate the shift of the resonance observed in the composite. The numerical value of the axisymmetric resonance (41) is f ≈ 1270 Hz whereas the lowest axisymmetric resonance of the composite is f ≈ 1060 Hz. Thus a reduction of the frequency of the axisymmetric elastic shell resonance by about 200 Hz is predicted. Equation (40) can also be reduced to the form of the resonant frequency (37) by assuming higher contrast in (39) between the fluid and elastic media.

Arrays of scatterers

Formulation

The method developed in the previous section can now be adapted to finite arrays of scatterers. We first formulate the method of multiple scattering that is described by the superposition of scattered wave fields of each element of the array [5] and boundary conditions (6) imposed on the surface of all scatterers in the array. As a result the solution p o of the Helmholtz equation (1) in the acoustic environment outside scatterers takes the following form

p o (r, θ) = H (1) 0 (kr) + N m=1 +∞ n=-∞ A m n H (1) n (kr m ) exp(in θm ), ( 42
)
where N is the number of scatterers in the array, variables rm (r, θ), θm (r, θ) are the polar coordinates with origin in the centre of scatterer of index m and A m n are unknown coefficients. The form of the solution inside of the m-th slit cylinder is the same as in the case of single scatterer problem that is equation (10).

Using Graf's addition theorem we can rewrite solution (42) in terms of polar coordinates (r m , θm ) of m-th scatterer that gives

p o (r m , θm ) = +∞ n=-∞
J n (kr m )H (1) n (kQ m )e -in(π+αm) + A m n H (1) n (kr m )

+ N p=1, p =m +∞ q=-∞ A p q J n (kr m )H (1)
q-n (kQ mp )e i(q-n)(αmp+π) e in θm , m = 1..N

where vector Q m = Q m (cos α m , sin α m ) is the radius vector to the centre of m-th scatterer, vector Q mp = Q mp (cos α mp , sin α mp ) defines the position of p-th scatterer with respect to m-th scatterer and Q m , Q mp > rm that is the requirement of addition theorem. Expansion (43) along with the inner solution (10) are subject to the boundary conditions (6). This results in the infinite algebraic A p n δ m,n 2πhH (1) n (kr o,p )

-F p n-m H (1) n (kr o,p ) -H (1) n (kr o,p )I p n + k 2 h 2π H (1) n (kr o,p ) ∞ j=-∞ F p j-m F p n-j I p j + N s=1, s =p ∞ n=-∞ ∞ v=-∞ A s n H (1)
n-v (kQ ps )e -i(n-v)(π+αps) δ m,v 2πhJ n (kr o,p )

-F p v-m [J n (kr o,p ) -J n (kr o,p )I p n ] + k 2 h 2π J v (kr o,p ) ∞ j=-∞ F p j-m F p v-j I p j = ( 4 4 )
∞ n=-∞

H (1) n (kQ p )e -in(π+αp) -δ m,n 2πhJ n (kr o,p )

+ F n-m [J n (kr o,p ) -J n (kr o,p )I p n ] - k 2 h 2π J n (kr o,p ) ∞ j=-∞ F p j-m F p n-j I p j , m ∈ Z, p = 1..N ,
where factors I p m and F p m are given by equations ( 17) and [START_REF] Morse | Theoretical Acoustics[END_REF], respectively. The superscript p relates these factors to p-th scatterer. To find unknown coefficients we again use truncated system of N (2M + 1) algebraic equations. Note that for the slit angles φ p l = 0, l ∈ Z, p = 1..N , equations (44) are transformed to the well-known algebraic system describing multiple scattering problem for an array of rigid cylinders [5].

Results

Solving the algebraic system of equations (44) enables us to compute the total wave field (42) and insertion loss (20). The latter is analysed in this section. Figure 4(a) shows measured and predicted insertion loss spectra for a 7 × 3 array of 4-slit cylinders. The band gap effect related to the cavity resonance (see Figure 2) and associated with a high insertion loss peak is predicted and observed in the vicinity of f = 1500 Hz. This effect is followed by the insertion loss peak related to the first Bragg band gap which is observed around f ≈ 2000 Hz.

In Figure 4(b) the insertion loss is computed for the array of composite scatterers consisting of concentric 4S rigid cylinders and latex elastic shells. The coupled resonances of this composite described in section 2.2 generate various band gaps that result in high insertion loss peaks. One of these peaks is observed around 1000 Hz and is related to the band gap due to the shifted axisymmetric resonance of latex elastic shell. The insertion loss peak associated with the shifted Helmholtz resonance is observed around 2600 Hz and it follows the first Bragg band gap peak observed around f ≈ 2000 Hz.

Low-frequency approximations

Here the model of equivalent fluid layer is applied to the problem of multiple scattering. Using the low-frequency approximations (29) and (30) together with the boundary conditions (31) one can derive the analogue of the algebraic system of equations (44) that is

A m n + N p=1, p =m ∞ q=-∞ A p q Z p q H (1) q-n (k o Q mp )e i(q-n)(π+αmp) = -H (1) n (k o Q m )e -in(π+αm) , (45)
where Z p q is given by (34) for the scatterer with index p = 1..N [5]. The truncated linear system of N (2M + 1) equations is solved to find the unknown coefficients A m n and, as a result, the insertion loss (20). The truncation number M is taken between 5 and 7. This results in substantial reduction in the computational time compared to that for the full problem (44).

In Figure 5 insertion loss is computed for the same arrays as considered in section 3.2. It is observed that replacement of the N-slit rigid cylinder by the equivalent fluid layer gives results accurate to within 5% in the frequency interval f ∈ (0, 2000) Hz. Therefore, for the given geometry of the array the proposed approximations accurately predict the presence of the band gaps below and including the first Bragg band gap. The band diagrams plotted in Figure 6(a) and (b) are aligned with the insertion loss spectra obtained previously in Figure 5. It is demonstrating that the band gaps coincide with the maxima of insertion loss. The diagrams also show that the first band gaps (at 1500 Hz in Figure 6(a) and at 1000 Hz in Figure 6(b)) are complete. This results in angular independent insertion loss peaks bounded in the frequency interval of the complete band gaps.

Laboratory measurements

Cylinder constructions

2 m long 0.25 mm thick Latex sheets have been formed into cylinders with outer diameter of 43 mm. This has been achieved by overlapping edges by a few mm and gluing them together. Slits have been gouged from the walls of 55 mm outer diameter PVC pipes. To ensure the structural integrity of the pipes, rather than making continuous slits along the complete length of the pipes, the slits were made in sections of approximately 20 cm length separated by about 2 cm (see Figure 7(a)). Concentric arrangements of pairs of latex and 4-slit PVC cylinders were formed as shown in Figure 7(b). To ensure the central location of each 27.5 mm outer diameter latex cylinder within the associated 4-slit PVC pipe, each end of the latex cylinder was secured on to a Perspex cap which fitted inside the PVC pipe. Smaller PVC pipes with inner diameters of 22 mm were used inside the 4-slit pipes for some measurements.

Measurement system and data analysis

The sound source was a Bruel & Kjaer point source loudspeaker controlled by a Maximum-Length Sequence System Analyzer (MLSSA) system enabling determination of impulse responses. Measurements were taken of the insertion loss (IL) spectra for single cylinders and arrays of cylinders in an anechoic chamber. Figures 8(a for the 2 m long cylinders were provided by holed wooden boards at the top and base of each array. The lattice constant (L) for the arrays of 4-slit cylinders and concentric cylinders was 80 mm. To maintain their shape and vertical orientation, the latex cylinders were slightly inflated above atmospheric pressure through a common pipe connecting to a small pump. During the array measurements, the receiver microphone was positioned 50 mm from the nearest face of the array but on the opposite side to the source. The loudspeaker was placed 1.5 m away from the array, such that the source-receiver axis was normal to the array orientation (see Figure 8(a)). Both source and receiver were 1.2 m above the floor of the chamber which had sections removed to reduce unwanted reflections (see Figure 8(b)). Insertion loss spectra were calculated by subtracting signals received without and with the cylinder array present but with the support structure in place on both occasions.

Concluding remarks

The acoustical properties of three forms of 2D resonators have been investigated: empty N-slit pipes, a concentric arrangement with rigid pipe inner and 4-slit pipe outer and a concentric arrangement of an elastic shell inner and 4-slit outer pipe. A theoretical formulation uses boundary conditions dependent on polar angle to represent arbitrary positioned slits. For simplicity, in the studied geometries the slits are positioned symmetrically. Jump boundary conditions imposed on the slit interface are used to represent the solution inside the slits. The versatility of the proposed method can be applied to various type of 2D resonators with concentric multilayered/solid cylinder inner.

It has been found theoretically and experimentally that increasing the number of slits in an empty pipe causes an increase in the frequency of the Helmhotz-type resonance. A low frequency approximation which models the slit cylinder by an equivalent fluid layer predicts that the frequency increase is proportional to the square root of the number of slits and this has been confirmed experimentally. The concentric arrangements result in resonances associated with both circular and annular cavities. With an inner elastic shell, an additional axisymmetric resonance of the shell is preserved but modified by the presence of the outer 4-slit pipe. Coupling between the components of the concentric arrangement results in shifts in the resonant frequencies corresponding to each element of the composite configuration. This is similar to the effect that is observed in massspring systems with multiple degrees of freedom. A low frequency approximation for the acoustical properties of the concentric arrangement with an elastic shell inner correctly predicts the observed frequency shift in the axisymmetric resonance.

When used in periodic arrays the concentric arrangements with inner elastic shells and outer 4-slit cylinders result in additional sound attenuation in the low-frequency range below the first Bragg band gap while still preserving the Bragg band gaps. This arrangement is more practical that the use of unprotected elastic shells so is potentially useful basis for a sonic crystal barrier design.

Figure 1 .

 1 Figure 1. (a) Cross-section of composite element consisting of a concentric arrangement of an outer 4-slit rigid cylinder and an inner elastic cylindrical shell. (b) Geometry of the slit.

  (a) Internal and external arc lengths of the rigid strip have the same length; (b) Internal and external arcs subtend the same angle. (c) The angle φ n subtended by the arc of n-th slit and its width d n are related by d n = φ n r o , n = 1..N (d) The angle ψ m subtended by the arc between slit (n-1) and slit n and its width D m are related by D m = ψ m r o , m = 1..N

2π 0 <

 0 • > exp(-imθ)dθ we arrive at infinite algebraic system of equations in A n , n ∈ Z, variables as follows ∞ n=-∞

Figure 2 .

 2 Figure 2. Semi-analytical predictions (----) compared with data (---) for single scatterers. Insertion losses are computed at the point (1.5775, 0) m. The rigid shell has outer radius ro = 0.0275 m, thickness h = 0.002 m and N identical slits of width dn ≈ 0.004 m, n = 1..N,. (a) Single slit rigid cylinder (1S); (b) N-slit rigid cylinder with N=4 uniformly distributed slits (4S); (c) Concentric arrangement with outer 4S and inner rigid cylinders with a 1 = 0.011m; (d) Concentric arrangement with outer 4S rigid cylinder and inner latex elastic shell.

Figure 3 .

 3 Figure 3. Solution (9)-(16) (---) for single N-slit scatterer compared with its approximations (----) defined by (29), (30) and (33). Insertion loss is computed at the point (1.5775, 0)m. (a) 1S; (b) 4S; (c) Concentric arrangement with outer 4S and inner rigid cylinders with a 1 = 0.011m; (d) Concentric arrangement with outer 4S rigid cylinder and inner latex elastic shell.

  Figure 3. Solution (9)-(16) (---) for single N-slit scatterer compared with its approximations (----) defined by (29), (30) and (33). Insertion loss is computed at the point (1.5775, 0)m. (a) 1S; (b) 4S; (c) Concentric arrangement with outer 4S and inner rigid cylinders with a 1 = 0.011m; (d) Concentric arrangement with outer 4S rigid cylinder and inner latex elastic shell.

  ) and (b) which show a shift of Helmholtz resonance frequency by a factor of 2 as N increases from 1 to 4.

Figure 4 .

 4 Figure 4. Semi-analytical predictions (----) compared with data (---) for 7 × 3 arrays of scatterers. Distances from array to the source and receiver are 1.5 m and 0.05 m, respectively. The scatterers are arranged in square lattice with lattice constant L = 0.086 m. (a) Empty 4-slit rigid cylinders (4S); (b) Concentric outer 4S cylinder and inner latex elastic shell.

Figure 5 .

 5 Figure 5. Semi-analytical predictions for 7 × 3 array of composites (---) compared with their low-frequency approximations (----). Distances from scatterer to the source and receiver are 1.5 m and 0.05 m, respectively. (a) 4S; (b) Concentric outer 4S rigid cylinder and inner latex elastic shell.

Figure 6 .

 6 Figure 6. Band diagram for the infinite periodic array of composites (left inset) compared with the low-frequency approximations of their finite 7 × 3 arrangement (right inset). All geometrical parameters are identical to those used in Figure 5. The enclosed shells in the infinite periodic structure are taken without viscous effects and their Young's modulus is identical to that used in [16], Figure 6. (a) 4S; (b) Concentric outer 4S rigid cylinder and inner latex elastic shell.

  ) and (b) show an example measurement arrangement. Supports

Figure 7 .Figure 8 .

 78 Figure 7. (a) Components of the resonating scatterer: (on the left) photograph of PVC pipe showing the discontinuous form of a (single) slit and (on the right) photograph of latex elastic shell. (b) Cross section of the concentric arrangement of inner latex cylinder and outer PVC cylinder with 4 symmetrically-placed slits

Acknowledgments

Acknowledgments This work was supported by the EPSRC research grants EP/E063136/1 and EP/E062806/1. Authors are grateful for this support.