Kevin Kendall

Aman Dhir

Chin Yong

Strength by Atomic Force Microscopy (A Kelly special issue)

Keywords: atomic force microscopy, fracture, indentation testing, molecular dynamic simulations AFM, molecular dynamics model, MgO, water squeezing

Localised strength testing of materials is often carried out in the Atomic Force Microscope (AFM) as foreseen by Kelly in his book Strong Solids. It is known that contamination during AFM indentation experiments can strongly influence the observed strength, but a major problem is the theoretical interpretation of the results. Here we use molecular dynamics computer modelling to describe the contact of NaCl and MgO crystal probes onto surfaces, comparable to an AFM experiment. Clean NaCl gave elastic, brittle behaviour in contact simulations at 300K whereas MgO was more plastic, leading to increased toughness. This paper goes on to consider the strength of an oxide substrate contaminated by water molecules tested by indentation with a pyramidal probe of oxide crystal. Recent theory on the effect of liquid contaminant layers on surface strength has been mainly focussed on Lennard Jones (LJ) molecules with some studies on alcohols and water, described by molecular dynamics, which allows the molecules to be squeezed out as the crystal lattice is deformed. In this work we have focused on water by studying the forces between a magnesium oxide (MgO) atomic force microscope (AFM) probe and an MgO slab. Force versus separation has been plotted as the AFM probed was moved towards and away from the substrate.

Simulation results showed that the water layers could be removed in steps, giving up to four force peaks. The last monolayer of water could not be squeezed out even at pressures where MgO deformed plastically. Interestingly, with water present, the strength was reduced but more in tensile

Introduction

In the book Strong Solids, Kelly and Macmillan [START_REF] Kelly | Strong Solids, 3 rd Edition[END_REF] discussed the highest stress that could be withstood by a perfect material such as a single crystal of NaCl or MgO following the original calculations of Zwicky [START_REF] Zwicky | [END_REF] and described the idea that a nanometre scale indenter would be likely to locate a region that was locally free of defects such as cracks or dislocations, leading to the idea of AFM strength testing of perfect crystals. Bowden and Gane [3] had pioneered such experiments using flat ended indenters inside a scanning electron microscope (SEM). Further work on spherical indentations in materials such as {100} MgO by Macmillan and Gane [4] showed that there was a sudden plastic indentation at a high load as the tip broke through the contaminant layer which is omnipresent inside SEM chambers [5,6]. They used the Hertz theory of elastic contact to interpret their results whereas Kelly had been interested in using molecular models to derive the maximum strength using a stress criterion of failure. Experimental AFM tips are normally supplied with a nominal tip radius down to 2nm, suggesting that the atomic facet making contact could be about 4x4 atoms, so that Hertz theory of contact is not applicable. Moreover, a dynamic atomic theory of the contamination effect did not exist at that time so it was not clear whether impurities on the surface should strengthen or weaken the material.

Recently there has been much interest in the molecular interpretation of AFM experiments, especially using molecular dynamics computer simulations to provide a more satisfactory theory of strength than the models given in chapter 1 of Strong Solids. For example, a flat-ended MgO probe contacting a clean MgO crystal was analysed by Yong, Smith and Kendall [7][8][9] and an NaCl probe contacting an NaCl slab was described by Kendall, Yong and Smith [10] using the DL_POLY code which has been much used in atomic scale simulations [START_REF] Smith | DL_POLY, version 2.12[END_REF]. The purpose of this paper is to review the AFM strength measurement using molecular simulations, then to look at the effects of probe size, shape and plasticity while finally considering the effect of contaminant layers on the crystal surfaces.

Four steps are necessary in the atomic model: a)

defining acceptable interatomic potentials b) adding potentials together for all the atoms in the system, assuming the total energy in the system to be built up from pairwise addition of constituent atom-atom interactions c) calculating thermal motions of the atoms to take the dynamic nature of the system into account.

d)

Working out the equilibrium condition for a range of AFM scenarios, including surface contamination.

Previous theories of solid surface interactions under water have been dominated by DLVO continuum theory (see Ref. [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF]) which deals with extended surfaces at relatively large separation, balancing ionic repulsions between charged double layers against the van der Waals attraction. The aim of this paper is to progress to atomistic modelling, ignoring double layer ionic effects in the first instance.

Model set-up; clean surfaces

The probe and slab set-up is shown in simulate increasing size of the AFM probe. Alternatively, pyramidal probes were constructed with varying contact areas from 2x2 to 10x10 atoms. These probes are small enough to ensure that they do not interact with its image neighbour in the periodic boundary conditions.

The layers of atoms were classified into three different groups as shown in Fig 1 . The first three layers of atoms, labelled free, near the contacting surfaces where commensurate contact was to take place were allowed to move freely to reach their equilibrium positions. No constraint was placed on these outer layers of atoms and this ensured that the natural outcome of the atomic configurations at the surface was entirely due to the interatomic adhesion forces.

However, the fourth and fifth layers of atoms were coupled to the Berendsen et al [12] heat bath to maintain the temperature of the whole system at 300K. The final layer of atoms, labelled rigid, was held fixed and gradually displaced to push the nanoprobe towards the surface. The initial gap between the surface slab and the probe was a vertical distance to the (100) surface plane, d z , of around 1.5nm and this was gradually closed, thereby increasing the probe interaction with the slab. The surface slab and probe were orientated such that they were in commensurate contact. In each step when the probe position is updated, the whole system is allowed to equilibrate before the forces experience by the rigid probe along the z-component, F z , is being sampled. This corresponds to an approach rate of 6. 10 -4 ms -1 for equilibration before the interaction force is calculated. The whole process was repeated for different distances so that a force profile of surface interactions was obtained. The potential chosen for this simulation contained three terms, w = -e 2 /r -C/r 6 + Aexp(-r/ρ)

Results for clean surfaces

where e is the electron charge, r is the distance between atoms and C relates to the secondneighbour van der Waals interaction. The parameters A and ρ relate to the repulsive ionic core interactions. These constants are fitted empirically to the experimentally determined properties of NaCl such as elastic constants and crystal lattice parameters [13]. The first term is the electrostatic attraction between Na and Cl ions, the second the van der Waals attraction and the third term is the repulsion between the ions, the last two terms comprising the well-known Buckingham potential.

For distribution equivalent to a temperature of 300K. The system was then allowed to equilibrate until a stable mean configurational energy was achieved, usually taking about 80 ps.

A typical force versus distance curve for the 8x8 NaCl probe is given in Fig 3 . The probe movement started on the right as shown by the bold line and it was clear that there was near zero force of attraction at this stage. The probe was then moved towards the slab until an attractive force was experienced.

Then a sharp change in F z was observed, the 'jump to contact'. Further displacement caused the force to approach zero, the equilibrium contact position. Subsequently, the probe was compressed as it was pushed into the slab giving elastic indentation to the left of Fig 3 . The withdrawal was essentially reversible but the 'jump from contact' or fracture of the bond was slightly delayed, indicating cohesion hysteresis, although the strength was essentially the same.

This indicates that strength should be calculated equally on both make and break of the bonds.

By changing the size of the nanoprobe, the average work of cohesion W in joules per square metre could be calculated from the model and compared with the breaking stress (ie pull-off force /contact area) in GPa as shown in Table 1. It can be seen that the work of cohesion remained almost constant, as assumed in the continuum theory of fracture [START_REF] Kendall | Molecular adhesion and its applications[END_REF] but the stress decreased as the AFM probe got bigger as expected from fracture mechanics theory. In other words, the stress criterion of fracture assumed by Kelly does not apply in this model. In addition, it may be seen that the energy is dominated by the electrostatic forces for NaCl; a van der Waals probe contact would have W around 0.1Jm -2 . The maximum strength was calculated to be 4.03 GPa compared to 2.4 GPa calculated in the book Strong Solids.

Elastic and plastic deformation

The 8x8 rectangular NaCl nanoprobe behaved elastically and separated cleanly despite the large attractive interatomic forces. More commonly, for MgO and also for pyramidal probes, the forces rose so high as the separation occurred that the atoms were displaced from their normal positions and irreversible plastic deformation took place. increased contact sizes of NaCl pyramidal probes [8], it was shown that small contacts eg 2x2

were brittle but larger contacts eg 8x8 became plastic as the nanoprobe was pulled off.

Plastic deformation could also be induced by preloading the probe in compression such that plastic collapse occurred at high force. This introduced defects into the probe, causing subsequent plastic flow when the probe was withdrawn [8].

Influence of plastic flow on strength

The influence of the plastic flow on the observed strength measurement could be studied by comparing the reversibility of the force curves during jump to contact then pull-off. (bottom). The jump to contact during approach caused no damage to the lattice structure and the strength calculated by dividing the total z-component force in the base of the probe by the area of the contact was 3.6 GPa. This was somewhat higher for the pyramidal probe than that calculated for a flat 4x4 NaCl probe seen in Table 1 which gave a strength of 3.1 GPa.

The most significant observation from Fig 5 was that the compressive collapse strength of the NaCl probe, 19.3 GPa, was very much higher than the tensile strength of 3.6 GPa. This is more like a hardness measurement in which plastic flow is required to produce an indentation, typically requiring three times more pressure that tensile yielding. Deleted: Force per ion vs displacement curve for pyramidal 4x4 NaCl probe. The jump to contact appears perfectly elastic but the pull-off dissipates energy without raising the force significantly. The stress required for compressive collapse was much higher. MgO showed very similar behaviour but at higher stresses and with more plastic deformation during pull-off. The 6x6 pyramidal MgO probe gave a jump-to-contact strength of 14.4 GPa without any lattice damage. On pull-off, the same stress was observed, demonstrating reversibility. But then substantial plastic deformation occurred and the jump-from-contact occurred in a series of steps with large dissipation (Fig 6) indicating increased toughness.

Compression of the probe into the slab eventually caused plastic collapse at a stress of 65 GPa, showing that the compressive strength was about five times higher than the tensile strength.

In conclusion, it is evident that a number of complex processes can influence nanoprobe strength including the nature of the crystalline material, the atomic contact, its loading history, the shape of the probe, tension or compression and the propensity for plastic flow.

Influence of surface contamination

The cohesive forces for clean NaCl surfaces were seen to be very large, leading to plastic effects and irreversibility, with adhesive energies around 0.47 Jm -2 . Of course, this can only be true for atomically clean surfaces, as tested in high vacuum experiments with baked out materials. In practical AFM testing, ubiquitous contamination such as nitrogen or water adsorbs on the surface and this has a very large effect of reducing the adhesive energy to much smaller values around 0.1 Jm -2 or even lower if thicker layers of contamination build up leading to energy and force oscillations first proposed by Kendall [START_REF] Kendall | [END_REF]. For example, Fig 7 shows a computational model of a pyramidal 6x6 MgO probe pushed towards an MgO slab Fig 6 Force vs displacement curve for 6x6 pyramidal MgO probe showing that the jump to contact was perfectly elastic. The pull-off was at the same force but dissipated energy plastically. The steps in the breaking curve were due to sudden motion of defects within the nanoprobe, which remained attached to the slab. interacting with the surface atoms with an energy parameter of ε = 0.02 eV and σ = 0.305 nm at a temperature of 100 K. Fig. 7c shows the response of the MgO probe in the presence of the contaminant, with the contact force increasing from right to left as the gap was squeezed. The contamination prevented the jump-to-contact in this case but there were four significant steps; first at 0.02nN where a slight bump could be seen, second at 0.1 nN as the third layer of xenon was penetrated at 3.5nm; next at 0.35 nN as the second xenon layer was squeezed out at 3.2nm; and finally at 0.9 nN as the probe collapsed plastically at high force because the first xenon layer could not be squeezed out by the probe.

Similar calculations had been carried out previously by Gao et al [16] but their probe shown in Fig 8 had not been suitable for experimental study because of its infinite length . However, they did show that the structuring phenomenon was moderated by the alignment of the solid surfaces, with commensurate crystals showing strong liquid structuring in the model. [16] showing the flat parallel plate configuration with liquid in between; the solid surfaces were infinite in the y direction. That simulation had shown squeezing out of the contaminant molecules in discrete steps, rather like the experimental results found in mica experiments by Horn, Pashley et al [17][START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF][START_REF] Pashley | [END_REF]. AFM experiments have not found this effect consistently, although Lim et al [20] showed stepwise squeezing out of organic molecules such as n-dodecanol from a graphite substrate, depending strongly on the size of the AFM tip, with finer tips (R tip =25nm) giving reliable force measurement.

They detected about 6 layers of molecules squeezing step-wise from the graphite, and also suggested the full tip contact with the graphite could be obtained using Si 3 N 4 AFM probes. Most previous experiments on this effect with water have been carried out with mica surfaces in the presence of ions [21] and this has been modelled by Leng [22] to show the repulsions associated with structured water layers. The purpose of our study was to use an AFM geometry to test the strength of MgO in the presence of water.

Simulation of water between MgO surfaces

To investigate the question of how water could be squeezed out from oxide surfaces in AFM strength testing, the model illustrated in The calculations showed that the jump to contact forces (Fig 9b ) were much reduced by the presence of the water molecules. As the probe moved towards the slab, some water molecules jumped across the gap to cover the probe surface in a monolayer, causing a dip in the approach force curve as the gap between the MgO surfaces was fully filled with water. Then one layer of water was squeezed out as shown by the increase in force followed by a drop. Two more layers were then squeezed out at increasing force, but the final layer could not be removed in this model.

On withdrawal, the force curves were not reversible but showed some hysteresis as the water penetrated back into the gap between probe and slab. It was clear that the MgO probe could now be removed from the MgO slab without any plastic deformation, which was impossible for the clean surfaces. The maximum tensile stress at pull-off was 1.57 GPa compared to 14.4 GPa for clean surfaces, illustrating the factor of nine weakening effect of water in tension.

A key question was whether the strength of the MgO was increased by the surface contaminant, as suggested in Strong Solids. To test this, the force on the probe was increased until the MgO collapsed plastically. With water present, the collapse stress for the 6x6 MgO probe was 9.6 GPa compared to 65 GPa for the clean probe, showing that the contamination had given a factor of seven weakening in compression.

One of the most interesting features of this model is that there are now several zeros in the force versus separation curve, depending on the number of water molecule layers separating the surfaces, intimating that the probe could appear to be in equilibrium at several separation values depending on the applied force. There were three discrete equilibria at gaps of 0.5 nm, 0.75 nm and 1nm. This raises the question of colloidal stability of MgO particles in water, where the kT of thermal energy of the particles will push the contacts into one of the energy wells by squeezing out one or more water layers [24].

Conclusions

Strength of perfect crystals can be measured in the Atomic Force Microscope if the theory of crystal fracture and plasticity is understood. Molecular dynamics models of NaCl and MgO probes pushed into perfect crystal substrates have been used to interpret the processes at the nanometer level for comparison with AFM tests.

Clean surfaces of NaCl or MgO attract each other and pyramidal probes give the expected 'jump to contact' followed by continued approach to equilibrium at zero force. This process was found to be elastic and reversible with NaCl flat probes, showing that the nanocrystals could be pulled apart without damage and that the strength of 4.03 GPa was reversible, ie equal on approach and pull-off, though slight 'cohesion hysteresis' was observed. NaCl followed an energy criterion of fracture in this case, showing essentially brittle fracture, with strength diminishing for larger probes as expected from fracture mechanics. Plastic flow of the NaCl probes could be encouraged by several means:-first by enlarging the contact region, second by changing the probe shape and thirdly by damaging the probe. MgO showed the same effects but differed from NaCl in that plastic flow always dominated the behaviour. Another measure of MgO strength was the probe collapse which was observed when increased compression of the probe onto the slab were made. Plastic collapse of the probe took place at a well defined compressive stress of 65 GPa, larger than the tensile strength of 14.4 GPa measured at jump to contact and pull-off.

When the MgO surfaces were contaminated with water molecules, there was a very significant change in behaviour. First, the jump to contact was now periodic as layers of molecules were squeezed out from the gap between probe and slab. Second, the process of detachment was now easier and the MgO probe could be detached at low tensile stress of 1.57 GPa without damage.

Third, the strength of the MgO probe as it collapsed on severe compression into the slab was altered. The compressive strength was diminished by a factor 7 in the presence of the water molecules.

Fig 1 .

 1 The slab consisted of 28x28 rows of either NaCl or MgO lattice. The Na + -Cl -inter-ionic spacing was 0.279 nm, while for MgO, the Mg 2+ -O 2- inter-ionic spacing was 0.2104 nm. Periodic boundary conditions were applied in the x-y (horizontal) direction to mimic an infinite crystal surface with the (100) surface oriented at the z (vertical) direction. A small rectangular probe of 6 layers thick was constructed and this was varied in cross-section from the smallest possible 2x2 shown in Fig 1 to 18x18 atoms to Formatted: Indent: Left: 36 pt Formatted: Font: Times, Superscript Formatted: Font: Symbol Formatted: Font: Times, Deleted: DL Deleted: molecular dynamics theory of Deleted: an Deleted: which is acceptable Deleted: NaCl Deleted: and the probe was 6 layers deep.

Fig 2 ,Fig 2

 22 Fig 2 Molecular dynamics model of a sodium chloride 8x8 atom probe approaching a slab of NaCl to form new bonds with the surface atoms

  Fig 4 shows the computed results for a stepped NaCl nanoprobe of contact size 8x8 in contrast to a stepped probe with contact size of 4x4. As stress was applied to the 8x8 contact (Fig 4 top), separation did not occur, but the step acted as a defect which promoted plastic deformation, allowing the NaCl crystal probe to extend in line with the applied force. The smaller contact (Fig 4 bottom) alsobegan to deform plastically, but separation then took place at the interface such that the nanoprobe was released but extra energy was expended in creating a plastic defect hole in the probe. This model shows that there should be an elastic/plastic transition for larger nanoprobes, such that material transfer from the probe to the slab occurred. By modelling

Fig 5 shows

 5 the force distance curves computed for the pyramidal NaCl probe shown in Fig4

Fig 4 (Fig 5

 45 Fig 4 (top) Major plastic deformation observed in the 8x8 contact model of the pyramidal NaCl nanoprobe. (bottom) Results for 4x4 NaCl contact show clean detachment with minor plastic deformation

Formatted

  covered by inert gas particle spheres with xenon-like van der Waals attributes at 100K. The 6 × 6 MgO probe was surrounded by the Lennard Jones (LJ) xenon-like spheres and it was immediately evident that these xenon molecules formed structured layers around the probe and upon the slab surface, with 3 distinct layers easily visible (Fig 7 a & b). Each sphere was weakly

Fig 7 aFig 8

 78 Fig 7 a),b) Stills from a movie of molecular dynamics model of 6x6 MgO pyramidal probe approaching an MgO slab coated with a monolayer of xenon molecules at 100K; c)force distance curve approaching from right to left, showing steps as molecules are squeezed out.

  Fig 9 was set up. The baseline comparison was with clean MgO which gave the results shown in Fig 6. The calculations for clean MgO were compared to the model with 3 layers of water molecules positioned on the slab as the probe was moved downwards (Fig 9a). The potential for water interacting with MgO was derived from de Leeuw and Parker [23].

Formatted

Fig 9a )

 9a Fig 9a) Pyramidal 6x6 MgO probe moving towards an MgO slab covered with 3 layers of water molecules; b) force distance curve revealing the filling of the gap then the oscillations as three layers of water were squeezed out, leaving one layer resisting removal, showing hysteresis

Fig 1

 1 Figs for 'Strength by Atomic Force Microscopy (AFM): Molecular dynamics of

Fig 2 Fig 3 Fig 5

 235 Fig 2 Molecular dynamics model of a sodium chloride 8x8 atom probe approaching a slab of NaCl to form new bonds with the surface atoms

Fig 6 Fig 8

 68 Fig 6 Force per 10x10 ion vs displacement curve for 6x6 pyramidal MgO probe showing that the jump to contact was perfectly elastic. The pull-off was at the same force but dissipated energy plastically. The steps in the breaking curve were due to sudden motion of defects within the nanoprobe, which remained attached to the slab.

  Fig 1 Schematic of the smallest 2x2 probe approaching a crystal slab. The probe was moved towards the slab by displacing the probe rigid layer downwards, keeping the slab rigid atoms fixed. The free atoms could then find their equilibrium positions while the bath atoms maintained the temperature at 300K.

	In several papers [7-9] the equilibrium was calculated for each incremental step as the contact
									Deleted: , small enough to avoid
	process occurred. A still picture from the movie for NaCl just before the point of complete	boundary effects which caused errors with large contacts eg 15x15.
	fracture is shown in Fig 2, but this does not reflect the true situation because all the atoms are	Formatted: Font color: Blue
	in	continuous	motion.	A	movie	is	available	at
	http://www.fuelcells.bham.ac.uk/media/nacl-md.shtml showing all the atoms
	vibrating around their equilibrium positions as the nanoprobe approaches and then retracts
	from the NaCl slab.
			F o r
			P
			e
				e r
					R				Formatted: Font: Times, Subscript
					e
					v i e		Formatted: Font: Times, Subscript
							w		Formatted: Superscript Deleted: At Deleted: ium
							O n l	Formatted: Superscript Deleted: was attained at each step and
								y

Table 1

 1 Work of cohesion W and average breaking stress for the different size NaCl probes

Contact size area/nm 2 W/ Jm -2 stress/GPa

  Force per ion versus distance curves for approach from right to left and withdrawal of the 8x8 NaCl nanoprobe showing the jump to contact and the equilibrium at zero force

	Deleted: Fig 3 Deleted: i r o F (8 × 8) 4.9818 0.45 2.42 (10 × 10) 7.7841 0.45 2.38 (14 × 14) 15.257 0.44 2.08 (18 × 18) 25.220 0.43 2.03
			P		Deleted: Deleted: article
			e		Deleted: the particle
			e r
				R
				e
				v i e
				w
	(2 × 2) (4 × 4)	0.3114 1.2455	0.47 0.48	4.03 3.06	O n l
	(6 × 6)	2.8023	0.44	2.55	y

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


Katarzyna Filipiak

Rafal Różański

Aneta Sawikowska

Dominika

Rafał Różański

Dominika Wojtera-Tyrakowska

Rafa L Różański

On the E-optimality of complete designs under an interference model

Keywords: Interference model, Information matrix, E-optimal designs, Spectral norm, Reducible matrices

In the experiments in which the response to a treatment can be affected by other treatments, the interference model with neighbor effects is usually applied. It is known, that circular neighbor balanced designs (CNBD) are universally optimal under this model. However, such designs cannot exist for each combination of design parameters. If in the class of block designs with a given number of blocks there exists CNBD, then such a design cannot exist in the class with one block less or more. The aim of this paper is to identify the E-optimal complete block designs under the interference model with neighbor effects over such classes of designs. The results are based on algebraic properties of information matrix.

Introduction

In the theory of experimental designs, the problem of determining optimal designs is often considered. If in an experiment the response to a treatment is affected by other treatments (for example in agricultural and horticultural experiments), then the optimality of designs under an interference model is studied. Recently some results on universal optimality of binary designs under this model were published. These results concern mainly optimality of circular neighbor balanced designs (CNBD) and orthogonal arrays of type I under the fixed and mixed interference models, where the observations are correlated or not (see e.g. [START_REF] Druilhet | Optimality of circular neighbour balanced designs[END_REF][START_REF] Filipiak | Optimality of circular neighbor balanced designs under mixed effects model[END_REF], 2004, 2005[START_REF] Filipiak | Optimal designs for a mixed interference model[END_REF]. It is known, however, that for some combinations of design parameters the universally optimal designs cannot exist. In such a case efficiency of some designs or optimality with respect to the specified criteria are considered. The efficiency of some cyclic designs under a fixed interference model was studied in [START_REF] Filipiak | Some properties of cyclic designs under an interference model[END_REF] In this paper we are interested in determining E-optimal designs under an interference model. This criterion is based on the eigenvalues of the information matrices of designs. In the literature this problem has not been discussed yet. Only the problem of determining E-optimal designs under the model with block effects as the only nuisance parameters over the class of designs is partially solved, see e.g. [START_REF] Constantine | Some E-optimal block designs[END_REF][START_REF] Jacroux | Some E-optimal designs for the one-way and two-way elimination of heterogeneity[END_REF][START_REF] Jacroux | On the E-optimality of block designs[END_REF][START_REF] Srivastav | On the E-optimality of certain class of block designs[END_REF]. However, there is no characterization of E-optimal designs under an interference model.

It is known, that in the class of block designs with the same number of treatments, t, as block sizes, k, CNBD cannot exist for the number of blocks b = t -2 and b = t. Thus, for such design parameters we determine the structure of left-neighboring matrix of E-optimal design under the interference model with left-neighbor effects and we give some examples of such designs.

It is worth observing that designs with t = k are often used in practice. For example in UPOV (The International Union for the Protection of New Varieties of Plants) research, complete block designs are recommended in experiments when the number of treatments is less than 16. The designs with the same number of blocks as number of treatments and units are also applied in clinical trials. This paper is organized as follows. First we present some general definitions and notation. In Section 3 we characterize the structure of the left-neighboring matrix of E-optimal design with t = k = b. In Section 4 we determine the structure of the left-neighboring matrix of E-optimal design with t = k = b + 2. In the last section we give the method of construction of E-optimal designs and some examples.

Let D t,b,k be the set of designs with t treatments, b blocks and k experimental units per block. An interference model with left-neighbor effects associated with the design d ∈ D t,b,k can be written as

y = T d τ + L d λ + (I b ⊗ 1 k )β + ε, (1)
where τ , λ and β are the vectors of treatment effects, left-neighbor effects and block effects, respectively. Here ε is a vector of random errors with E(ε) = 0 and Cov(ε) = σ 2 I bk , where σ 2 is an unknown constant. The matrix I b denotes the identity matrix of order b, 1 k is the k-vector of ones and ⊗ denotes the Kronecker product. Let T du be the design matrix of treatment effects in block u, 1 ≤ u ≤ b. Further, define T d = (T d1 : • • • : T db ) as the design matrix of treatment effects. For each u we define

L du = H k T du , where H k is a k × k matrix of the form H k = 0 k-1 1 I k-1 0 k-1 . (2)
Then, L d = (I b ⊗ H k )T d is the design matrix of left-neighbor effects. This form of the matrix H k follows from the assumption that each treatment has a left neighbor. This situation may occur if each block of a design has the form of a circle. If plots in blocks are arranged in linear forms, we can obtain the effect of circularity by adding border plots at the beginning of each block, where the treatment at the border plot is the same as the treatment at the opposite end of the block (for more details see e.g. Druilhet, 1999, Filipiak and[START_REF] Filipiak | E-optimal designs under an interference model[END_REF]. Border plots are not used for measuring the response variables. Under the interference model (1), the information matrix for estimating the treatment effects, C d , can be expressed as

C d = T d Q I b ⊗1 k T d -T d Q I b ⊗1 k L d (L d Q I b ⊗1 k L d ) -L d Q I b ⊗1 k T d , (3)
where Q I b ⊗1 k is the orthogonal projector onto the orthocomplement of the column span of

I b ⊗ 1 k and A -denotes a generalized inverse of A. It is easy to see that Q I b ⊗1 k = I b ⊗ E k , where E k = I k -1 k 1 k 1 k . Because of the form of E k and since H k is orthogonal, it can be seen that T d Q I b ⊗1 k T d = L d Q I b ⊗1 k L d .
We are interested in determining E-optimal design, i.e. such design, that the maximal variance among all best linear unbiased estimators of normalized linear contrasts is minimal under this design. The condition of E-optimality can be expressed in terms of eigenvalues of the information matrix as follows. For a design [START_REF] Constantine | Some E-optimal block designs[END_REF].

d ∈ D t,b,k let λ max (C d ) = λ 1 (C d ) ≥ λ 2 (C d ) ≥ • • • ≥ λ t-1 (C d ) ≥ λ t (C d ) = 0 be the eigenvalues of its information matrix C d . A design d * ∈ D t,b,k is called E-optimal over D t,b,k if λ t-1 (C d * ) ≥ λ t-1 (C d ) for all designs d ∈ D t,b,k (cf.
Observe, that in formula (3) it is involved a generalized inverse of a matrix, which depends on a design, i.e. which changes itself with the arrangement of treatments on experimental units. It makes the determination of E-optimal design difficult. Therefore in this paper we consider experiments with t = k (t > 2). The class of designs with t = k, in which each treatment occurs at most once in each block (the class of complete binary designs) we will denote by B t,b,t . In such a case the information matrix has the form

C d = bE t - 1 b K d K d ,
where

K d = T d Q I b ⊗1t L d .
Since the vectors T d 1 t and L d 1 t are in the column space of I b ⊗ 1 t , the matrix K d has zero row and column sums, i.e.

K d 1 t = K d 1 t = 0 t
, where 0 t is a t-vector of zeros.

Observe that the matrix E t is symmetric, idempotent and t -1 of its eigenvalues are equal to 1. Moreover, E t 1 t = 0 t and K d K d 1 t = 0 t . Thus, under the class B t,b,t , the E-optimality condition can be rewritten as follows

λ 1 (K d * K d * ) ≤ λ 1 (K d K d ) for all d ∈ B t,b,t . (4)
Since for a nonnegative-definite matrix A, λ 1 (AA ) = σ 1 (A), where σ 1 (A) denotes the spectral norm of A, to determine E-optimal design d * we have to find the matrix K d * with minimal spectral norm over matrices K d , d ∈ B t,b,t . We use some special properties of K d .

Using the form of Q I b ⊗1t , the matrix K d can be expressed as

T d L d -b t 1 t 1 t .
The (i, j)th entry of T d L d denotes the number of occurrences of treatment i with treatment j as left neighbor in a design d ∈ B t,b,t . Therefore we will call this matrix a left-neighboring matrix and we will denote it by S d with entries s ij , i, j = 1, 2, . . . , t. It is easy to observe, that in the class B t,b,t , the diagonal entries of S d are equal to 0 and the off-diagonal entries belong to the set {0, 1, . . . , b}. Moreover, S d 1 t = S d 1 t = b1 t . Thus, the set of all matrices K d for d ∈ B t,b,t can be defined as

K b = K d = (k ij ) 1≤i,j≤t : K d 1 t = K d 1 t = 0 t , k ij ∈ - b t , 1 - b t , . . . , b - b t , k ii = - b t .
Let P t be the set of all permutation matrices of order t and let P t ⊂ P t be the set of permutation matrices with zero diagonal (the set of derangement matrices). We show that the matrix K d * with minimal spectral norm over K b belongs to the set

K t = K d ∈ K t : K d = P d -I t , P d ∈ P t if b = t and K t-2 = K d ∈ K t-2 : K d = 2 t 1 t 1 t -I t -P d , P d ∈ P t if b = t -2.
We prove that design d * which is E-optimal over B t,b,t , i.e. design for which K d * has minimal spectral norm over K b , is also E-optimal over D t,b,t .

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

In the paper reducibility of matrices, will play an important role. Reducibility is defined as follows.

Definition 2.1 [START_REF] Horn | Matrix Analysis[END_REF]. An n × n matrix A is said to be reducible if either of the following conditions is satisfied 1. n = 1 and A = 0, 2. n ≥ 2 and there is a permutation matrix P ∈ P n and an integer r with 1 ≤ r ≤ n -1 r) and

such that P AP = B C Θ D , where B ∈ C r×r , D ∈ C (n-r)×(n-r) , C ∈ C r×(n-
Θ ∈ C (n-r
)×r is a zero matrix.

A matrix is called irreducible if it is not reducible.

The matrix P AP, where P ∈ P n , is called the matrix permutationally similar to A. Recall, that the eigenvalues of A and a matrix permutationally similar to A are the same.

Observe, that every irreducible derangement matrix is permutationally similar to the matrix H t , defined in (2). Moreover, all reducible derangement matrices are permutationally similar to the block diagonal matrices with derangement matrices on the diagonal. It is worth to note, that the eigenvalues of each matrix permutationally similar to H t are equal to ω k , k = 0, 1, . . . , t -1, where ω k are tth roots of unity. It follows from circularity of H t (for more details see e.g. Horn andJohnson, 1985, John, 1987).

In this paper we consider b = t -2 and b = t, because for the remaining numbers of blocks determining the spectral norm of K d requires the spectral norm of two (or more) different permutation matrices, which is very difficult to solve in general.

E-optimal designs over D t,t,t

In this section we determine the structure of the left-neighboring matrix of E-optimal design over D t,t,t . First we determine, which structure of K d from K t gives the minimal spectral norm over K t . Then we show the structure of the left-neighboring matrix of E-optimal design over the class of binary designs B t,t,t and, finally, over the class of all designs D t,t,t .

Observe that a matrix K d ∈ K t is reducible if and only if it is permutationally similar to the block-diagonal matrix diag(H t 1 , H t 2 , . . . , H tn ) -I t , where n j=1 t j = t and n = 1. It is obvious that for n = 1 we obtain an irreducible matrix.

Since the sets of eigenvalues of permutationally similar matrices are the same, the problem of finding the matrix from K t with minimal spectral norm can be reduced to studying the spectral norms of H t -I t or diag(H t 1 , H t 2 , . . . , H tn ) -I t .

Consider a matrix [START_REF] John | Cyclic Designs[END_REF]. It is easy to see that λ max (K d K d ) is obtained for k = t 2 for even t and k = t±1 2 for odd t and

K d = H t -I t . Then K d K d = 2I t -(H t + H t ) = 2I t -(H t + H -1 t ) is circular and its eigenvalues are equal to 2 -ω k -ω -1 k = 2 -2 cos 2kπ t , k = 0, 1, . . . , t -1 A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT (cf.
λ max (K d K d ) = 2 -2 cos π = 4 for even t 2 -2 cos (t±1)π t = 2 + 2 cos π t for odd t.
(5)

Moreover, it is easy to see that for odd t

lim t→∞ cos π t = 1, (6)
and

λ max (K d K d ) → 4 -if t → ∞. Since σ 1 (K d ) = λ max (K d K d ), for K d = H t -I t we get σ 1 (K d ) = 2 for even t and σ 1 (K d ) → 2 -if t is odd and t → ∞.
Recall, that the spectrum of the block-diagonal matrix is the sum of the spectra of diagonal blocks (see e.g. [START_REF] Marcus | A Survey of Matrix Theory and Matrix Inequalities[END_REF]. From the block-diagonal form of the reducible matrices from K t it is easy to see that the matrix K d K d is also block-diagonal and we can prove the following theorem.

Theorem 3.1 The matrix K d * which minimizes the spectral norm over K t is permutationally similar to the matrix:

(i) H t -I t if t = 2, 7; (ii) I 2 ⊗ H 2 -I 4 or H 4 -I 4 if t = 4; (iii) I m ⊗ H 3 -I t if t = 3m, m ∈ N;
(iv) diag(I i ⊗ H 3 , I j ⊗ H 5 ) -I t if t = 5 or t ≥ 8 and t = 3m, m ∈ N, where t = 3i + 5j for some i ∈ N ∪ {0}, j ∈ N.

Moreover, in case (iii) there is σ 1 (K d * ) = √ 3 and the spectrum of K d * K d * has the form {0, 0, . . . , 0 m times , 3, 3, . . . , 3 2m times }, and in case (iv) the spectral norm of 5) it is easy to see that the matrix with minimal spectral norm, K d * , has to be permutationally similar to the block-diagonal matrix with diagonal blocks of odd order as small as possible. Hence we obtain immediately case (iii) and respective form of K d * for t = 7. For t = 2 and t = 4 such a situation is impossible. Thus, we obtain respective forms of K d * . For case (iv) it is enough to observe that considered t's can be expressed as 3i + 5j, where i ∈ N ∪ {0}, j ∈ N. Again from ( 5) and ( 6) condition ( 4) is satisfied for

K d * is σ 1 (K d * ) = 5+ √ 5 2 . Proof. Let K d ∈ K t . From (
K d * in a given form. Let K d ∈ K t \ K t .
It is known that the spectral norm of the nonnegative-definite matrix is not smaller than the maximal diagonal entry of this matrix (see e.g. [START_REF] Marshall | Inequalities: Theory of Majorization and Its Application[END_REF]. We will show, that the maximal diagonal entry is not smaller than the upper bound of the spectral norm of the matrix from K t (obtained from the triangle inequality). For simplicity we denote the entries of K d K d by k ij , i, j = 1, 2, . . . , t. We obtain

λ 1 (K d K d ) ≥ max i k ii
and it can be observed that max

1≤i≤t k ii ≥ 4. Hence σ 1 (K d ) ≥ 2. Observing σ 1 (K d * ) = σ 1 (P d * -I t ) ≤ σ 1 (P d * ) + σ 1 (I t ) = 2 we obtain the thesis.
The second part of the theorem is obtained directly from the spectral norm of (H 3 -I 3 )(H 3 -I 3 ) and H 5 -I 5 , respectively.

According to the Theorem 3.1 we can formulate the following corollary, which gives the structure of the left-neighboring matrix of E-optimal design. Corollary 3.1 If there exists design d * with the left-neighboring matrix S d * permutationally similar to the matrix

(i) 1 t 1 t -I t + H t for t = 2, 7; (ii) 1 4 1 4 -I 4 + I 2 ⊗ H 2 for t = 4; (iii) 1 t 1 t -I t + I m ⊗ H 3 for t = 3m, m ∈ N; (iv) 1 t 1 t -I t + diag(I i ⊗ H 3 , I j ⊗ H 5 ), for t = 5 or t ≥ 8 and t = 3m, m ∈ N,
up to the order of diagonal blocks, where t = 3i + 5j for some i ∈ N ∪ {0} and j ∈ N,

then d * is E-optimal over B t,t,t .
Note, that in case (ii) of Corollary 3.1 there is only one form of the left-neighboring matrix, while in Theorem 3.1 there are two forms of K d * . It follows from the fact, that the design with the left-neighboring matrix S d * = H 4 -I 4 + 1 4 1 4 does not exist (cf. [START_REF] Azaïs | A catalogue of efficient neighbour-designs with border plots[END_REF]. Now we prove, that the designs with left-neighboring matrices given in Corollary 3.1 are E-optimal over the class of all designs D t,t,t . For d ∈ D t,b,k , the information matrix C d is of the form (3). Since C d is nonnegative-definite, the following inequality is satisfied:

λ t-1 (C d ) ≤ λ t-1 (T d Q I b ⊗1 k T d ). (7)
Observe, that for d ∈ D t,b,t

T d Q It⊗1t T d = L d Q It⊗1t L d = T d T d - 1 t b j=1 r dj r dj = R d - 1 t b j=1 r dj r dj , (8)
where r j is the t-vector of replications of treatments in the jth block and R d is the diagonal matrix with the number of replications of the ith treatment in the design on the diagonal, i.e. R d = diag(r d1 , . . . , r dt ). We can prove the following.

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

Theorem 3.2 If design d * satisfies conditions of Corollary 3.1, then d * is E-optimal over D t,t,t .

Proof. Let d * satisfies conditions of Corollary 3.1 and let d ∈ D t,t,t \B t,t,t . We prove the theorem in two steps. First we show the thesis for equireplicated designs and then for designs which are not equireplicated. (i) If d is equireplicated, then R d = tI t in (8) and

λ t-1 (T d Q It⊗1t T d ) = t - 1 t λ max t j=1
r dj r dj .

Using the fact that the diagonal elements of a nonnegative-definite matrix are majorized by the characteristic roots, we have

t - 1 t λ max t j=1 r dj r dj ≤ t - 1 t max 1≤i≤t t j=1 r 2 dij ,
where r dij is the number of occurrences of the ith treatment in the jth block in the design d. Since design d is not binary, we get

max 1≤i≤t t j=1 r 2 dij > 4.
According to ( 7) and ( 5) we obtain 8) is not greater than t -1. From [START_REF] Constantine | Some E-optimal block designs[END_REF]

λ t-1 (C d ) ≤ t -4 t ≤ t -1 t λ max (K d * K d * ) = λ t-1 (C d * ). (ii) If d is not equireplicated, then minimal diagonal element in R d given in (
λ t-1 (T d Q It⊗1t T d ) ≤ (k -1)t k(t -1) min 1≤i≤t r di , (9)
and thus in our case

λ t-1 (C d ) ≤ min 1≤i≤t r di ≤ t -1 ≤ t - 4 t ≤ λ t-1 (C d * )
for t ≥ 4. For t = 2, 3 the thesis is obtained by direct calculations.

4 E-optimal designs over D t,t-2,t

In this section we determine the structure of the left-neighboring matrix of E-optimal design over D t,t-2,t . Similarly as in previous section first we determine, which structure of K d from K t-2 gives the minimal spectral norm over K t-2 . Then we show the structure of the left-neighboring matrix of E-optimal design over the class of binary designs B t,t-2,t and, finally, over the class of all designs D t,t-2,t .

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

ACCEPTED MANUSCRIPT

µ i (K d * K d * ) = 0 for i = 0 2 + 2 cos 2iπ t for i = 1, 2, . . . , t -1
and we obtain the maximal eigenvalue of

K d * K d * for i = 1 or i -t -1, i.e. λ max (K d * K d * ) = 2 + 2 cos 2π t = 4 cos 2 π t . ( 10
) Let P d = diag(H t 1 , . . . , H tn ), t 1 + • • • + t n = t. Then the unordered, nonzero eigenvalues of K d K d are µ m (K d K d ) = 2 + 2 cos 2iπ t j for m = 1, 2, . . . , t -1, i = 1, 2, . . . , t 1 -1, 0, 1, . . . , t 2 -1, . . . , 0, 1, . . . , t n -1, j = 1, . . . , n,
and we get λ max (K d K d ) for i = 0 and finally

λ max (K d K d ) = 4 > 4 cos 2 π t = λ max (K d * K d * ). Let K d ∈ K t-2 \ K t-2 .
It is known, that the spectral norm of a matrix is not smaller than the spectral norm of any submatrix of this matrix (cf. [START_REF] Horn | Topics in Matrix Analysis[END_REF]. We use a 2 × 2 submatrix of K d K d . Observe, that in the class K t-2 \ K t-2 the matrix K d has at least two entries equals to 2 -t-2 t , and these two entries are in different rows, say in positions (i 1 , j 1 ) and (i 2 , j 2 ), i 1 = i 2 , j 1 = j 2 (it follows from zero row and column sum of K d ). In such a situation two diagonal entries of K d K d are equal to 4 -4 t and the (i 1 , j 2 )th entry of K d K d is equal to -4 t or 1 -4 t or 3 -4 t . Let A be the symmetric submatrix of K d K d of order 2, which contains entries (i 1 , j 1 ) and (i 2 , j 2 ) on the diagonal and (i 1 , j 2 ) as off-diagonal entries. Then

σ 1 (K d ) ≥ λ max (A) =        √ 3 for t = 2, 3, 4, 5 -8 t for t = 5, 6, 7, 2 for t ≥ 8.
It is easy to calculate that the above spectral norm of K d is not smaller than the spectral norm σ 1 (K d * ) = 2 cos π t obtained from (10). The proof is complete.

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

Corollary 4.1 If there exists design d * with the left-neighboring matrix S d * permutationally similar to the matrix 1 t 1 t -I t -H t , then d * is E-optimal over B t,t-2,t .

Now we prove, that the designs with left-neighboring matrices given in the above corollary are E-optimal over the class of all designs D t,t-2,t .

Theorem 4.2 If design d * satisfies conditions of Corollary 4.1, then d * is E-optimal over D t,t-2,t .

Proof. Let d * be an E-optimal design over B t,t-2,t and let d ∈ D t,t-2,t \B t,t-2,t . We prove the theorem in two steps. First we consider equireplicated designs and in second step we study the class of not equireplicated designs. (i) If d is equireplicated, then R d = (t -2)I t in (8) and

λ t-1 (T d Q It⊗1t T d ) = t -2 - 1 t λ max t-2 j=1
r dj r dj .

Using the fact that the diagonal elements of a nonnegative-definite matrix are majorized by the characteristic roots, we get According to ( 7) and (10) we have

λ t-1 (C d ) ≤ t -3 ≤ t -2 - 1 t -2 • 4 cos 2 π t = λ t-1 (C d * ).
(ii) If d is not equireplicated, then minimal diagonal element in R d given in ( 8) is not greater than t -3. Using (9) we obtain

λ t-1 (C d ) ≤ min 1≤i≤t r di ≤ t -3 ≤ t -2 - 1 t -2 • 4 cos 2 π t = λ t-1 (C d * ).

  is the number of occurrences of the ith treatment in the jth block in the design d. Since design d is not binary,

  Theorem 4.1 The matrix K d * which minimizes the spectral norm over K t-2 is permutationally similar to the matrix 2 t 1 t 1 t -I t -H t .Proof.Let K d ∈ K t-2 . Then K d K d = 2I t -4 t 1 t 1 t + (P d + P d ), P d ∈ P t . Observe, that for P d * = H t the unordered eigenvalues of K d * K d * are of the form

Acknowledgements

The authors are thankful to Professor Charles R. Johnson for helpful discussion on the topic of this paper during the conference MAT-TRIAD 2005.

1 Research partially supported by the Interdisciplinary Joint Research Project 2004/2005 of Adam Mickiewicz University and Agricultural University of Poznań.

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

Construction of E-optimal designs

In accordance to Corollary 3.1 and 4.1 we can identify the left-neighboring matrix S d * of an E-optimal design. In this section we give the method of construction of E-optimal designs over D t,t,t and D t,t-2,t .

From the form of the matrix T d we can observe that each block of complete binary design can be written as a one-cycle permutation matrix. Thus, if we decompose the left-neighboring matrix to the sum of b one-cycle permutation matrices, we can construct a design. Recall, that each permutation is a cycle or a product of disjoint cycles (see e.g. [START_REF] Birkhoff | A Survey of Modern Algebra[END_REF]. For example, the matrix H 3 is a one-cycle permutation matrix, which presents the cycle: 1 → 2, 2 → 3, 3 → 1 (indexes of elements of H 3 that equal 1). Observe that, if we write this cycle in the form [123], we can regard it as a treatment sequence in a single block of design.

Consider the class D t,t,t . From Corollary 3.1. it follows that the left-neighboring matrix of an E-optimal design has zero diagonal, exactly one off-diagonal element equal to 2 in each row and column, and the remaining off-diagonal elements are equal to 1. This implies that an E-optimal design has two classes of ordered pairs of neighboring treatments: t(t -2) pairs which occur exactly once in a design and t ordered pairs which occur exactly twice in a design. The main problem in the construction of E-optimal designs is determining which pairs occur once in a design and which twice.

Since for different t we have different forms of S d * , the methods of construction of E-optimal designs are different.

Observe that the left-neighboring matrix of circular neighbor balanced design (CNBD), S CNBD , is of the form S CNBD = 1 t 1 t -I t . From Corollary 3.1 it can be noticed that for t = 3, 5, 7 the matrix S d * is obtained by adding one-cycle permutation matrix to S CNBD . Thus, in such a case E-optimal designs can be constructed from CNBD by repeating one arbitrary block.

Example. Let t = 3, 5 or 7. The following designs are E-optimal over D t,t,t : for t = 3 : for t = 5 : for t = 7 :