E Daniau

M Bouchez
email: marc.bouchez@mbda.fr

O Herbinet
email: olivier.herbinet@ensic.inpl-nancy.fr.

P M Marquaire

N Gascoin

P Gillard

Fuel reforming for scramjet thermal management and combustion optimization

Keywords: Combustion Chamber DMR, Dual Mode Ramjet PAH, Polycyclic Aromatics Hydrocarbons λ, thermal conductivity

It is common knowledge that one of the main issues of hypersonic flight is the thermal management of the overall vehicle and more specifically the cooling of the engine, since even composite materials can't withstand the large heat load found in a Scramjet combustion chamber. Another critical point is that mixing and combustion should be sufficiently fast in order to avoid long combustion chamber caused by supersonic internal flow and short residence time. Cryogenic fuels are a logical choice but their lack of storability and low density make them second choice compared to liquid hydrocarbons for small vehicle application. Researches are currently conducted in order to determine if the cooling can be achieved by the endothermic thermal decomposition of the fuel itself circulating trough the engine. The other benefit of this decomposition is the expected shift in the fuel mole fraction, from heavy hydrocarbons (with long induction delays), to light species (mainly H 2 , CH 4 and C 2 H 4 ). In order to quantify the heat transfer in the cooled structures and the composition of the cracked fuel entering the combustor, an accurate predictive model of the thermal decomposition of the fuel is required. For this purpose, an experimental and modeling study of the thermal decomposition of generic molecules (typically long-chain alkanes and polycyclic non-aromatic complex molecules) that could be good surrogates of real fuels, has been started at the DCPR laboratory. The thermal decomposition of those generic molecules is studied in an isothermal jet-stirred reactor. Products of the thermal decomposition are heavy species like aromatics and poly-aromatics and lighter species like hydrogen, methane, ethylene and other middle-weight alkenes. Those results allowed the validation of the kinetic models that have already been developed. During the same time, MBDA-F launched collaborative project named COMPARER, focusing on system analysis to identify one or two characteristic parameters (able to be measured) needed to understand and control the complex phenomena involved in the presented cooling technology and to evaluate some associated sensors. Laminar flame speed for different fuel composition were computed in order to be able to tune the mixture for the lowest possible ignition delay. If at first glance a high mole fraction of hydrogen seems to be interesting, this generally also lead to a high mole fraction of methane (a molecule well known for its long ignition delay and an un-welcomed product) and high mole fractions of middle-weight species (such as ethylene, acetylene, butadiene, cyclopentadiene...) seems to be more promising from a combustion point of view.

Introduction

Hypersonic flight is expected to be achieved with dual-mode-Ramjet (Ramjet under Mach 6 and Scramjet beyond) because of its high specific impulse and its capability to be reusable (especially interesting for space transportation) [START_REF] Falempin | French hypersonic propulsion program PREPHA -results, lessons & perspectives[END_REF] , but one of the main issues at these flight conditions is the thermal management of the engine and the vehicle. Different cooling strategies have been evaluated by MBDA-France (calculations, material tests). Metallic panels have been tested as composite ones (C/SiC) [START_REF] Bouchez | PTAH-SOCAR fuel cooled composite material structure for dual mode ramjet and liquid rocket engine -2005 status[END_REF] , which seem to be promising. But even composite materials could not withstand such large heat load (for example, total temperature of external air reaches 4950 K at Mach 12). Consequently an active cooling system has to be used but not a dedicated one because it would increase the vehicle weight. Furthermore, another issue occurs under theses flight conditions. The time allocated to mix the injected fuel with inlet air, to ignite the combustion and to complete it before the chamber outlet is about 1 ms. These two points lead to the so-called "regenerative cooling" solution : using the fuel to cool down the engine's wall and then burning it in the CC. The fuel is injected in a composite channel (which surrounds the engine) near the outlet of the CC, it flows to the injection on the opposite way of the burned gases. A heavy hydrocarbon fuel is chosen here because of its high density compared to cryogenic fuels (800 kg.m -3 instead of 70/80 kg.m -3 for cryogenic hydrogen, with a specific impulse of liquid hydrocarbon halved) [START_REF] Giraudot | Dual-fueled advanced high-speed ramjets : students paper[END_REF] . When heated and pyrolysed, it produces lighter hydrocarbons species that are both more energetic and easier to ignite. This point allows responding to rapid phenomenon in the CC. But this cooling system requires knowing firstly how the fuel is decomposed and ensures the cooling and secondly how it will burn in the CC (to manage the thrust). It has to be noticed that due to the expected high pressure in the cooling loop (>3 MPa) the fluid becomes supercritical in the channel, which leads to some modelling difficulties (fluid properties and flow rate measurement) for the cooling study. The injected flow rate is expected to be slightly different from the one pumped out of the tank because a carbon deposit (coke) could appears at high temperature (above 1000 K) and because transpiration cooling is planned to be used through the mastered porosity of the C/SiC wall. This phenomenon will also change the carbon/hydrogen ratio of the fuel in the channel. This point needs to be studied because it influences the combustion and a priori the thrust.

American Institute of Aeronautics and Astronautics

II.

Overview of MBDA-France activities

In addition to the long established cooperation with ONERA 4 5 on this topic, MBDA launched two small-scale programs in collaboration with French laboratories and universities.

A. The COMPARER project COMPARER is the French acronym for " COntrol and Measure of PArameters in a REacting stReam". The aim of this project is to identify one or two characteristic parameters (able to be measured) needed to understand and control the complex phenomena involved in the presented cooling technology and to evaluate some associated sensors.

The different actors of this program are MBDA-F, the "Laboratoire Energétique, Explosions et Structures" (LEES, explosions dynamics and reactive systems laboratory), the "Laboratoire de Vision et de Robotique" (LVR, vision and robotics laboratory) and the "Pôle Capteur" (excellence centre for sensors), all located near Bourges, in France.

The used strategy is to first determine the different parts of the engine involved in the cooling. Then, steady simulations have been conducted with a specific code [START_REF] Daniau | Scramjet Active Cooling Analysis Using n-Dodecane as a Generic Endothermic Fuel[END_REF] developed thanks to the existing in-house MBDA program (called NANCY). We use a pyrolysis mechanism [START_REF] Dahm | Experimental and modeling investigation of the thermal decomposition of n-dodecane[END_REF] for the n-dodecane written in CHEMKIN II format. The code allows taking into account the endothermic effects on the wall cooling [START_REF] Gascoin | Mesures et contrôle de paramètres pertinents pour la régulation d'un statoréacteur mixte refroidi[END_REF] .

In order to identify the possible and interesting control measurements, some generic engines were studied in steady conditions, and the different parts of the coupling cooling/bruning loop was analyzed with pluridisciplinary, simplified but unsteady approach.

This sensitivity analysis allowed us to pinpoint critical parameters and to develop strategies to implement the measurement of those parameters in a real engine.

The next step of this program is to build a research test bench (Figure 1) that will enable the development and calibration of specific sensors that could be used for the on-board regulation of a DMR. The fuel is flowing through a metallic or C/SiC reactor installed in a high temperature oven (1800 K) and then passes through a measurement prototype block, then is cooled and burnt in an academic flame.

COMPARER measurement bloc

This working bench will be used to validate the unsteady model of the engine and to evaluate different real time measurement techniques on the decomposed fuel : mass flow, decomposition level and capacity to burn.

The identified techniques can be tested on laboratory level on the fuel, at the exit of the reactor (located as "COMPARER measurement bloc" in Figure 1), with steady or unsteady conditions (temperature, fuel mass flow rate, etc). At least for steady operating points, some characterization is planned both on the decomposition of the fuel and on the academic burner, in order to analyze the COMPARER real time measurements.

Due to security and sizing consideration, the fuel mass flow rate in this reactor will be very low (0.1 g.s -1 ) This way, the burner power is limited to 5 kW and a real time spectroscopy could be used as the main diagnostic device.

The main diagnostic method on this test bench will be the Gas Phase Chromatography coupled to a Mass Spectrometer. Those two methods will be used to cross-checked the results obtained with the COMPARER sensor.

American Institute of Aeronautics and Astronautics

B. MBDA-F -DCPR collaboration

The DCPR is a worldwide-known laboratory heavily involved in the study of fuel oxidation and thermal decomposition and located in Nancy, France. In collaboration with MBDA-F, some kinetic models were developed and validated against experimental data obtained in a perfectly stirred quartz reactor (Figure 2).

The maximum possible temperature in a quartz reactor is near 1300 K. The fuel models include a full scheme primary mechanism, a lumped secondary mechanism and a Polycyclic Aromatic Hydrocarbons (PAH) and soot formation mechanism. A typical mechanism involve more than 1500 elementary reactions and 300 compounds.

Figure 2. A perfectly stirred quartz reactor.

III. Fuel kinetics models

A. Cooled engine modeling

To compute the active cooling of such an engine, the semi-empirical NANCY code is extensively used by MBDA FRANCE. This code deals with stationary heat exchanges between two fluids flowing upstream. Implemented fluids are air, water and hydrogen, but it is possible to input other fluids (such as kerosene or endothermic hydrocarbons [START_REF] Ser | Carburants liquides endothermiques : problématique du refroidissement et projet d'études expérimentales à l'ONERA[END_REF] ) by defining their thermal characteristics versus their temperature and pressure. The code calculates step by step walls and fuel temperatures of each section and thus the cooling of the panels. Moreover it estimates thermal and mechanical stresses in the coolant channel material and computes the pressure drop.

Further investigation is currently taking benefit of more complex modelling such as 3D thermalhydraulic computations 10 or hydrocarbon fuel kinetic modelling, under current development for industrial use by MBDA FRANCE with some research institutes.

Another interest is the capability to check the unsteady behaviour of such a cooled structure.

For the PTAH technology, several system "paper" studies have been performed, mainly on dual-mode ramjets (with hydrogen or hydrocarbon as fuel). Preliminary investigation has also been performed on Liquid-Rocket Engines. B. Fuel selection and modeling Two molecules were selected for the development of a kinetic fuel model. One is a long-chain alkane (type 1) and the other is a non-aromatic, poly-cyclic compound (type 2). Those two very different molecules were selected in order to obtain a better understanding of the thermal behavior of linear and cyclic hydrocarbons.

The primary model is sufficient to describe the decomposition of the initial molecule, but the secondary model is needed to gain the knowledge of the products distribution.

The Figure 3 shows some validation cases of the primary model of the type 1 molecule. This model is relatively simple and contains 500 reactions and 100 compounds.

American Institute of Aeronautics and Astronautics

The Figure 4 shows the conversion profile of the type 2 molecule (cyclic) as a function of the temperature and for a residence time inside the reactor of 1 second.

The dots represent the experimental points and the line is the result of the kinetic primary model. This model is more complex than the previous one, with 1700 reactions and 350 compounds.

Both models shows the utmost importance of the temperature in the decomposition of an endothermic fuel.

In the case of a long-chain alkane, the conversion at 950 K for a very short residence time of 0.05 s is very low (less than 5 %). If we double the residence time (0.1 s), the conversion increases to a modest 5 % and if we double again the residence time (0.2 s), the conversion reaches only 12 %. In contrast, a modest increase of the temperature by 50 K (from 950 K to 1000 K, a small 5 % increase) leads to 17 % conversion (instead of less than 5 %). A further increase by 50 K (at 1050 K) increases the conversion rate between 40 % and 55 %.

C. Fluid composition in a cooling channel

Different test cases have been calculated for various heat wall fluxes, channel length and configuration, fuel flow rate, pressure condition.

The Figure 5 shows the fuel remaining mole fraction as a function of the fuel temperature.

We could conclude that for a limited range of residence time (from tenth of seconds to seconds), the pyrolysis of long-chain alkane is principally temperature dependent except in a small area where the conversion rate is maximal (slightly above 1000 K).

Previous computations taking into account the fuel mass flow rate and the expected heat fluxes show that the fuel temperature in our DMR engine is expected to reach 1300 K and maybe up to 1500 K in some cases, such high temperatures leading to a clear 100 % conversion degree. With the decomposition of the initial molecule, a lot of secondary products could be found in the fuel (Figure 6). For the chained alkane, the main products are ethylene (C 2 H 4 ), ethane (C 2 H 6 ) methane (CH 4 ) and hydrogen at high temperature (H 2 )).

The case of hydrogen is very interesting because it's production rate in the primary mechanism is very low and that any large mole fraction of H 2 imply also a large number of lowly hydrogenated components such as benzene (C 6 H 6 ), toluene (C 7 H 8 ), PAH and even coke.

American Institute of Aeronautics and Astronautics

From Figure 6, it appears that the pyrolysis of long chain alkanes occurs in two step. The first step involve the decomposition of the initial molecule and the formation of a wide range of middle-weight compounds, mainly alkenes and alkanes.

After this step, and if the temperature increases again, the middle-weight compounds start their pyrolysis with the formation of hydrogen, methane and coke For the type 2 (polycyclic) molecule, the products distribution exhibits some notable differences (Figure 7).

In this case, hydrogen is a major product, followed by ethylene and methane (even at "low" temperature). This high hydrogen production is directly linked to the production of benzene and other aromatics like toluene, styrene, indene and naphtalene.

Due to the ring shape of this molecule, a relatively large amount of cyclopentadiene is also produced.

It should be noticed that benzene is produced directly in the primary mechanism, another difference between type 1 and type 2 molecules.

The amount of PAH produced by the cyclic molecule is higher than that produced by the chained alkane, but the PAH formation is not linked to the primary mechanism and is not a direct function of the fuel conversion rate.

Figure 8 shows the remaining mass fraction of the cyclic compound as a function of the residence time and the temperature. The PAH mass fraction is also represented.

For a fuel temperature of 933 K, the conversion of initial fuel is less than 30 %, even after 5 seconds. In this case, the PAH mass fraction is 2 %. The same degree of fuel conversion could be achieved if the fuel temperature is raised to 973 K with a residence time of less than 1 second, but the PAH mass fraction is only 1 %. So if we want to maximise the fuel conversion and minimise the PAH formation, one should select a high fuel temperature with a low residence time (and a short cooling loop with high fuel speed). American Institute of Aeronautics and Astronautics

IV.

Fuel combustion

A. Predicting the fuel burning capability Giving the fact that temperature is the first parameter governing the decomposition of hydrocarbons fuel, this parameter could be sufficient to know the pyrolysed mixture composition before injection in the combustion chamber. But to control the scramjet in terms of thrust, it is needed to predict how the injected fuel can burn. A criterion has to be found to represent the capacity for the mixture to burn. The inflammation delay could be chosen although other criteria like the fundamental flame velocity or the activation energy are of great interest. It can be chosen that the ignition delay should not be greater than 0.1 ms. It corresponds roughly to a tenth of the combustion chamber length as it is confirmed in the literature [START_REF] Buriko | Influence of Active Radical Concentration on Self-Ignition Delay of a Propane/Air Mixture[END_REF] . As we do not have a full combustion mechanism involving all and every species for our main fuel and all its pyrolysis products, we are not able to calculate such delay with PREMIX [START_REF] Kee | A Fortran Program for Modeling Steady Laminar One-Dimensional Premixed Flames[END_REF] and SENKIN [START_REF] Kee | CHEMKIN-II: A Fortran Chemical Kinetics Package for the Anlysis of Gas-Phase Chemical Kinetics[END_REF] but we could at least calculate the delay of a large number of alkanes and alkenes.

Because production of species is not a linear phenomenon, it is possible to expect an "ideal" composition of mixtures for a burning application.

As an example, Figure 6 shows the different species mole fraction in a cooling channel. As some species are more able to burn than others, it is interesting to determine an "ideal" composition which could respond to the chosen criterion of 0.1 ms. Davidenko et al. [START_REF] Davidenko | Kinetic mechanism validation and numerical simulation of hypersonic combustion of methane-hydrogen fuel[END_REF] estimates that a 20 % in mass of hydrogen is needed in a binary mixture of methane and hydrogen for the combustion in scramjet. It should be reminded that hydrogen -air mixtures typically exhibit the smallest ignition delay of all gaseous hydrocarbons and that methane -air mixtures, on the other hand, exhibit the longest ignition delay.

Due to the lack of data, we need to find a way to determine whether a mixture agree with this criterion. We adopt these of , considering that all "good" species in terms of combustion are as good as hydrogen. Then, this choice allows us to compute inflammation delays with PREMIX or SENKIN for any compound found in our fuel mixture, as for example hydrogen, methane, ethane and so on.

The Figure 9 shows some laminar flame speed for different fuel composition. On the left side of the diagram we found the results obtained for two reference DMR and the associated experimentations (with two different specific conditions), on the right the comparison between the computed laminar flame speed of pure compounds and some references, and between the two a generic mixture of methane and hydrogen. It is too early to surely associate a criterion to control an actual regeneratively-cooled hydrocarbon-fueled DMR, but this analyses allows to identify some relationship between the composition or the characteristics (physical, optical, chemical, electrical, ...) of the decomposed fuel and its general capability to burn easily.

We decided to use combustion mechanism found in the literature (like those of San Diego, GRI-Mech) and specially made for certain species. Thanks to these ignition delays for different species or mixtures (function of initial fluid temperature), we have classified compounds in "good" or "bad" in terms of ignition. Products with higher inflammation delay than methane and hydrogen (20% in mass) are considered as bad species. This classification permits to determine whether a mixture contains enough good species to respond to the 0.1 ms criterion or not.

If we look at the Figure 6 again, we see that near the 1 m abscissa the fuel is mainly composed of ethylene, ethane, methane, propene and hydrogen.

American Institute of Aeronautics and Astronautics

Ethylene have a lower laminar flame speed (0.73 m.s -1 ) than hydrogen (2.65 m.s -1 ) or acetylene (1.46 m.s -1 ) but is still one f the easiest burning hydrocarbons. Its molar weight of 28 g is also higher than that of methane (16 g) or hydrogen (2 g), so a high mole fraction of ethylene leads also to a high weight fraction.

On the same figure, we could see that when the temperature increase the fuel products distribution narrowed and if the temperature increase too much, the main products are only methane and hydrogen. In that case, it is also clear that a lot of coke should be produced and that most of the fuel mass flow stays in the cooling loop (forming gums and tar).

So it is clear than even if there is more hydrogen in the mixture at the abscissa of 2.2 m than at the abscissa of 1 m, the mixture will be less efficient in term of burning capability, which is another key parameter (with the cooling capacity, of course) when studying regenerative cooling.

B. Measuring fuel burning capability

In order to be able to control the engine and the fuel ability to burn quickly (composition), some fuel diagnostic are needed onboard the engine. After a review of the possible measurement methods, passive IR spectroscopy was selected as one of the promising and usable ways. Some computations of the fuel IR spectrum, as a function of the fuel temperature and composition, have been done with the HITRAN software [START_REF] Rothman | The Hitran Molecular Spectroscopic Database And Hawks (Hitran Atmospheric Workstation): 1996 Edition[END_REF] . This work enable us to identify some interesting wavelengths. For example, one is located at 4170 cm -1 and its intensity is directly related to the methane mole fraction inside the mixture (Figure 10).

The COMPARER test bench will prove to be a useful tool for testing the application of narrow-band, real-time IR spectroscopy to fuel analysis for engine control.

V. Conclusion

In order to gain some additional knowledge on fuel thermal decomposition for regenerative cooling DMR application, MBDA-France funded two small-scale studies.

With the collaboration of the French DCPR laboratory, a complete model for the pyrolysis of long-chain alkane was developed and validated in a broad range of pressure and temperature. No supercritical phase fluid effect was explicitly taken into account and some discrepancy in pyrolysis products composition was found between computations and literature experimental results, but the general agreement between numerical and experimental results is quite good. This kinetic model was then integrated in a in-house semi-empirical thermal model (NANCY, used by MBDA-France) and a specific package was developed in order to calculate the cooling fluid physical properties as a function of the fluid chemical composition.

This integrated software should be very useful for the comprehension of the mechanism and the design of a regenerative cooled scramjet combustion chamber which are currently under investigation by MBDA.

The success of this approach lead to the study of a second kind of molecule that could be found in a real aviation fuel. The primary kinetic model was developed and experimentally validated. A secondary model, including PAH and coke formation is under development.

During the same time, the COMPARER project focused on the characterization of hot fuel associated to the unsteady behavior of a cooling loop. Some promising measurement are under further analysis and a small but complete specific test bench is currently under development for the testing of those methods. American Institute of Aeronautics and Astronautics

  separation volume -gas analysis w/ M.S and G.P.C. -liquid analysis w/ G.P.C. Temperature gage G.P.C. G.P.C.

Figure 1 .

 1 Figure 1. COMPARER test bench.

Figure 3 .

 3 Conversion of type 1 molecule as a function of temperature and residence time.

Figure 4 .

 4 Figure 4. Conversion of type 2 molecule as a function of the temperature for a 1 s residence time.

Figure 5 .

 5 Figure 5. Main effect of the temperature on the chained alkane pyrolysis (different elementary reactors with various residence times from 0,1 s to 1 s)

1 Figure 6 .Figure 7 .

 167 Figure 6. Repartition of components concentrations in a cooling channel of a generic DMR (long chain alkane).

Figure 8 .

 8 Figure 8. Type 2 molecule remaining mass fraction and PAH formation as a function of fuel temperature and residence time.

propaneFigure 9 .

 9 Figure 9. Laminar flame speed for different pure hydrocarbons and mixtures.

Figure 10 .

 10 Figure 10. intensity of the 4170 cm -1 wavelength as a function of the temperature and the CH 4 mole fraction.

Acknowledgments

The authors wish to thank particularly F. Battin-Leclerc and R. Fournet of the DCPR laboratory and X. Rocourt and S. Bernard of the LEES laboratory.

The COMPARER project is funded by the Conseil Général du Cher (18), the Conseil Régional du Centre, the European Union (FRED, FEDER, FSE) and MBDA-France.


Maxime Clusel

Jean-Yves Fortin
email: fortin@lpt1.u-strasbg.fr

Vladimir N Plechko
email: plechko@thsun1.jinr.ru

Alternative description of the 2D Blume-Capel model using Grassmann algebra

Keywords: numbers: 02.30.Ik, 05.50.+q, 05.70.Fh

We use Grassmann algebra to study the phase transition of the 2D Blume-Capel from a fermionic point of view. This model presents a phase diagram, with a second order critical line which becomes first order through a tricritical point, and was used to model the phase transition in liquid mixtures of He 3 -He 4 . In particular, we are able to map the spin-1 system onto an effective fermionic action from which we obtain the exact mass of the theory. This effective action is actually an extension of the free fermion Ising action with an additional quartic interaction term. The effect of this term is to render the excitation spectrum of the fermions unstable at the tricritical point.

Introduction

The Blume Capel (BC) model is a classical spin 1 model, originally introduced to qualitatively explain the phase transition in a mixture of He 3 -He 4 [1,2,3]. Below a concentration of 67% in He 3 , the mixture undergoes a so-called λ transition: the two components separate through a first order phase transition and only He 4 is superfluid. On a 2D lattice representing an Helium film, He atoms are modelled by a spin-like variable, according to the following rule: an He 3 atom is associated to the value 0, whereas a He 4 is represented by a classical Ising spin taking the value ±1. Within this framework, all the lattice sites are occupied either by an He 3 or He 4 atom. The Blume-Capel model describes the behaviour of this ensemble of spins {σ}. In addition to the usual nearest-neighbour interaction, it includes the term ∆ 0 m,n σ 2 mn = ∆ 0 N s to take into account a possible change in vacancies number. ∆ 0 can be thought as a chemical potential or crystal field. A simple analysis of the BC hamiltonian shows that this model presents a rather complex phase diagram in the plane (T, ∆ 0 ), where T is the temperature in the canonical ensemble [START_REF] Cardy | Scaling and Renormalization in Statistical Physics[END_REF]. In the limit ∆ 0 → -∞, the values σ mn = 0 are excluded and the standard 2D Ising model is recovered, with its wellknown second-order critical point at (T, ∆ 0 ) = (T c = 1/arctanh( √ 2 -1) ≃ 2.27, -∞). At zero temperature, a simple energy argument shows that the ground state is the ordered state |σ| = 1 if ∆ 0 < 2J, and |σ| = 0 else. There is therefore a first order phase transition at (T, ∆ 0 ) = (0, 2J), suggesting a change in the order of the transition at some tricritical point of the critical line. Mean field theory confirms this behaviour and provides a transition line in the plane (T, ∆ 0 ) which is of second order for large negative ∆ 0 . Beyond the tricritical point, the transition becomes first order. Precise numerical simulations have been performed to study precisely the phase diagram and to locate the tricritical point [START_REF] Silva | [END_REF]. From a theoretical backgrounds, several approximations have been used such as mean field theory [6] and high temperature expansion [7]. Using correlation identities and Griffith's and Newman's inequalities, rigourous upper bounds for the critical temperature have also been obtained by Braga et al [8]. The aim of this article is to present a different analytical method, originally proposed by one of us for the 2D classical Ising model in the case of free fermions [9] and since then used to treat various problems around the 2D Ising model, such as finite size effects and boundary conditions [10], quenched disorder [11], boundary magnetic field [12,13]. This method is based on the introduction of Grassmann variables directly in the partition function, in contrast with usual uses of Grassmann variables in this context [14]. It turns out to be particularly efficient to deal with the nearest-neighbour interactions in the 2D plane. As the additional crystal field term in the BC hamiltonian is local, we hoped that the method will be applicable as well in this context. We will see in the following that though it is not possible to compute exactly the partition function and thermodynamics quantities, our approach allows to derive in a controlled way the underlying fermionic lattice field theory. In the continuum limit the effective theory can be analysed in the low energy sector, leading to the exact equation of the critical line and to the effective interaction between fermions that is responsible for the existence of a tricritical point. An approximate method such as the Hartree-Fock method can be used to locate this point. It is interesting to note that a phase diagram with first order transition and tricritical point can be described not only with a bosonic Ginzburg-Landau theory [START_REF] Zinn-Justin | Quantum field theory and critical phenomena[END_REF], where the order parameter is a simple scalar, but also with fermionic variables. The article is organized as follow: after presenting the BC Hamiltonian and the related partition function, we apply the fermionization procedure, leading to the exact fermionic action on the lattice. Then, from this result, we derive the effective action in the continuum and extract the exact mass. We then give the physical interpretation for the existence of a tricritical point on the phase diagram, by studying the fermionic stability of the critical line, and compare our results with recent numerical Monte Carlo simulations.

The 2D Blume-Capel model

Hamiltonian and partition function

The 2D Blume-Capel (BC) model is defined, on a square lattice of linear size L, via the following hamiltonian:

-βH = L m=1 L n=1 K 1 σ mn σ m+1n + K 2 σ mn σ mn+1 + ∆ L m=1 L n=1 σ 2 mn . (1)
In the previous expression K 1,2 = J 1,2 /T represent the nearest-neighbours interactions along the bonds in the two different directions. The crystal field parameter ∆ = -β∆ 0 plays the role of a chemical potential, responsible for the level splitting between σ = 0 and σ = ±1. The partition function Z of the BC model is obtained by summing up all possible configurations of the ensemble {σ mn = 0, ±1}, Z = σ=0,±1 e -βH = Tr {σ} e -βH . Using the property σ 2 mn ∈ {0, 1}, it is easy to develop each Boltzmann factor appearing in the trace formula:

exp (K i σσ ′ ) = 1 + λ i σσ ′ + λ ′ i σ 2 σ ′2 , (2)
with

λ i = sinh K i , λ ′ i = cosh K i -1. (3)
We thus obtain a representation of the partition function as a product of polynomials:

Z = Tr {σmn=0,±1} L m=1 L n=1 e ∆σ 2 mn (1 + λ 1 σ mn σ m+1n + λ ′ 1 σ 2 mn σ 2 m+1n ) × (1 + λ 2 σ mn σ mn+1 + λ ′ 2 σ 2 mn σ 2 mn+1 ) . (4)
This expression of the partition as a product of spin polynomials is the starting point of the fermionization procedure using Grassmann variables. As we will see below, there are different possible ways of introducing Grassmann variables directly in order to disconnect the previous interactions between the local spins. This method is useful to sum up the spin variables and obtain the fermionic action underlying the BC model. The partition function is then given by a path integral of this action, as we will see below.

Local spin decomposition

In the following we will use an alternative interpretation that turns out to be convenient for our derivation. Let us consider a 2D lattice of size L on which we put N s Ising spins S mn = ±1 on each site of coordinates (m, n). If N s = L 2 = N all available sites are occupied (2D Ising model). If N s < L 2 , vacancies appear in the lattice: we can then define on each lattice site a new variable, ρ, that represents the occupation number taking two possible values ρ = 0, 1, so that N s = L m,n=1 ρ mn . The effective spin variable on each lattice site (m, n) is then

σ mn = ρ mn × S mn ∈ {0, +1, -1}.
(

) 5
The usual 2D Ising model is a particular case of this model, whom the configurations space is the subspace where ρ mn = 1 for all (m, n) sites.

In using this decomposition to sum up the spin degrees of freedom, we have to weight the various contributions in order not to count twice the contribution σ = 0. This comes from the fact that there are two different decompositions giving the same σ = 0. While the original averaging over the individual spin means σ=0,±1

f (σ) = f (0) + f (+1) + f (-1), (6)
the naive summation over independent states in the factorized form gives ρ=0,1;S=±1

f (ρS) = f (0) + f (0) + f (+1) + f (-1). ( 7
)
This may be corrected by introducing in the definition of the averaging the weight 1 2 at ρ = 0, or factor 2 ρ-1 , which gives ρ=0,1;S=±1

2 ρ-1 f (ρS) = f (0) + f (+1) + f (-1). ( 8
)
The decomposition scheme ρ mn S mn is close to the one of the site dilute Ising model [11,[START_REF] Plechko | Proceedings of the 6th International Conference on Path Integrals from peV to[END_REF]. The case (7) means simply that there is a spin S mn = ±1 at site ρ mn = 0, which is not interacting with its nearest-neighbours. This empty site, by flipping over two states ±1 under temperature fluctuations, will give however a contribution to the entropy, ln 2 by empty site. The case (8) means that the site ρ mn = 0 is really dilute, or empty, with no spin degree of freedom at it, even disconnected. This interpretation in term of dilute site variables is however quite different from the derivation of a quenchedimpurity model where the set ρ mn is fixed, and the average over these variables is perform after the logarithm of the partition function is taken. The present approach is closer to the annealed disorder case.

Another decomposition that we will not take into consideration in the following computations, but which appear to be convenient as well, relies on the Z 2 gauge symmetry of the model. Indeed, only even terms of the polynomial for Z, which include variables 1, σ 2 mn , but not σ mn , will actually contribute to the sum, since σ mn → 0 under the averaging. Therefore, it is easy to see that Z will not be changed if we multiply each spin variable under the sum by a fixed factor taking values either 1 or -1, which we denote as S mn = ±1, so that the gauge symmetry σ mn → σ mn S mn makes invariant Z[σ mn ] → Z[σ mn S mn ]. Now, though the above invariance under σ mn → σ mn S mn holds already for a fixed set of S mn = ±1, we can also average Z[σ mn S mn ] over S mn = ±1 at each site, with a factor 1 2 to compensate the doubling of the terms, leading to

Z = σmn=0,±1 Z[σ mn ] = σmn=0,±1 Z[σ mn S mn ] = 1 2 L 2 Smn=±1 σmn=0,±1 Z[σ mn S mn ]. (9)
This symmetry outlined in the above transformations can also be used in what follows instead of the occupation number leading, of course, to the same result. As we will see later, it is however essential to use one of these two transformations to simplify the algebra .

Fermionization and lattice fermionic field theory

The expression of the partition function as a product of polynomials, as given by equation ( 4), is the starting point of the fermionization procedure. This procedure has first been introduced by one of us on the 2D Ising model [9]: it relies on interpreting each spin polynomial in (4) as the result of an integral over a set of Grassmann variables. Thought the fermionization is slightly more difficult in the BC case, the idea remains the same. Before going into details, we remind in the following subsection basic features about Grassmann variables and integration/derivation definitions.

Grassmann variables

We remind here the essential features about Grassmann variables that are needed in the rest of the paper. More details can be find in [START_REF] Nakahara | Geometry, topology and physics[END_REF][START_REF] Zinn-Justin | Path Integrals in Quantum Mechanics[END_REF]. A Grassmann algebra A of size N is a set of N anti-commuting objects {a i } i=1,N satisfying:

∀ 1 ≤ i, j ≤ N, a i a j = -a j a i ,
which implies a 2 i = 0. Functions defined on such an algebra are particularly simple, they are polynomials with a finite degree. It is possible to define the notion of integration [START_REF] Nakahara | Geometry, topology and physics[END_REF][START_REF] Berezin | The methods of second quantification[END_REF] with the following rules: da a = 1, da 1 = 0, (10) and for any function f (a),

da f (a) = ∂f (a) ∂a . (11)
With these definitions, Gaussian integrals are expressed by

N i=1 da * i da i exp N i,j=1 a i M ij a * j = det M. (12)
In the following we will also introduce trace operator acting on functions defined on a Grassmann set of variables. In order to easily identify them, such operators will be noted as Tr.

Fermionization trick

In the same spirit than for the 2D Ising model [9], we introduce two pairs of Grassmann variables, (a mn , āmn ), and (b mn , bmn ), to factorize the polynomials appearing in (4). Namely we use the relations

1 + λ 1 σ mn σ m+1n + λ ′ 1 σ 2 mn σ 2 m+1n = dā mn da mn e (1+λ ′ 1 σ 2 mn σ 2 m+1n ) a mn āmn ×(1 + a mn σ mn ) (1 + λ 1 āmn σ m+1n ), (13)
and

1 + λ 2 σ mn σ mn+1 + λ ′ 2 σ 2 mn σ 2 mn+1 = d bmn db mn e (1+λ ′ 2 σ 2 mn σ 2 mn+1 ) b mn bmn ×(1 + b mn σ mn ) (1 + λ 2 bmn σ mn+1 ). ( 14
)
For sake of simplicity we can introduce the following link factors

A mn = 1 + a mn σ mn , Ām+1n = 1 + λ 1 āmn σ m+1n , (15)
B mn = 1 + b mn σ mn , Bmn+1 = 1 + λ 2 bmn σ mn+1 , (16)
and the Grassmann local trace operators which associate to any function f on the Grassmann algebra

Tr (amn) [f ] = dā mn da mn e (1+λ ′ 1 σ 2 mn σ 2
m+1n ) a mn āmn f (a mn , āmn ), ( 17)

Tr (bmn) [f ] = d bmn db mn e (1+λ ′ 2 σ 2 mn σ 2 mn+1 ) b mn bmn f (b mn , bmn ). ( 18
)
The Boltzmann weights then reads

1 + λ 1 σ mn σ m+1n + λ ′ 1 σ 2 mn σ 2 m+1n = Tr (amn) [A mn Ām+1n ], (19) 1
+ λ 2 σ mn σ mn+1 + λ ′ 2 σ 2 mn σ 2 mn+1 = Tr (bmn) [B mn Bmn+1 ]. (20)
Introducing those Grassmann factors in the expression (4) of the partition function, we then obtain a mixed representation containing both spins and Grassmann variables. One should note that as A mn , Āmn , B mn , Bmn are neither commuting or anti-commuting objects, the order in which it appears in the product is therefore important.

We then define the following trace operator :

Tr a,b [f ] = L m=1 L n=1
dā mn da mn d bmn db mn e ∆σ 2 mn f (a mn , āmn , b mn , bmn )

× exp

L m=1 L n=1 (1 + λ ′ 1 σ 2 mn σ 2 m+1n )a mn āmn + (1 + λ ′ 2 σ 2 mn σ 2 mn+1 )b mn bmn .
Note that all even power of the spin variables are put into the general measure of integrals, including the chemical potential term which is only local. The partition function is given by

Z = Tr {σ} Tr a,b   -→ L n=1   -→ L m=1 (A mn Ām+1n )(B mn Bmn+1 )     . ( 22
)
At this stage the partition function appears as a double trace, over the spin degrees of freedom, with Tr {σ} , and over the Grassmann variables, with Tr a,b .

Spin summation

Up to now we only add extra variables to obtain the mixed expression [START_REF] Mattuck | A guide to Feynman Diagrams in the Many-Body Problem[END_REF]. Further algebraic manipulations are necessary to simplify this expression, in particular to move the link factors through the product, in order to put together those having the same indices (or same spin variable), i.e. A mn , Āmn , B mn and Bmn . Then, using the decomposition σ mn = ρ mn S mn , a partial sum over the variables S mn is possible, since the general trace Tr a,b depends only on the squares σ 2 mn = ρ 2 mn S 2 mn = ρ 2 mn , therefore only on the ρ mn s. To do so, we use the mirror-ordering introduced in the 2D Ising model [9], to move together the different link factors containing the same spin. Up to boundary terms irrelevant in the thermodynamic limit, and using the notation of the ordered products [9], this leads to

Z = Tr σ Tr a,b L m=1 L n=1 (A mn Ām+1n )(B mn Bmn+1 ), , = Tr σ Tr a,b -→ L n=1 -→ L m=1 Bmn A mn Ām+1n • ←- L m=1 B mn , = Tr {σ} Tr a,b   -→ L n=1     -→ L m=1 Āmn Bmn A mn   •   ←- L m=1 B mn       . (23)
In principle we should pay attention to boundary terms, which can actually be treated rigourously [10,12,13]. However, as we are interested only in the continuous limit of the BC model, this would lead to a more complicate and unnecessary analysis : we can neglected those boundary terms since they do not contribute to the thermodynamic limit.

In the case of the 2D Ising model, as σ 2 mn = 1, we can explicitly perform the trace over the spins recursively at the junction of two products, directly from the last relation in (23). The situation is slightly different in the BC case, since σ 2 mn = 0, 1 instead. As we said previously, the trace operator (21) contains terms σ 2 mn = ρ 2 mn = ρ mn that are coupled to the neighbouring sites. Therefore it is not possible to trace over the individual spins σ mn directly, but only over the variables S mn in a first stage. In other words, the ordering procedure on the link variables allows us to eliminate one part of the degrees of freedom. At the junction of the two products in (23), we begin to perform the trace for m = L or, in general, 1 2

Tr Smn Āmn Bmn A mn B mn = 1 2

Smn=±1

Āmn Bmn A mn B mn ,

= 1 + a mn b mn + λ 1 λ 2 ām-1n bmn-1 + (λ 1 ām-1n + λ 2 bmn-1 )(a mn + b mn ) +λ 1 λ 2 ām-1n bmn-1 a mn b mn ρ mn , = exp a mn b mn + λ 1 λ 2 ām-1n bmn-1 + (λ 1 ām-1n + λ 2 bmn-1 )(a mn + b mn ) ρ mn .

As this term is quadratic in Grassmann variables, it commutes with all other terms, and can therefore be move out from the ordered product. We then repeat the calculation for m = L -1 and so on. The partially traced partition function finally reads

Z = 2 N Tr ρ L m,n=1
dā mn da mn d bmn db mn exp ∆ + a mn b mn + λ 1 λ 2 ām-1n bmn-1

+(λ 1 ām-1n + λ 2 bmn-1 )(a mn + b mn ) ρ mn + (1 + λ ′ 1 ρ mn ρ m+1n )a mn āmn +(1 + λ ′ 2 ρ mn ρ mn+1 )b mn bmn . ( 25
)
The trace over the ρ variables can be performed by decoupling them in the exponential term. Several methods are possible. One of them is to introduce another set of Grassmann link variables as we previously did to decouple the factors ρ mn ρ m+1n and ρ mn ρ mn+1 . It is possible however to avoid the introduction of new fields by using instead the following rescaling trick. It is achieved by a change of variables which leaves the integral invariant :

a mn → a mn ρ mn , b mn → b mn ρ mn and (26)
da mn → ρ mn da mn , db mn → ρ mn db mn .

The variable ρ mn disappears in some places inside the exponential and the factors that were previously linked ρ mn ρ m+1n and ρ mn ρ mn+1 are now decoupled. Also, the seemingly singular expressions like ρ mn exp(a mn āmn /ρ mn ) are to be understood as ρ mn exp(a mn āmn /ρ mn ) = ρ mn (1 + a mn āmn /ρ mn ) = ρ mn + a mn āmn . After shifting some indices, we obtain

Z = 2 N Tr ρ L m,n=1
dā mn da mn d bmn db mn (ρ mn + a mn āmn )(ρ mn + b mn bmn ) (28)

× exp ρ mn ∆ + λ 1 λ 2 ām-1n bmn-1 + λ ′ 1 a m-1n ām-1n + λ ′ 2 b mn-1 bmn-1 × exp a mn b mn + (λ 1 ām-1n + λ 2 bmn-1 )(a mn + b mn ) .
Now we can locally perform the sum over ρ mn , using the rule (8), since we have to pay attention not to count twice the contribution σ mn = 0. In that case we have, for the part depending on ρ mn , Tr ρmn=0,1 (ρ mn + a mn āmn )(ρ mn + b mn bmn ) (29)

× exp ρ mn ∆ + λ 1 λ 2 ām-1n bmn-1 + λ ′ 1 a m-1n ām-1n + λ ′ 2 b mn-1 bmn-1 = 1 2 a mn āmn b mn bmn + 2e ∆ e Gmn ,
where G mn is the quadratic form

G mn = a mn āmn + b mn bmn + λ 1 λ 2 ām-1n bmn-1 + λ ′ 1 a m-1n ām-1n + λ ′ 2 b mn-1 bmn-1 . (30)
The last term in equation ( 29) can be easily written as an exponential :

a mn āmn b mn bmn + 2e ∆ e Gmn = 2e ∆ e Gmn 1 + 1 2 a mn āmn b mn bmn e -∆-Gmn ,

= 2e ∆ exp G mn + 1 2 a mn āmn b mn bmn e -∆-Gmn .

Just to simplify the comparison with the 2D Ising model, we can rescale some Grasmann variables using the following transformation:

(1 + λ ′ 1 )ā mn → āmn , (1 + λ ′ 2 ) bmn → bmn . ( 32
)
Introducing

t i = tanh K i = λ i /(1 + λ ′ i )
, we obtain the final result of this section, namely the representation of the 2D BC model as a fermionic integral :

Z = 2 N (cosh K 1 cosh K 2 ) N e N ∆ L m,n=1
dā mn da mn d bmn db mn (33) exp a mn āmn + b mn bmn + a mn b mn + (t 1 ām-1n + t 2 bmn-1 )(a mn + b mn ) +t 1 t 2 ām-1n bmn-1 + g 0 a mn āmn b mn bmn e -γ 1 a m-1n ām-1n -γ 2 b mn-1 bmn-1 -t 1 t 2 ām-1n bmn-1 ,

where we have introduced the following constants

g 0 = e -∆ 2 cosh K 1 cosh K 2 , t i = tanh K i , γ i = 1 - 1 cosh K i = 1 -1 -t 2 i . (34)
We can recognize in (33) the Ising action

S Ising = L m,n=1
a mn āmn + b mn bmn + a mn b mn

+(t 1 ām-1n + t 2 bmn-1 )(a mn + b mn ) + t 1 t 2 ām-1n bmn-1 , (35)
with an additional interaction part which is a polynomial of degree 8 in Grassmann variables (after expanding the exponential):

S int = g 0 L m,n=1
a mn āmn b mn bmn e -γ 1 a m-1n ām-1n -γ 2 b mn-1 bmn-1 -t 1 t 2 ām-1n bmn-1 .

The BC model differs from the Ising model by the interaction term (36), which is not quadratic, therefore it is not solvable in the sense of free fermions. When g 0 = 0, i.e. in the limit ∆ → ∞, the gap between the two degenerate states σ = ±1 and the singlet σ = 0 becomes so large that the model reduces effectively to the 2D Ising model. In the following, for simplification, we will only consider the homogeneous coupling case, where

t 1 = t 2 = t and γ 1 = γ 2 = γ = 1 - √ 1 -t 2
, and we will note S = S Ising + S int the total action.

Partial bosonization

The previous action contains two pairs of Grassmann variables per site. This can not be reduced to one pair (minimal action) unlike the Ising model, where half of the variables are irrelevant in the sense they do not contribute to the critical behaviour. Here the 2 pairs are coupled together by Eq. (36). However, as we will see in the following, it is possible to recover the minimal Ising action [11] using bosonic variables. In the previous interaction action (36), it is indeed tempting to replace the products a mn āmn and b mn bmn , which look similar to occupation number operators, by the commuting variables ρ mn = a mn āmn , η mn = b mn bmn , (37)

ρ mn ρ m ′ n ′ = ρ m ′ n ′ ρ mn , and ρ 2 mn = 0. ( 38
)
This new variables are nilpotent (as Grassmann variables) but commuting: that is why in the following we will abusively call them (hard core) "bosons". This allows us to reduce the degree of the polynomials in Grassmann variables by a factor 2 each time the replacement is performed, even if terms like ām-1n bmn-1 in (36) can not be replaced. We will show below that we can write down and action containing one pair of Grassmann variables and one pair of bosonic ones per site. To do so, we introduce the following Dirac distribution for any function f of a mn āmn or b mn bmn :

f (a mn āmn ) = dρ mn dρ mn f (ρ mn ) exp [ρ mn (ρ mn + a mn āmn )] ,

and f (b mn bmn ) = dη mn dη mn f (η mn ) exp ηmn (η mn + b mn bmn ) .

Applying (39,40) directly in (33), we finally obtain the following expression for the action

S = m,n
a mn āmn + b mn bmn + t 2 ām-1n bmn-1 + a mn b mn + t(ā m-1n + bmn-1 )(a mn + b mn ) +g 0 ρ mn η mn 1 -γ(ρ m-1n + η mn-1 ) + γ 2 ρ m-1n η mn-1 -t 2 ām-1n bmn-1 +ρ mn (ρ mn + a mn āmn ) + ηmn (η mn + b mn bmn ).

We can now integrate over the a mn s and b mn s, and replace formally, for convenience, the variables āmn by a mn and bmn by āmn . We obtain

S = L mn=1 t 2 a m-1n āmn-1 + t(a m-1n + āmn-1 )[ā mn (1 + ηmn ) -a mn (1 + ρmn )] (42)
+ā mn a mn (1 + ηmn )(1 + ρmn ) + ρmn ρ mn + ηmn η mn +g 0 ρ mn η mn 1 -γ(ρ m-1n + η mn-1 ) + γ 2 ρ m-1n η mn-1 -t 2 a m-1n āmn-1 .

The next operation is to notice that we can further integrate over the variables η mn and ηmn , using the integration rules (40), since dη mn dη mn exp ηmn (η mn -tā mn (a m-1n + āmn-1 ) + āmn a mn (1 + ρmn )) f (η mn ) = f [-tā mn (a m-1n + āmn-1 ) + āmn a mn (1 + ρmn )].

We could also have chosen to integrate over the ρ mn s and ρmn s instead. We obtain finally S = āmn a mn + t 2 a m-1n āmn-1 + t(a mn -āmn )(a m-1n + āmn-1 ) (44) +ρ mn ρ mn + ρmn āmn -t(a m-1n + āmn-1 ) a mn

+g 0 ρ mn Q mn 1 -γ(ρ m-1n + Q mn-1 ) + γ 2 ρ m-1n Q mn-1 -t 2 a m-1n āmn-1 , with Q mn = āmn [a mn (1 + ρmn ) -t(a m-1n + āmn-1 )]. (45)
We recognize in the first line of the previous expression (44) the minimal action for the Ising model with only one pair of Grassmann variables [11] S Ising = āmn a mn + t 2 a m-1n āmn-1 + t(a mn -āmn )(a m-1n + āmn-1 ). ( 46)

The rest of the action describes the interaction between the fermions and bosons.

S int = ρmn ρ mn + ρmn āmn -t(a m-1n + āmn-1 )

a mn + g 0 ρ mn Q mn 1 -γ(ρ m-1n + Q mn-1 ) +γ 2 ρ m-1n Q mn-1 -t 2 a m-1n āmn-1 . (47)
It is easy to check that when g 0 = 0 the boson variables can be integrated out.

We would like to end this section commenting the previous results. We finally obtained a lattice field theory containing the same number of "fermions" (a, ā) and "bosons" (ρ, ρ). Physically, this means that it is indeed possible to describe the system with fermionic variables for the states σ = ±1 and bosonic ones for the third state σ = 0. In the limit ∆ 0 → -∞, the system is completely described in terms of fermions. While lowering ∆ 0 , an interaction between fermions is added. Beyond a value ∆ 0t , fermions form bosonic pairs : in the limit ∆ 0 → +∞, all the fermions condense into bosons, leading to an purely bosonic systems. In this interpretation, the tricritical point may be interpreted as the particular point on the critical line where the interaction is such that an additional symmetry between fermions and bosons appears. This would correspond to super-symmetry appearing in the conformal field theory describing tricritical Ising model [START_REF] Ph | Conformal Field Theory[END_REF]. To our knowledge there is no evidence of super-symmetry derived directly form a lattice model : action (44) could be a good way to see how super-symmetry emerges from a lattice model. Of course all we said so far is only speculative : we are currently studying it in more details, to confirm or infirm this hypothesis.

Effective action in the continuum limit

In the 2D Ising model case, the fermionic action on the lattice is quadratic, and the Grassmann integrals equivalent to expression (33) with g 0 = 0 can actually be computed exactly. The situation for the 2D BC model is less simple, as the interaction (36) contains terms of order up to 8th. The Grassmann integral leading to the partition function can no longer be decomposed into elements of size independent of the system size : in this sense the 2D BC model is not integrable. However it is still possible to extract physical information by taking the continuous limit of the action (33) and analysing it using tools from field theory.

Effective 2nd order fermionic field theory

We would like to obtain an effective theory up to order 2 in momentum k, from the previous calculations, in order to analyse the critical behaviour of the model. In the Ising model, the critical behaviour is given, in the continuous limit, by a mass-less theory.

In the following, we will see how to compute the mass of the model. In the infrared limit, the spectrum is given by expanding the action in second order in the momentum k. The "stiffness" of the Ising model, or the coefficient in front of the term k 2 , is always a strictly positive coefficient. Here in the BC model, we will show that it can vanish at the tricritical point, rendering the spectrum unstable and changing the nature of the singularity.

Instead of integrating over the variables η mn and ηmn as in Eq. ( 44), we proceed by integrating first over ηmn and ρmn from Eq. ( 42)

Z ∆ = 2 N cosh(K) 2N e N ∆ m,n
dc mn dc mn dρ mn dη mn [c mn c mn + ρ mn q mn + η mn qmn + ρ mn η mn ]

× exp(S Ising + S int ),
where

S int = g 0 m,n ρ mn η mn (1 -γρ m-1n )(1 -γη mn-1 ) -t 2 c m-1n cmn-1 ,
and qmn = cmn c mn + tc mn (c m-1n + cmn-1 ), q mn = cmn c mn -tc mn (c m-1n + cmn-1 ). ( 48)

It is also useful to notice that q mn qmn = 0. The action S Ising is given by the relation (46). The previous integral involves the product of monomials such as cmn c mn + ρ mn q mn + η mn qmn + ρ mn η mn which can not be written as a single exponential. However, when integrating over ρ mn and η mn , it is easy to realize that these terms roughly impose the following substitution in the action S int :

ρ mn η mn → cmn c mn (49)
ρ mn → qmn (50)

η mn → q mn . (51)
In one sense, this product can be considered as an approximate Dirac function for the variables ρ mn and η mn . However it is not exact: when expanding the exponential of S int it appears terms that couple to each other to give cmn c mn and not qmn q mn = 0 as given by the substitution rules. For example terms such as (g 0 ρ m+1n η m+1n γρ mn ) × (g 0 ρ mn+1 η mn+1 γη mn ) (52) leads to a contribution equal to

g 2 0 γ 2 cmn c mn cm+1n c m+1n cmn+1 c mn+1 . ( 53
)
Therefore there are more terms in the final effective action S eff (c mn , cmn ) than in the one resulting from the substitution rules. However, we would like to obtain an effective action up to order 2 in momentum k, in order to study the stability of the free fermion spectrum which gives basically a log(k 2 ) contribution in the free energy for the Ising model at criticality. The terms that contribute to the second order in momentum are, from the free quadratic part Eq. 46 and in the continuous limit, combinations of c∂ x c or c∂ x c, and c∂ y c or c∂ y c. We indeed expect that the effective action will contain quartic contributions such as cc∂ i c∂ j c, with i, j = x, y at the lowest order. This term is order 4 in Grassmann variables and 2 in derivatives. The expansion of the exponential of such terms will give corrective coefficients to the k 2 behaviour, and thus change the order of the transition if such renormalized coefficient vanishes. We also have to consider not only the direct substitution of the variables with the rules given above, but also the possible corrective terms that can lead to the same order. Also we should drop terms which contain a ratio of number of derivatives to the number of Grassmann variables higher strictly than 1/2. After some algebra, these following terms contribute to the effective action :

S eff = S Ising + g 0 m,n cmn c mn (1 -γ qm-1n )(1 -γq mn-1 ) -t 2 c m-1n cmn-1 +g 2 0 γ 2 m,n cmn c mn cm+1n c m+1n cmn+1 c mn+1 + . . . (54)
In the continuous limit we replace c mn by c = c(x, y) and cmn by c = c(x, y), with the rules

c m-1n = c -∂ x c + ∂ xx c/2! + . . . For example, q mn → q = cc(1 -t) + tc(∂ x c + ∂ y c) and qmn → q = cc(1 -t) -tc(∂ x c + ∂ y c).
The free Ising part gives simply

S Ising = dxdy (t 2 + 2t -1)cc + t(t + 1)c∂ x c -t(t + 1)c∂ y c -tc∂ x c + tc∂ y c . (55)
In the continuous limit, the last term of Eq. (54) give cc∂ x c∂ x c∂ y c∂ y c, which is order 4 in derivatives and 6 in Grassmann variables. The ratio of these numbers is 2/3 which is higher than 1/2, and therefore this term can be dropped, as explained above. The term in factor of g 0 in Eq. (54) contain qm-1n and q mn-1 which need to be expanded up to the order 2 in derivatives, with ∂ xx q = 2(1 -t)∂ x c∂ x c, and ∂ yy q = 2(1 -t)∂ y c∂ y c. The effective action finally can be written in the continuous limit as

S eff = S Ising + dxdy g 0 cc + g 0 cc[t(t + 2γ)∂ y c∂ x c + γ(1 -t)(∂ x c∂ x c + ∂ y c∂ y c)] + • • • (56)
In the following we use this effective action to obtain information on the phase diagram of the BC model.

Spectrum analysis and phase diagram

In this section we analyse the critical properties of the effective action (56) and the low energy spectrum. In particular we develop a physical argument for the existence of a tricritical point on the phase diagram from a fermionic action.

Critical line

From Eqs. (55) and (55), we can easily identify the mass M 0 = -t 2 -2t + 1 + g 0 in factor of the term cc, and which vanishes on the critical line M 0 = 0,

tanh 2 1 T + 2 tanh 1 T -1 = e ∆ 0 T 2 cosh 2 1 T . (57)
On this line, only derivative contributions remain in Eq.( 55) which give the free fermion kinetic terms, and the ones in Eq.( 56) give the interaction between the singlet level and the Ising doublet at the quartic level. The critical line given by ( 57) is plotted in figure 1, and compared with recent Monte Carlo simulations by da Silva et al. [START_REF] Silva | [END_REF]. The agreement between numerical simulations and our results is very good, the mass of the system (57) being exact. Note that our results is also compatible with exact upper bond obtained by Braga et al. [8].

Tricritical point : Hartree-Fock analysis

The main physical feature of the 2D BC model is the existence of a tricritical point, on the critical line. Below this point, the phase transition goes from second order to first order: the tricritical point is characterized by a change in the nature of the singularity. This change should be seen in the BC spectrum from (56). In this section, we analysis the effect of the quartic terms, on the critical line g 0 = t 2 + 2t -1 by considering the effect of the quartic interaction on the stability of the free fermion spectrum. The Ising part can be easily written in the Fourier space, after having defined the following transformations

c(r) = 1 L k c k exp(ik.r) and c(r) = 1 L k ck exp(-ik.r). ( 58
)
Using these transformations, the Ising part is almost diagonalized, in particular

S Ising = k∈S it(t + 1)(k x -k y )(c k c k -c-k c -k ) + 2itk x c k c -k + 2itk y ck c-k , ( 59
)
where S is the set of Fourier modes that correspond to half of the Brillouin zone, i.e. we choose k in S such that -k is not included in S (to avoid repetition of modes in the different sums above), and such that the couples of modes (k, -k) fill up the Brillouin zone exactly once. The quartic term can be written in the Fourier space as

S int = 1 L 2 k 1 +k 2 =k 3 +k 4 V (k 2 , k 4 )c k 1 c k 2 ck 3 ck 4 , (60)
with the potential

V (k 2 , k 4 ) = -αk x 2 k y 4 + β(k x 2 k x 4 + k y 2 k y 4 ) α = g 0 t(t + 2γ) and β = g 0 γ(1 -t). (61)
Up to now we only expressed the action in the Fourier space, without further approximations. In order to see if a tricritical point appears in this approach, we use a mean field like approximation, similar to the Hartree-Fock method [START_REF] Thouless | The quantum mechanics of many-body systems[END_REF][START_REF] Mattuck | A guide to Feynman Diagrams in the Many-Body Problem[END_REF]. To do so, we decompose the fourth order interacting terms into sums of quadratic terms with coefficients to be determined self-consistently. These coefficients are actually two point correlation functions. The interaction can be decoupled in different ways. For example, considering the terms contributing to the Ising action, we may take account of the averages ck c k , c-k c -k , c k c -k and ck c-k . There are also three different ways to decouple the interacting term, since c k 1 can be paired with either c k 2 , ck 3 or ck 4 . For example,

c k 1 c k 2 = c k 1 c k 2 + (c k 1 c k 2 -c k 1 c k 2 ) ≡ c k 1 c k 2 + δ c 1 c 2 , (62)
where δ c 1 c 2 is assume to be a small fluctuation. In this case, from Eq. (59), the average is non zero only for k 1 = -k 2 = k or -k, with k ∈ S. We can pair the other terms by writing the action in the three different possible ways that are compatible with Eq. (59), and by using the fermionic rules

S int = 1 L 2 g k 1 +k 2 =k 3 +k 4 V (k 2 , k 4 ) ( c k 1 c k 2 + δ c 1 c 2 )( ck 3 ck 4 + δ c3 c4 ) (63) -( c k 1 ck 3 + δ c 1 c3 )( c k 2 ck 4 + δ c 2 c4 ) + ( c k 1 ck 4 + δ c 1 c4 )( c k 2 ck 3 + δ c 2 c3
) ,

with g = 3. The next step is to discard terms that are proportional to the squares of fluctuations δ 2 , and keep the others. After some algebra, we obtain the mean-field quadratic operator for the interaction term

S int = 1 L 2 g k,k ′ ∈S 4c k c -k ck ′ c-k ′ V (k, k ′ ) + 4c k c-k c k ′ c -k ′ V (k ′ , k) +c k c k ck ′ c k ′ v(k, k ′ ) + c-k ′ c -k ′ v(k, -k ′ ) +c -k c -k ck ′ c k ′ v(-k, k ′ ) + c-k ′ c -k ′ v(-k, -k ′ ) , (64)
where we have defined the potential

v(k, k ′ ) = -V (k, k) -V (k ′ , k ′ ) + V (k, k ′ ) + V (k ′ , k). (65)
In the previous expressions, there are three different kinds of quantities that contribute to the action, for example sums involving quantities like

k ck c k , k ck c k k i or k ck c k k i k j ,
with i, j = x, y. The first term gives a contribution to the total mass, the second one corresponds to current operators, and the third one can be thought as a dispersion energy tensor. Considering the symmetries of the Ising part, and the fact that the action should be invariant by a dilation factor at criticality, we may only take account of the current operators. We define therefore the following unknown parameters, for fermions and pairs of fermions or "Cooper pairs" (i = x, y)

t i = i 2L 2 k∈S ( ck c k -c-k c -k )k i , u i = i L 2 k∈S c k c -k k i , ūi = i L 2 k∈S ck c-k k i . (66)
From the previous discussion, we can drop the first two terms in the potential v(k, k ′ ) defined in Eq. ( 65), since we already assume that only currents are kept as parameters along the critical line. In this case, it is easy to rewrite, from the property v(k, k

′ ) = -v(-k, k ′ ) = -v(k, -k ′ )
, the effective mean field action (64) as :

S int = 1 g k∈S 4ic k c -k [(αū y -β ūx )k x -β ūy k y ] + 4ic k c-k [-βu x k x + (αu x -βu y )k y ] +2i(c k c k -c-k c -k )[(αt y -2βt x )k x + (αt x -2βt y )k y ]. (67)
We make the further assumption that, by symmetry invariance in the space, there exists a solution satisfying ūy = u x , ūx = u y and t x = -t y :

S int = 1 g k∈S 4ic k c -k [(αu x -βu y )k x -βu x k y ] + 4ic k c-k [-βu x k x + (αu x -βu y )k y ] -2i(c k c k -c-k c -k )(k x -k y )(α + 2β)t x . (68)
The total effective action can finally be written as

S eff = k∈S i t(t + 1) - 2 g (α + 2β)t x (k x -k y )(c k c k -c-k c -k ) (69) +i 4 g g 2 t + (αu x -βu y ) k x -βu x k y c k c -k + i 4 g -βu x k x + g 2 t + (αu x -βu y ) k y ck c-k ,
or in a more compact form as

S eff = k∈S ic(k x -k y )(c k c k -c-k c -k ) + 2i(ak x -bk y )c k c -k + 2i(-bk x + ak y )c k c-k , (70)
with the following coefficients

a = t + 2 αu x -βu y g , b = 2β u x g , c = t(t + 1) -2t x α + 2β g . (71)
The partition function can then be written as a product over the Fourier modes Z = k∈S Z k , with

Z k = k 2 [A + B sin 2θ], (72)
θ being the angle of the vector k, and

A = c 2 -4ab, B = -c 2 + 2(a 2 + b 2 ). (73)
We assume that |A| is larger than |B| on the second order critical line, until the tricritical point is reached, where eventually A 2 = B 2 . Indeed, the expression (72) is valid only if the elements A + B sin 2θ are all positive, which is the case only if A 2 > B 2 . This will be checked using numerical analysis. Beyond this point, the effective action is unstable and has to be modified to incorporate further corrections. In a ϕ 6 Ginzburg-Landau theory describing a first order transition, this will be equivalent to say that at the tricritical point, not only the mass coefficient of the ϕ 2 term vanishes, but also the coefficient of the ϕ 4 term. The parameters t x , u x and u y are determined self-consistently from the definitions Eqs. (66). In the continuous limit, these reduce to

t x = c 4π π 0 dθ 1 -sin 2θ A + B sin 2θ , u x = 1 2π π 0 dθ a sin 2θ -b A + B sin 2θ , u y = 1 2π π 0 dθ a -b sin 2θ A + B sin 2θ . ( 74
)
After computing the trigonometric integrals, we obtain the relations

t x = c 4B -1 + (A + B) sign(A) √ A 2 -B 2 , u x = 1 2B a -(aA + bB) sign(A) √ A 2 -B 2 , u y = 1 2B -b + (bA + aB) sign(A) √ A 2 -B 2 . ( 75
)
To find the tricritical point, we proceed the following way. For each temperature T , we solve the consistency equations for t x , u x and u y with the value of ∆ 0 given by the critical line (57). The self-consistent values are then plug into the coefficients A(T ) and B(T ).

We then plot A(T ) 2 -B(T ) 2 as a function of T , as shown on figure 2. By doing so we find a tricritical point approximatively located at (T t , ∆ 0,t ) ≃ (0.42158, 1.9926). Monte Carlo simulations [START_REF] Silva | [END_REF] gives a T t ≃ 0.61. As expected the mean-field like treatment of the underlying field theory underestimates the fluctuations, rendering the second order critical line more stable at lower temperature, as we approach (T = 0, ∆ 0 = 2). Stronger fluctuations can be simulated by lowering the value of g, which increases (lowers) the value of T t (∆ 0 ) respectively. Instead of g = 3, taking g = 2.5 for example leads to a T t ≃ 0.48, closer to the Monte Carlo result. This can be achieved precisely by incorporating more diagrams in the computation of the effective free energy [START_REF] Mattuck | A guide to Feynman Diagrams in the Many-Body Problem[END_REF]. Also, due to the fact we are in a region near (T = 0, ∆ 0 = 2) where the change in temperature is large compare to the change of ∆ 0 (the slope is vertical at this point as seen in figure 1), it is more difficult to obtain a precise value of T t with a mean-field treatment.

Conclusion

In this paper, we have considered the physics of the BC model as a fermionic field theory. Using Grassmann algebra, we have shown that the model can be transformed as an exact fermionic action on the discrete lattice. This action can be reduced, in the continuum limit, to an effective field theory which is a modified Ising action with a quartic interaction. From there we have extracted the exact mass of the model, and analysed the effect of the quartic term on the stability of the free fermion spectrum.

In particular, the stiffness of the excitation spectrum (the coefficient in front of the k 2 term as we expand the dispersion relation in momentum variable) vanishes at a singular point we have identified as the tricritical point. A Hartree-Fock analysis gave an approximate location of this point on the phase diagram which can be compared to the numerical results of Monte Carlo simulations. The precise location of the point could be achieved by taking into account further diagrams contributing to the effective free energy. However, we have shown the existence of a singular point by studying the stability of the kinetic spectrum where the nature of the transition is changed. The main result of this paper is the possibility to study precisely first-order transition driven systems from a fermionic point of view using Grassmann algebra. That method can be useful for systems where the effective field theory is described by a similar action such as Eq. ( 56) which is one of the simplest form that can be written out from a unique pair of Grassmann variables. Application of this method to other extensions of the BC Hamiltonian, such as the Blume-Emery-Griffiths model [3], is also possible. Finally, a partial bosonization of the system leads to a mixed representation of the model not only in term of fermions but also in term of hard core bosons, as written explicitely in Eq. 44. This result could be useful to see the tricritical point as a special location in the phase diagram where a possible hidden symmetry between fermion and boson variables may appear.

Figure 1 .

 1 Figure 1. Comparison between the critical line (57) (plain line) and numerical results from Monte Carlo simulations by da Silva et al. [5] (dots). The last numerical point (arrow) is the tricritical point identified by da Silva et al..

2 Figure 2 .

 22 Figure 2. Stability of the spectrum: coefficient A(T ) 2 -B(T ) 2 as a function of T . The temperature where A(T ) 2 -B(T ) 2 = 0 gives the location of the tricritical point.


Alois Würger

Viscoelastic properties of confined molecular layers

We study the viscoelastic properties of a film of n layers of spherical molecules confined between two walls. We find that the dynamic response arises from two competing contributions: the effective stiffness of n+1 springs in series and softening due to strain fluctuations. In particular, the latter are of the origin of the oscillatory behavior of stiffness and damping coefficient. Dissipation is strongest at the minima of the stiffness; the inverse behavior may occur for a modulated relaxation time. As a corollary we show that confined molecular layers cannot be described as Maxwell fluids.

I. INTRODUCTION

Molecular liquids confined at a nanometer scale undergo striking structural changes. Their viscoelastic properties differ significantly from those of the bulk and cannot be understood in terms of continuum fluid mechanics [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF]. Surface-force apparatus (SFA) and atomic-force microscopy (AFM) studies on molecular liquids report oscillatory solvation forces and give evidence for the layered structure of the film [START_REF] Horn | [END_REF][3][4][5][6], in agreement with numerical simulations of the molecular dynamics [7]. When reducing its thickness L one observes a sequence of solid-like states, which are separated by rather unstable configurations related to the expulsion of a molecular layer. These alternating states explain the periodic modulation of the elastic properties: The maximum stiffness occurs for the commensurate case with the film thickness corresponding to an integer number of molecular layers [8][9][10][11][12][13][14][15][16][17][18][19].

A more complex picture, however, emerges regarding the modulation of the viscous forces exerted by a confined film on a vibrating AFM tip. Experimental findings of different groups disagree on whether or not the dissipative response shows an oscillatory behavior as the film thickness is reduced. Several studies on water [10][11][12], octamethyltetracyclosiloxane (OMTCS) [13,14], and dodecanol [17] reported a modulation of viscoelastic response, with a period corresponding to the molecular size; the maxima of stiffness and dissipation may coincide, or be shifted by half a period. These discrepancies have been related to a cross-over from liquid to solid behavior due to "jamming" [13], or could partly be due to the complex data analysis [15,16]. Very recent experiments confirm the oscillations of the damping coefficient: For a small AMF tip, its maxima coincide with minimum stiffness [18]; a more complex situation occurs for a larger tip with a rough surface [19].

In the present Letter we derive the viscoelastic linear response function. The free energy of n layers is calculated in harmonic approximation for the molecular interactions; local fluctuations are included in a grandcanonical approach with variable n. The dynamic stiffness and damping coefficient are obtained from Kubo behavior observed in experiments.

II. DYNAMIC RESPONSE

The thermodynamic potential Ω of a confined molecular liquid depends on the film thickness L. For thick films the solvation force

F = -dΩ/dL (1)
is constant and simply related to the hydrostatic pressure. Structuration effects in nanoscale films, however, lead to an oscillatory behavior of both the force F and the static stiffness

K 0 = d 2 Ω dL 2 . ( 2
)
Viscoelastic properties are obtained from the dynamic response function

φ(t -t ′ ) = {F (t), F (t ′ )} , (3)
which is given by the Poisson brackets of the fluctuating part of the solvation force [START_REF] Kubo | Statistical Physics II, Non-equilibrium Statistical Mechanics[END_REF]. The time dependent force F (t) is realized by vibrating an AFM tip at frequency ω and amplitude A, as illustrated in Fig. 1. The linear response approximation is justified if A is significantly smaller than the molecular size. The real parts of the Fourier transform φ(ω) = φ ′ (ω) + φ ′′ (ω) gives the dynamic contribution to the film stiffness

K(ω) = K 0 + φ ′ (ω), (4)
whereas the imaginary part is related to the damping coefficient,

γ = φ ′′ (ω) ω . (5)
Experimentally these quantities are extracted from the measured frequency shift and the attenuation of the vibration of the AFM tip.

In the following we evaluate the grand potential Ω in terms of the n-layer free energy in harmonic approximation. Since we are mainly interested in structuration effects, we consider the response function φ due to local fluctuations of the number of layers n.

III. HARMONIC APPROXIMATION

As suggested by the force experiments and confirmed by numerical work [7], the molecules in a confined film arrange in a regular array, as illustrated in Fig. 1. We start from a close packed structure [START_REF] Kittel | Introduction to Solid State Physics[END_REF] of n layers parallel to the confining walls; other symmetries may occur [START_REF] Schmidt | [END_REF][23][24][25] but would not change the essential aspects. The potential energy of the ordered domain then reads

E = i,j v(r ij ) + i,w v 0 (r iw ) , (6)
where the first sum runs over 12 nearest neighbors.

The second term accounts for the interaction of the upper and lower layers with the molecules of the confining walls. Expanding the molecular Lennard-Jones potentials to quadratic order in the distance, we have v(r) = 1 2 k (r -a) 2 . A similar expression with a constant k 0 arises for the interaction v 0 with the confining walls.

Elastic theory for an infinite system usually relies on the strain tensor ε αβ which is defined by deformation of the distance vector r ij with respect to that of the perfect lattice [START_REF] Kittel | Introduction to Solid State Physics[END_REF]. In the present case the strain arises from the boundary conditions imposed by the confining walls. The vertical mismatch with respect to the thickness of n unstrained layers reads as

ℓ n = L -nã, ã = 2 3 a, ( 7
)
where ã is the effective spacing in a cubic crystal and where we have discarded an irrelevant constant due to the boundary layers. The mismatch is related to the diagonal strain element through ℓ n = ε xx nã. Expanding (6) in powers of the vertical deviations xi from the equilibrium position of molecule i, we obtain

E = 1 2 i,j k ij (x i -xj ) 2 .
The k ij are finite for nearest neighbors only and take the values 2 3 k and 2 3 k 0 . As a consequence, the canonical partition function of a n-layer configuration reduces to a Gaussian form. Summing over the 3 neighbors in each adjacent layer and integrating over the xi , one obtains the canonical free energy

Ω n /N = k n 2 ℓ 2 n + 1 2 k B T ln ñ. ( 8
)
where N is the number of molecules per layer. The first term is the elastic energy, which depends on the square of the vertical mismatch ℓ n and the effective stiffness

k n = 2k ñ , ñ = (n -1) + 2k k 0 . (9)
The physical origin of k n becomes obvious when writing 2/k n = (n -1)/k + 2/k 0 : For n molecular layers, there are n -1 inter-layer couplings and the interaction of the outer layers with the solid boundaries. The case n = 1 describes a single layer coupled to both confining walls. For a few layers, the overall stiffness k n is determined by the softer of the springs. We use k 0 = k throughout this paper. Yet note that a strong coupling to the confining walls would tightly bind the outer molecular layers to the surfaces, whereas the opposite case (k 0 ≪ k) leads to a more complex alteration of the elastic response.

In Fig. 2 we plot Ω n as a function of the spacing L for several values of n. The elastic energy results in a parabola centered at ℓ n = 0, with a constant given by the second term in Eq. ( 8). If the vertical mismatch is much smaller than the molecular size, |ℓ n | ≪ ã, adding or subtracting a layer would cost a considerable amount of energy; thus the n-layer configuration is stable. A different situation arises if the values of the free energy corresponding to n and n + 1 are comparable. In terms of the mismatch this means ℓ n ≈ -ℓ n+1 . Then there are local fluctuations of small groups of particles between configurations n and n + 1, as illustrated in Fig. 2 for n = 3. The resulting grandcanonical potential

Ω = -N k B T ln n e -Ωn/NkBT (10)
is minimum for stable configurations and maximum at the intermediate states. This form is obtained under the assumption that the film is drained bead by bead through independent relaxation processes. If, on the contrary, a layer of N molecules was expelled as a whole, one would have to suppress both factors N occurring in (10).

The number of molecules in one ordered layer is given by the ratio N = A/A 1 of the domain size A and the area per molecule A 1 = ã2 3/4 [START_REF] Kittel | Introduction to Solid State Physics[END_REF]. In order to obtain an estimate for the former, we note that the spacing h between tip and solid support varies as

h = L + R - √ R 2 -r 2
, where r is the lateral distance from the tip axis. Existence of an ordered state requires that h -L be small as compared to the molcular size; neglecting a numerical factor we find the domain radius r 0 ≈ √ Ra. Thus the area A ≈ πRa and the corresponding number of molecules N ≈ πR/a increase linearly with the tip radius, and so does the free energy (8).

When varying the film thickness L, the grand potential Ω oscillates with period ã, and so do its derivatives. We start with the solvation force (2) exerted by the film on the walls. Taking the derivative, one finds F as the mean of the force in an n-layer configuration,

F/N = f n , f n = -k n ℓ n , (11)
where the average over the number of layers reads explicitly ... = n (...) e (Ω-Ωn)/NkBT .

(Here and in the following, we consider the force per molecule in one layer, F/N, which is identical to the disjoining pressure multiplied with the area per molecule.) From the curve Ω in Fig. 2 it is clear that its slope F oscillates as a function of L, in perfect agreement with various measurements on molecular liquids in the range n = 2...7. The envelope function decreases exponentially, F ∼ e -L/ξ ; for large stiffness the characteristic length ξ is proportional to k. The oscillatory solvation force has been reproduced by density-functional theory hardsphere systems [26]. The overall picture is confirmed by the accord of experimental and numerical studies of the drainage kinetics [28][29][30].

IV. THE ROLE OF FLUCTUATIONS

At maximum mismatch, two adjacent states of similar free energies Ω n ≈ Ω n+1 are equally probable, and the resulting fluctuations strongly affect the elastic response of the film. The underlying mechanism is best exposed in terms of the stiffness, or static susceptibility (3). From the grand potential one readily calculates

K 0 N = k n - δf 2 k B T . (12)
These two contributions are of rather different physical origin. The first one is given the mean stiffness k n ; it decreases monotonously with the inverse film thickness, as is clear from (9). The second one arises from the fluctuations of the solvation force,

δf 2 = f 2 n -f n 2 . ( 13
)
If n layers correspond to optimal filling ℓ n ≈ 0, adding or subtracting a layer requires a significant amount of energy; then the states n ± 1 are hardly occupied, δf 2 is small, and the stiffness is positive. If, on the other hand, the system size L lies between two values of optimal filling, fluctuations n ⇄ n + 1 are important and soften the elastic response of the confined liquid such that K 0 < 0. The curves labeled ∞ and 0 the upper panel of Fig. 3 correspond to N k n and K 0 , respectively; they show that the leading terms of k n and δf 2 /k B T vary linearly with the inverse film thickness and cancel each other. The oscillating behavior of K 0 is due to the fluctuations δf 2 ; the envelope decreases exponentially ∼ e -L/ξ with ξ ≈ 2ã. So far we have considered the static response of the confined film, where the change of the external control parameter L is much slower than the characteristic time τ of a configuration change. Dynamic effects are accounted for by the retarded response

φ(t -t ′ ) = N {f n (t), f n (t ′ )} , (14)
which is given by the Poisson brackets of the fluctuating part of the solvation force [START_REF] Kubo | Statistical Physics II, Non-equilibrium Statistical Mechanics[END_REF]. When assuming the fluctuations δf 2 to occur on a single time scale τ, one has for its Fourier transform

φ(ω) N = δf 2 k B T iωτ 1 + iωτ . ( 15
)
The complex susceptibility of a confined film is given by the sum K 0 + φ(ω). Its real part is the frequencydependent stiffness K = K 0 + φ ′ , which after rearranging the fluctuation parts reads as

K(ω) N = k n - δf 2 k B T 1 1 + ω 2 τ 2 . ( 16
)
Since at finite frequencies the fluctuation-driven part diminishes, the overall stiffness increases with ω, as illustrated in the upper panel of Fig. 3. The modulation as a function of the film thickness is due to the mean square force δf 2 . The quantity K(ω) describes the linear response to a periodic variation A cos(ωt) of the thickness. In the limit of a rapid perturbation (ωτ ≫ 1), where the duration of one cycle is much shorter than the mean relaxation time τ , the film does not have the time to adjust its configuration. Then the stiffness is given by the quasiinstantaneous response, that is, by the mean spring constant k n . In the case ω = 0 the external perturbation is slow, allowing the system to fully relax; the fluctuations reduce the stiffness to its static value (12). Note that the high-frequency response is always positive, whereas at low ω it may take either sign. The curves for small ωτ agree rather well with experiment [9][10][11][12][13][14].

The imaginary part of (15) gives the confinement-

induced damping coefficient γ = φ ′′ /ω, γ(ω) N = δf 2 k B T τ 1 + ω 2 τ 2 . ( 17
)
Assuming a constant relaxation time one finds that γ decreases with increasing frequency. The modulation as a function of L shown in the lower panel of Fig. 3, is due to the force fluctuations; damping minima occur at optimal filling where δf 2 is small, and the maxima in the intermediate states.

V. SHEAR RESPONSE FUNCTION present paper deals with motion perpendicular to the film. Formally our results apply equally well to shear motion, where the external driving force is parallel to the solid support. Experimentally this could be realized by horizontal small-amplitude vibrations of the AFM tip of Fig. 1. Here we address several aspects that may lead to important differences in the shear or draining response function.

The molecular ordering in a thin film is anisotropic. As shown in Fig. 1, the vertical confinement imposes a unique mismatch ℓ n ; the corresponding static stress is homogeneous throughout the confined volume. This is quite different in the parallel directions. If the surfaces were structured with the period corresponding to the molecular size, they would pin the first liquid layer and result in an elastic response function similar to (16). Commonly used materials such as silicates and metal coatings, however, provide rather flat surfaces; the molecules move freely along the solid boundaries and thus rapidly relax upon a static shear. In other words, there is no permanent lateral mismatch between the first and last layer. As a consequence, the static shear stiffness vanishes, and the response function consists of the dynamic suseptibility φ S only. Thus the stiffness reduces to the relaxator contribution

K S (ω) N = δf 2 S k B T ω 2 τ 2 S 1 + ω 2 τ 2 S , (18)
whereas the damping coefficient is similar to (17), albeit with parameters δf S and τ S . In particular this results in a linear frequency dependence of the ratio of elastic and viscous forces, φ ′ S /φ ′′ S = ωτ S , which is characteristic for the liquid state [27]. For the force fluctuations and the relaxation time one expects a behavior similar to that discussed above for the normal response function.

Many shear studies apply amplitudes well beyond the molecular size, and thus do not satisfy the limitations of the linear-response approximation used in the present work. Similarly, squeezing out successive molecular layers probes the rheology on large scales [28][29][30], which is not directly related to the small-amplitude vibrations. Like any complex many-body system, molecular films show a broad relaxation spectrum. Each experiment probes a particular motion, and the resulting time constants may be quite different.

VI. DISCUSSION

The oscillatory behavior of the solvation force is well established experimentally, and has been related to the the number of molecular layers in the film [START_REF] Israelachvili | Intermolecular and Surface Forces[END_REF][START_REF] Horn | [END_REF][3][4][5][6][7]. The role of fluctuations is illustrated in Fig. 2; the force F and stiffness K 0 are readily obtained from the local slope and curvature of the grand potential Ω, according to (1) and (2). The situation is less clear regarding structuration effects in the damping constant [8][9][10][11][12][13][14][15][16][17][18][19]. Depending on the system parameters, some studies reported oscillations as a function of the film thickness, others didn't. Moreover, experiments that did find oscillations, do not agree whether the damping maxima occur at the maxima or the minima of the stiffness. The present work was largely motivated by this latter question.

According to ( 16) and ( 17), elastic and viscous forces show opposite behavior: The maxima of the stiffness correspond to minima of the damping function and vice versa. These oscillations are related to fluctuations of the number of molecular layers; the fluctuation strength δf 2 is largest close at half-integer filling, and almost vanishes at optimal filling. Since fluctuations weaken the stiffness and enhance dissipation, these quantities show op-posite behavior as a function of the film thickness. These features agree with recent experimental findings on an OMTCS film confined by a AFM tip of radius R = 15 nm [18]. So far our discussion relied on the relaxation time approximation with constant τ . One should keep in mind that, due to the broad many-body relaxation spectrum, the quantity τ in general depends on frequency. More importantly, the relaxation time varies with the film thickness and is largest in the ordered state. This means that in Eq. ( 17) the factors δf 2 and τ show an opposite oscillatory behavior. If the modulation of the relaxation time dominates, maxima of the damping coefficient occur in the ordered state where L/ã is an integer. As an illustration, we compare in Fig. 4 the damping coefficient for constant τ with that obtained from an Arrhenius type relaxation time τ = τ 0 e E L /k B T , with a modulated activation energy E L = E 0 [1+cos(2πL/ã)]ã/L. For E 0 = 0 the damping maxima occur at intermediate states, whereas a high enough activation energy shifts them to optimal filling where L ≈ nã. Such a change of the relaxation behavior could be related to the size and surface roughness of the AFM tip as discussed in [19], or to "jamming" at higher drift velocity of tip [13].

We conclude with a remark on the formal equivalence of the response function (15) to Maxwell's relaxation model for viscoelastic fluids. The present work adds two features that are essential for confined films: First, the modulation of the prefactor δf 2 is at the orign of the oscillatory behavior of both elastic and dissipative response. Second, the response function contains, beyond the Maxwell relaxation term φ, a static contribution K 0 which completely changes the elastic behavior with respect to a Maxwell fluid. In view of these findings, the viscoelastic response of molecular layers arise from a subtle interplay of elastic energy and fluctuations.

Helpful discussions with A. Maali are gratefully acknowledged.

FIG. 1 :

 1 FIG.1: Schematic view of the confinement-induced ordering of a liquid film. The ordered domain is indicated by darker beads; its area A ≈ aR is given by the curvature radius R of the AFM tip and the molecular radius a.

FIG. 2 :

 2 FIG.2: Grand potential Ω and canonical excess free energies Ωn in units of kBT . Minima of Ω correspond to stable ordered configurations where the film thickness L is close to a multiple integer of ã. We put k = 0.1 N/m, ã = 0.8 nm, and N = 30. The intermediate states are highly strained and result in strong fluctuations of the local configuration.

3 FIG. 3 :

 33 FIG.3: Dynamic stiffness K(ω) and reduced damping coefficient γ/τ for different values of ωτ . Increasing ωτ augments the stiffness but reduces its modulation amplitude. The damping coefficient γ/τ is maximum in the static limit. The modulation is entirely due to the force fluctuations of Eq. (13). Note that damping maxima occur at minimum stiffness. Parameters are k = 0.06 N/m, ã = 0.8 nm, and N = 30.

FIG. 4 :

 4 FIG.4: Damping constant γ for a non-uniform relaxation time τ = τ0e E L /k B T , with a modulated activation energy EL = E0[1 + cos(2πL/ã)]ã/L. The curve for E0 = 0 is identical to those of Fig.3, with damping maxima in the intermediate states. A finite activation energy shifts the oscillations; for E 0 > k B T the maxima of γ occur in the ordered state, in phase with those of the stiffness. For the sake of clarity, we have shifted the vertical scale of the upper curves. Parameters as in Fig.3with τ0 = 10µsec and ωτ ≪ 1.